Quaternions to Euler, small numbers, NaNs, and the Decimal module

I’m trying to get better results for quaternion to euler conversions than I seem to get from the native Blender conversions. It’s turning out to be tricky. Apparently eulers are wily beasts. :slight_smile:

I have had the most success with the following code, but it fails significantly with the gimbal lock zones.

 
def quat_to_euler2(q1):
  #adapted from Maths - Conversion Quaternion to Euler - Martin Baker
  #http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/index.htm
  sqw = q1.w * q1.w
  sqx = q1.x * q1.x
  sqy = q1.y * q1.y
  sqz = q1.z * q1.z
  unit = sqx + sqy + sqz + sqw # if normalised is one, otherwise is correction factor
  test = q1.x * q1.y + q1.z * q1.w
 
  if (test > 0.499999*unit) and (test < 0.500001*unit):  # singularity at north pole
    heading = 2 * math.atan2(q1.x,q1.w)
    attitude = math.pi/2
    bank = 0    
  elif (test < -0.499999*unit) and (test > -0.500001*unit):  # singularity at south pole
    heading = -2 * math.atan2(q1.x,q1.w)
    attitude = -math.pi/2
    bank = 0
  else:     
    heading = math.atan2(2 * q1.y * q1.w-2 * q1.x * q1.z , sqx - sqy - sqz + sqw)
    attitude = math.asin(2 * test/unit)
    bank = math.atan2(2 * q1.x * q1.w-2 * q1.y * q1.z , -sqx + sqy - sqz + sqw)
  return r2d(bank),r2d(heading),r2d(attitude)
 

At the suggestion of someone on another message board, I adapted the following from another source. It returns somewhat worse results with the gimbal zones.

 
#This produces too many gimbal errors
def matToEuler_v2(mat): #mat can be 3X3 or 4X4
  #adapted from The Matrix and Quaternions FAQ
  #www.j3d.org/matrix_faq/matrfaq_latest.html
  angle_y = math.asin(mat[0][2])    #Calculate Y-axis angle
  C = math.cos(angle_y)
  angle_y = r2d(angle_y)
  if math.fabs(C) > 0.005:                #Gimball lock?
    x1 = mat[2][2]/C                #No, so get X-axis angle
    y1 = -mat[1][2]/C
    angle_x = r2d(math.atan2(y1,x1))
    x1 = mat[0][0]/C                #Get Z-axis angle
    y1 = -mat[0][1]/C
    angle_z = r2d(math.atan2(y1,x1))
  else:                             #Gimball lock has occurred
    angle_x = 0                     #Set X-axis angle to zero
    x1 = mat[1][1]                  #And calculate Z-axis angle
    y1 = mat[1][0]
    angle_z = r2d(math.atan2(y1,x1))
  return angle_x, angle_y, angle_z
#Radians conversion borrowed from xsi export, by "Elira"    
def r2d(r): #radians to degrees
  return round(r*180.0/math.pi,4) #Approx. 57.295828 before rounding
 
#This works
def quatToMat3(quat):
  w,x,y,z = quat
  m00 = 1 - 2*y**2 - 2*z**2
  m01 = 2*(x*y) - 2*(w*z)
  m02 = 2*(x*z) + 2*(w*y)
  m10 = 2*(x*y) + 2*(w*z)
  m11 = 1 - 2*x**2 - 2*z**2
  m12 = 2*(y*z) - 2*(w*x)
  m20 = 2*(x*z) - 2*(w*y)
  m21 = 2*(y*z) + 2*(w*x)
  m22 = 1 - 2*x**2 - 2*y**2
  mat = [[m00,m01,m02],[m10,m11,m12],[m20,m21,m22]]
  return mat

When I asked at this other board about tightening up the conditions for the gimbal/singularity handling, I received the following response:

The solution is probably dependant on the behaviour of asin and cos in the problem region.
To get some idea of the evaluation shortcuts and argument reduction that takes place inside these functions you can refer to the routines from fdlibm hosted at netlib.
I’d suggest confirming that your are getting values that produce what your expecting when passed to other the functions.
It may be that you’ll need to add more special case processing to deal with some NaNs.

I’ve looked into Nans, and Python’s built-in method for handling them seems to be the new Decimal module in Python 2.4. But this isn’t available with Blender’s normal installation. So I thought I should ask if there is any known alternate way to guarantee good results with very small numbers and such. Or if the Decimal module could be considered for inclusion as part of the Blender distribution, perhaps. It seems like it might be useful.

Does anyone have any thoughts or suggestions on any of this?

More. I tried adapting all of this from the Blender C source, to see if I could achieve more controlled results from it than BPy normally gives me. But it didn’t really work. I’m wondering whether I have the value of FLT_EPSILON wrong, or if the handling of such values differs significantly between Python and C++.

 
#I'm assuming that FLT_EPSILON == FLT_EPSILON10 divided by ten.  FLT_EPSILON is never defined in the Blender source.
#Maybe the right answer needs to be derived using cmath?
FLT_EPSILON10 = 1.19209290e-06
print FLT_EPSILON10/10
FLT_EPSILON = FLT_EPSILON10/10
#FLT_EPSILON = 1.1920929e-007
print 16.0*FLT_EPSILON
 
#This makes a horrible mess... it may need the normalization routine from the source, or FLT_EPSILON may be wrong....
def Mat3ToEul(mat): #adapted from arithb.c in Blender source
  FLT_EPSILON = 1.1920929e-007
  eul = [0,0,0]
  cy = math.sqrt((mat[0][0]*mat[0][0]) + (mat[0][1]*mat[0][1]))
  
  if cy > 16.0*FLT_EPSILON:
    eul1 = [0,0,0]
    eul2 = [0,0,0]
    
    eul1[0] = r2d(math.atan2(mat[1][2],mat[2][2]))
    eul1[1] = r2d(math.atan2(-mat[0][2],cy))
    eul1[2] = r2d(math.atan2(mat[0][1],mat[0][0]))
    eul2[0] = r2d(math.atan2(-mat[1][2],-mat[2][2]))
    eul2[1] = r2d(math.atan2(-mat[0][2],cy))
    eul2[2] = r2d(math.atan2(-mat[0][1],-mat[0][0]))
    #if abs(eul1[0])+abs(eul1[1])+abs(eul1[2]) > abs(eul2[0])+abs(eul2[1])+abs(eul2[2]):
    if math.fabs(eul1[0])+math.fabs(eul1[1])+math.fabs(eul1[2]) > math.fabs(eul2[0])+math.fabs(eul2[1])+math.fabs(eul2[2]):
      eul = eul1
    else: eul = eul2
  else:
    eul[0] = r2d(math.atan2(-mat[2][1],mat[1][1]))
    eul[1] = r2d(math.atan2(-mat[0][2],cy))
    eul[2] = 0.0
  return eul
#This gets most of it right, but is worse than average, producing too many wholly incorrect rotations.
def Mat3ToEuln(tmat): #adapted from arithb.c in Blender source
  FLT_EPSILON = 1.1920929e-007
  sin1 = -tmat[2][0]
  cos1 = math.sqrt(1 - sin1*sin1)
  eul = [0,0,0]
  if math.fabs(cos1) > FLT_EPSILON:
    sin2 = tmat[2][1]/cos1
    cos2 = tmat[2][2]/cos1
    sin3 = tmat[1][0]/cos1
    cos3 = tmat[0][0]/cos1
  else:
    sin2 = -tmat[1][2]
    cos2 = tmat[1][1]
    sin3 = 0.0
    cos3 = 1.0
  eul[0] = r2d(math.atan2(sin3,cos3))
  eul[1] = r2d(math.atan2(sin1,cos1))
  eul[2] = r2d(math.atan2(sin2,cos2))
  
  return eul
#untested
def QuatToMat3(q): #adapted from arithb.c in Blender source
  M_SQRT2 = 1.41421356237309504880
  m = [[1,0,0],[0,1,0],[0,0,1]]
  
  q0= M_SQRT2 * q[0]
  q1= M_SQRT2 * q[1]
  q2= M_SQRT2 * q[2]
  q3= M_SQRT2 * q[3]
  qda= q0*q1
  qdb= q0*q2
  qdc= q0*q3
  qaa= q1*q1
  qab= q1*q2
  qac= q1*q3
  qbb= q2*q2
  qbc= q2*q3
  qcc= q3*q3
  m[0][0]= (1.0-qbb-qcc)
  m[0][1]= (qdc+qab)
  m[0][2]= (-qdb+qac)
  m[1][0]= (-qdc+qab)
  m[1][1]= (1.0-qaa-qcc)
  m[1][2]= (qda+qbc)
  m[2][0]= (qdb+qac)
  m[2][1]= (-qda+qbc)
  m[2][2]= (1.0-qaa-qbb)
  
  return m

If anyone can shed any light on this… or point out where I’ve erred (I’m sure I probably have), please do…

Sorry for the diversion, but exactly why are Euler angles needed? I’ve had very good luck with interconversions using matrix notation.

myEuler.unique() should generate a unique Euler rotation that avoilds gimbal lock. Unless your rotations are fairly simply It may behoove you to use quaternions instead.

Sorry for the diversion, but exactly why are Euler angles needed? I’ve had very good luck with interconversions using matrix notation.

I guess I should have clarified the purpose of all this. I’m trying to convert Blender armature rotations for an export format. I have to export euler angles for that format, and the normal toEuler doesn’t seem to give the results needed. The quat to euler code I posted seems better for this purpose, spcifically with the gimbal zones, but it still needs a bit of improvement.

myEuler.unique() should generate a unique Euler rotation that avoilds gimbal lock. Unless your rotations are fairly simply It may behoove you to use quaternions instead.

I rely on euler.unique a great deal in my overall project, but, sadly, it actually worsens things for the purpose of this export process. The frustrating thing is that Blender’s conversions, unique or not, are perfectly fine within Blender. They just fall short when the rotations are to be applied in the destination application.

Blenders conversion was improved recently, might be worth seeing if that meets your needs.

Blenders conversion was improved recently, might be worth seeing if that meets your needs.

Cool. :slight_smile: How recently? I’ve tested with 242 initial release and with the 9/7/06 CVS build from the JMS site.

I’m downloading the 9/15 build from the same site, now.

Crossing my fingers. :slight_smile:

1 or 2 weeks ago? - Im not even sure that python API makes use of this area of the code.

Im not even sure that python API makes use of this area of the code.

Do you mean the code which you say was changed, or the code I was adapting from the source, up above? :confused:

I tested with the 9/15/06 CVS build, and got the same problematic results I was getting with previous tests of toQuat. :frowning: So I’m still stuck with my initial inquiries.

Cambo, this matter might be of some interest to you, actually. I’m trying to export Poser poses. Wouldn’t these same problems obtain with BVH export, however? That’s one of your projects, right?

Thanks for translating the math to Python. It solved an old problem for me. I have adopted your first code example to create pydrivers for my model to use instead of the standard X and Z rotation drivers which drive shape keys to correct deformation of a shoulder.

I can verify Blender’s Euler conversions are sometimes wrong. If I twist (or roll) a bone, say 45° on its Y axis, the X and Z coordinate system changes. If I rotate constraining on the bone’s local X axis, the 3D viewport’s number input shows correct values. If I rotate on the local Z axis, I get the correct values using your code, but the values I get from the GUI and from toEuler() make no sense.

The shoulder of my model now deforms the way it should. Even when the arm twists, because I can make keys for the isolated Y rotation.

could you post a sample blend to the bug tracker that shows the numbers being wrong in Blender?

LetterRip

Glad that helped, tolobán. It looks to me like the gimbal conditions aren’t catching anything in the code I posted above, however… I monkeyed around with the conditions, trying to refine them, but I think I actually broke them. You might want to check that in your own implementation. In the meantime, I’m still fighting with this…

could you post a sample blend to the bug tracker that shows the numbers being wrong in Blender?

Are you directing that at me, Algorith, or at tolobán? I’ll see if I can put something together which illustrates the dilemma. The trouble is that the converted rotations look just fine in Blender, using Blender’s native conversions. They simply aren’t what is needed on the receiving end, unfortunately. Apparently one can resolve toEuler conversions in various different ways, ending up with variations on the same rotation. Poser seems to need one specific variation. And I have no idea how to produce that variant… :frowning: :confused:

Basil,

is the rotation order that you are applying them correct? Eulers are order dependent (ie rotating x, then y, then z, is different from rotating z,y,x; z,x,y; x,z,y; y,z,x; or y,x,z). I’d try rotations in all 6 orders and see if one of them is right. Also I think you might need to do make changes if the xyz orientation is different in the other program. (Ie Blender is Z+ is the top of the screen, Y+ is into the screen, and X+ is to the right of the screen. There are again a bunch of different combinations that it could be).

LetterRip

LetterRip, I am actually parsing the target rotation orders for Poser from a Poser file, so the order issues are controlled. I’ve tried the various combinations of orders for both import and export. Varying orders upon import into Blender seems unnecessary - it actually worsens things. Varying orders upon export from Blender makes some difference in Poser (incorrect orders usually create more problems), but rotation order seems to be only one aspect of the matter. For instance, I test in Poser with two figures which are derived from the same .cr2 source. They have the same rotation order and joint centers, but they respond differently to the same poses. So something odd is still happening… The problems, when importing into Poser, occur with rotations which should fit the gimbal conditions: one or more axes rotated at or near 90 degrees. But none of the gimbal handling tweaks have really resolved this. I think this problem needs someone with more math smarts than I have… :slight_smile: I’m in over my head.

LetterRip, your comment on the world axis order might present a new, er, angle on the problem. Poser uses y up, compared to Blender’s z up. I don’t think that is the direct cause of any of the problems - too much already comes out correct for something so broad to be the trouble, I would think. But it may be reflected somehow in the quat to euler code I adapted from the Martin Baker site (the code which tolobán has used).

I have adapted the code from the following link:

http://www.euclideanspace.com/maths/geometry/rotations/conversions/quaternionToEuler/index.htm

I didn’t really change much, aside from altering his singularity testing conditions, but I do return the elements in a different order than his code intends. I change his heading, attitude, bank to bank, heading attitude. I am now wondering whether this represents anything equivalent to a world axis rotation. :confused: It looks like the standard he’s applying expects the z-y ordering which Belnder submits. But the fact that I’ve had to alter the order of his reults makes me wonder whether his code could be better adapted by changing the elements in the actual calculations. I’ve looked at his proofs, and the matrix math he uses is over my head. I’m wondering whether a re-ordering of the matrix elements he uses in deriving his conversion might give a better result for this purpose.

Does anyone with a better grasp of math have any ideas?

BF,

After “digging around in the basement”, do you have any insight on a user level on what the Y values represent in the IPO editor for bone rotations ? and/or what the “W” curve represents.

Mike

After “digging around in the basement”, do you have any insight on a user level on what the Y values represent in the IPO editor for bone rotations ? and/or what the “W” curve represents.

Mike - probably not, since I’m not quite sure what you’re asking. “Y” and “W” in the Ipo curve constants are, of course, components of the quaternion representing the bone’s rotation. But you presumably know that, so I assume you’re asking something else. Can you clarify your question? :confused:
(Like I keep saying - ol’ Bas is a little slow, sometimes.) :frowning:

 
icu_quatw = ipo[Ipo.PO_QUATW] #25
icu_quatx = ipo[Ipo.PO_QUATX] #26
icu_quaty = ipo[Ipo.PO_QUATY] #27
icu_quatz = ipo[Ipo.PO_QUATZ] #28

Mike_S:

Quaternions are basically an angle around an axis. The ‘W’ component represents the angle of rotation, around the axis going through the origin and the point defined by ‘x’ ‘y’ ‘z’ components. I’m not too familiar with them yet myself, so I might be wrong on some aspect of this.

Aligorith

Well, the quaternions have a real aspect and an imaginary aspect. I can never remember whether the W is the real or the imaginary part…
I kind of feel like I’ve misunderstood the whole question…

<H2>Basic Quaternion Arithmetic

Quaternions form what is called a non-commutative division ring. This space can be visualised in a similar manner to complex numbers. Just as the complex numbers can be represented by a two-dimensional co-ordinate on the complex plane, quaternions can be represented by a co-ordinate in 4-D space. This is written:
q = w + xi + yj + zkThe numbers i, j and k are related to the real numbers and each other by the following equalities:
i2 = j2 = k2 = -1
ij = kQuaternions can be viewed as a vector space over the real numbers. In this way they are given the notation q = (w, x, y, z), with w, x, y and z defined as above. When adding or subtracting quaternions you perform the same operation as for vectors and matrices: add and subtract component-wise.
As with many other vector spaces, the quaternions come equipped with a norm N, which is function (usually regarded as a “length”) that returns a single real number given a quaternion. It is usually defined as:

</H2>
http://www.sjbrown.co.uk/