Converting HPB to Euler - exactly
-
THE POST BELOW IS MORE THAN 5 YEARS OLD. RELATED SUPPORT INFORMATION MIGHT BE OUTDATED OR DEPRECATED
On 20/07/2007 at 15:22, xxxxxxxx wrote:
User Information:
Cinema 4D Version: R8.2-R10
Platform: Windows ; Mac ; Mac OSX ;
Language(s) : C++ ;---------
It took literally years (!) but I have found a stable and consistent method for converting Cinema 4D HPB or matrix/quaternion equivalent into exact ordered Euler angles - not equivalents (like Eberly, Shoemake, et al).First, I'm going to pay homage to Greg Abbas for the only viable solution that I've ever encountered ever anywhere. The article that led to the solution is at this link:
http://www.allyourpixel.com/post/well-behaved-euler-angles/
Now, this is not a direct solution. It is a heuristic solution. That is, it expects smooth transitions (deltas/changes) in rotation to allow the conversion to Euler angles by way of the previous Euler angles as comparative equivalents. So, this won't work if you just type in any rotation and then get the matrix and convert - it requires that there is a previous state (x-dx, y-dx, z-dx) where dx,dy,dz are relatively small angle deltas. This makes it very good for user interaction with the mouse - but not good for random values entered into the HPB coords (etc.).
This is exactly the idea that I had - but without the mathematical rigor needed to develop the heuristic algorithm.
The code considers rotation order (XYZ, XZY, YXZ, YZX, ZXY, ZYX) and is taking a previous representation in a left-hand (-z) and applying it to results in a right-hand systems (note the z negation at the end). The equivalent angles are extracted from the matrix using David Eberly's routines first. Then they are made 'well-behaved' using Greg Abbas' algorithm. This is not a general solution for rotations greater than 360d. I leave it up to the user to generalize for that case.
Note that the general solution must be considered both for +/-360d integrals and for the 180d identity similarly.static CHAR* RotOrdStr[7] = { "", "XYZ", "XZY", "YXZ", "YZX", "ZXY", "ZYX" }; // IPPDial.HPBToPoser - Convert HPB to Poser ordered-Euler angles //*---------------------------------------------------------------------------* void IPPDial::HPBToPoser(LONG rotorder, Matrix mat, Vector prot) //*---------------------------------------------------------------------------* { Vector v1; Vector v2; Real sx = 1.0f, sy = 1.0f, sz = 1.0f; Real x, y, z; Real ox = prot.x, oy = prot.y, oz = prot.z; // Take into account Rotation Order and extract Euler angles switch (rotorder) { case ROTORDER_XYZ: v1 = mat.v1; v2 = mat.v2; y = asin(v1.z); if (y < pi05) { if (y > -pi05) { x = Angle(-v2.z, mat.v3.z); z = Angle(-v1.y, v1.x); } else { x = -Angle(v2.x, v2.y); z = 0.0f; } } else { x = Angle(v2.x, v2.y); z = 0.0f; } sy = -sy; break; case ROTORDER_XZY: v1 = mat.v1; v2 = mat.v3; z = asin(-v1.y); if (z < pi05) { if (z > -pi05) { x = Angle(v2.y, mat.v2.y); y = Angle(v1.z, v1.x); } else { x = -Angle(-v2.x, v2.z); y = 0.0f; } } else { x = Angle(-v2.x, v2.z); y = 0.0f; } sz = -sz; break; case ROTORDER_YXZ: v1 = mat.v1; v2 = mat.v2; x = asin(-v2.z); if (x < pi05) { if (x > -pi05) { y = Angle(v1.z, mat.v3.z); z = Angle(v2.x, v2.y); } else { y = -Angle(-v1.y, v1.x); z = 0.0f; } } else { y = Angle(-v1.y, v1.x); z = 0.0f; } sx = -sx; break; case ROTORDER_YZX: v1 = mat.v3; v2 = mat.v2; z = asin(v2.x); if (z < pi05) { if (z > -pi05) { y = Angle(-v1.x, mat.v1.x); x = Angle(-v2.z, v2.y); } else { y = -Angle(v1.y, v1.z); x = 0.0f; } } else { y = Angle(v1.y, v1.z); x = 0.0f; } sz = -sz; break; case ROTORDER_ZXY: v1 = mat.v1; v2 = mat.v3; x = asin(v2.y); if (x < pi05) { if (x > -pi05) { z = Angle(-v1.y, mat.v2.y); y = Angle(-v2.x, v2.z); } else { z = -Angle(v1.z, v1.x); y = 0.0f; } } else { z = Angle(v1.z, v1.x); y = 0.0f; } sx = -sx; break; case ROTORDER_ZYX: v1 = mat.v1; v2 = mat.v3; y = asin(-v2.x); if (y < pi05) { if (y > -pi05) { z = Angle(mat.v2.x, v1.x); x = Angle(v2.y, v2.z); } else { z = -Angle(-v1.y, v1.z); x = 0.0f; } } else { z = Angle(-v1.y, v1.z); x = 0.0f; } sy = -sy; break; } Vector d1, d2; Real dc, dp, dn; Real adc, adp, adn; // Find minimum relation in full revolutions (factors of 2*pi or 360d) // r10 = r +/- (n*pi2) // r11 = r +/- (n*pi2) // r12 = r +/- (n*pi2) // - X dc = x; dp = x + pi2; dn = x - pi2; adc = Abs(x - ox); adp = Abs(dp - ox); adn = Abs(dn - ox); d1.x = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // - Y dc = y; dp = y + pi2; dn = y - pi2; adc = Abs(y - oy); adp = Abs(dp - oy); adn = Abs(dn - oy); d1.y = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // - Z dc = z; dp = z + pi2; dn = z - pi2; adc = Abs(z - oz); adp = Abs(dp - oz); adn = Abs(dn - oz); d1.z = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // Make relative to hemi revolution (factor of pi or 180d) // r20 = pi + r10 +/- (n*pi2) // r21 = pi - r11 +/- (n*pi2) // r22 = pi + r12 +/- (n*pi2) // sx,sy,sz handle subtraction on second rotation (pi - r + (n*pi2)) // - X dc = pi + (sx*x); dp = pi + (sx*x) + pi2; dn = pi + (sx*x) - pi2; adc = Abs(dc - ox); adp = Abs(dp - ox); adn = Abs(dn - ox); d2.x = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // - Y dc = pi + (sy*y); dp = pi + (sy*y) + pi2; dn = pi + (sy*y) - pi2; adc = Abs(dc - oy); adp = Abs(dp - oy); adn = Abs(dn - oy); d2.y = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // - Z dc = pi + (sz*z); dp = pi + (sz*z) + pi2; dn = pi + (sz*z) - pi2; adc = Abs(dc - oz); adp = Abs(dp - oz); adn = Abs(dn - oz); d2.z = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // Find closest to reference triplet (prot) Real d1a, d2a; d1a = Abs(d1.x - ox); d2a = Abs(d2.x - ox); x = (d1a < d2a) ? d1.x : d2.x; d1a = Abs(d1.y - oy); d2a = Abs(d2.y - oy); y = (d1a < d2a) ? d1.y : d2.y; d1a = Abs(d1.z - oz); d2a = Abs(d2.z - oz); z = -((d1a < d2a) ? d1.z : d2.z); GePrint("Rotation ("+String(RotOrdStr[rotorder])+") = "+VectorToString(Vector(ToDeg(x),ToDeg(y),ToDeg(z)))+" (xyz)"); }Notes: Angle() is a 'smart' atan2() wrapper:
// Return Real angle (in radians) whose tangent = y/x //*---------------------------------------------------------------------------* Real Angle(Real y, Real x) //*---------------------------------------------------------------------------* { if (Abs(x) < 0.000001f) { if (Abs(y) < 0.000001f) return 0.0f; if (y > 0.0f) return (RAD_90); return (RAD_270); } return atan2(y,x); }And ToDeg() should be obvious. It is an inline function to convert radians to degrees:
inline Real ToDeg(Real r) { return ((r)*57.295779513082320876798154814105f); } -
THE POST BELOW IS MORE THAN 5 YEARS OLD. RELATED SUPPORT INFORMATION MIGHT BE OUTDATED OR DEPRECATED
On 23/07/2007 at 21:02, xxxxxxxx wrote:
Now I'm really peeved! These tests were done with a pristine rotation matrix concatentated in the form:
Vector prevRot = bc- >GetVector(IPP_ROTATION);
// Any of x, y, z were changed on my plugin's dials
Matrix rot = MatrixRotZ(z) * MatrixRotY(y) * MatrixRotX(x);
// Extract the exact Euler into x, y, z
HPBToPoser(order, rot, prevRot);
bc->SetVector(Vector(x, y, z));
obj->SetRot(MatrixToHPB(rot));
obj->Message(MSG_UPDATE);And off it went and converted the matrix back verbatim to the original x,y,z rotation values every time under many, many, many tests. This is the same matrix (as you can plainly see) applied to the object.
I figured there would be no difference between this matrix and the one from HPBToMatrix(obj->GetRot()) or from obj->GetMl() (considering all influences like orientation, scale, translation which are removed). But, NO! Cinema 4D seems to be employing its anti-gimbal-lock magick on the rotation values here in a way that defies even this rather robust extraction process (and it is very ROBUST in its limited scope of application - heuristically).
The results are definitely more reliable - depending upon how big the HPB rotations are. But at some stage (hmmm, lemme guess, somewhere around where at least two rotations are near +/-90d and so on) the conversion flips out. This could ONLY be because your matrices are reflecting the ole swap-the-axes trick to avoid gimbal lock and it is causing the incorrect angles to be associated during the extraction (say, because the H axis turned into the B axis and vice versa).
If I can't get a straight answer on how to avoid this (say, maybe, possibly, GetOptimalAngle(), maybe?) I throw up my hands, surrender, and move onto another software's plugin development.
I'm practically an expert on rotation extraction from matrices. I've talked to DAVID EBERLY (you might have seen his books, papers, website). I used the only available source I know, chanced upon in one of my many varied searches, for exact Euler rotation extraction and yet Cinema 4D HPB (nice and proprietary) still confounds even this vast armada of weaponry!
If I can't have information on this, interPoser Pro dies here.
Bye,
-
THE POST BELOW IS MORE THAN 5 YEARS OLD. RELATED SUPPORT INFORMATION MIGHT BE OUTDATED OR DEPRECATED
On 24/07/2007 at 07:33, xxxxxxxx wrote:
I try to get an answer on this.
cheers,
Matthias -
THE POST BELOW IS MORE THAN 5 YEARS OLD. RELATED SUPPORT INFORMATION MIGHT BE OUTDATED OR DEPRECATED
On 24/07/2007 at 07:50, xxxxxxxx wrote:
CINEMA uses the same code that is in the API for HPBToMatrix and MatrixToHPB (the source is in c4d_tools.cpp). Also GetOptimalAngle is in c4d_tools.cpp. When you call SetRot CINEMA uses HPBToMatrix to store the rotation and its matrix. When you call SetMl then it uses MatrixToHPB with GetOptimalAngle of the current HPB.
-
THE POST BELOW IS MORE THAN 5 YEARS OLD. RELATED SUPPORT INFORMATION MIGHT BE OUTDATED OR DEPRECATED
On 24/07/2007 at 11:11, xxxxxxxx wrote:
At first I thought the flipping out might be due to the IK (three bones being generally rotated together might be too much for the conversion to do in a timely fashion), but I then tested by rotating one bone using the Rotate tool and watched it flip out eventually as well.
Interestingly, it takes a while for the flip out to occur most of the time. Almost as if inaccuracies are mounting. They are assuredly not mounting in that code above - it is simply addition/substraction and the atan2 for extraction (and always worked faithfully with the concatenated rotation matrix no matter how long I twiddled using my plugin's rotation tool). And I've since converted the method to be using doubles (i.e.: LMatrix, LVector, LReal) to maintain a higher degree of accuracy during the conversion process and only convert to Real x,y,z at the very end.
I'm certain that the orientation is the only other rotation consideration to be removed in the local rotation - it may be part of the problem (said with much doubt) but that would mean that the stored orientation matrix is incorrect and it is definitely not as it is used in all rotations (including mesh deformations).
The code that imparts the anti-gimbal lock might be right there in GetHPB(), but how would that help? I'm relying upon the local rotations in matrix form to extract Euler angles in particular orders. I tried the method above without the Eberly extraction code and just using the HPB angles returned by obj- >GetRot() - again minus the orientation. That was disastrous! The body parts were already twitching as soon as IK was enabled.
The problem looks like 'dueling rotation orders'. That is, I set the rotations from my Eulers using the matrix concatenation in a particular order - this always works in that direction (Euler sliders to C4D rotation matrix). But now I'm dealing with rotations set on Cinema 4D HPB and creating its own matrix with no regard for my Euler rotations (HPB/Ml to C4D rotation matrix to Euler sliders). Extracting the rotations from there in the expected rotation order mostly works (mostly isn't good enough) as long as correspondences in rotation value are maintained. But it seems that once Cinema 4D does the anti-gimbal lock, this changes the related rotation order in the matrix (beyond my control) and the extraction results go haywire. What I'm saying is that the ordered rotation matrix concatentation is always stable and the C4D rotation matrix is always stable but the two diverge if only the latter is considered in not only rotations but seemingly in rotation order. And I say this because these changes must change the Eberly (and Shoemake) extraction results enough that the extracted rotations cannot be reconciled with the previous rotations for the heuristic determination.
Everyone seems to be able to offer 'answers', but has any of you actually successfully done this? I don't just write the code and complain - I endlessly test the code and watch it falter and fail inextricably and research over and over. There is an IEEE paper done by Robotic engineers trying to do motion capture from a live subject and impart it onto a Euler-based Robotic simulation (on the computer). Their similar stance is that Shoemake and other matrix extraction routines can't return exactly representable Euler angles to the original ones and had instead to go directly Euler to Euler and forgo the intermediate matrix. I can only note that these guys are problably Engineering PhDs with more maths under their belts than all of us combined (!). And I don't have this option unfortunately (either PhD or Euler-Euler).
The biggest point of this exercise is to avoid rolling my own IK system (after having rolled my own bone+deformation, morph, master-slave, tool, and other systems for the plugin suite) and instead use C4D IK which would make everyone (esp. me) very happy. As noted, it works rather well but those problems can't be excused in a commercial plugin. I can sit here and print HPB and Euler rotations and see where and when the problems begin and try it for various HPB angles and Euler orders and attempt to make a generalization on the process. Wouldn't it be simpler if I just had a direct way to know when this might happen. I don't see anything in GetHPB() standing out and saying 'Hey, this is where we're screwing up your code - yeah, right here- >."
-
THE POST BELOW IS MORE THAN 5 YEARS OLD. RELATED SUPPORT INFORMATION MIGHT BE OUTDATED OR DEPRECATED
On 24/07/2007 at 11:56, xxxxxxxx wrote:
Here's an alternative idea. Maybe someone can comment on its efficacy?
Instead of all of this mumbo jumbo which doesn't seem to work with Cinema 4D HPB, how about two direction vectors built from the two rotations, VectorToHPB() both, and then GetOptimalAngle() or just go straight for VectorAngle()? The good thing here is that I could equivocate both the previous and new HPB rotations to get rotation deltas using this and apply them to the previous Eulers. The idea is that the direction vectors are at least non-denominational - with respect to angles. Does VectorAngle() always return the direct angle between vectors or migh it return the opposite?
Am I desperate enough yet?
-
THE POST BELOW IS MORE THAN 5 YEARS OLD. RELATED SUPPORT INFORMATION MIGHT BE OUTDATED OR DEPRECATED
On 24/07/2007 at 13:31, xxxxxxxx wrote:
Okay. The direction vector method has the same problem - except that this method seems to neglect twist (Bank) because the two vectors can only provide axial rotation and not rotation along the vector length (I suppose).
If either of you would like to see this for yourself, I'd be glad to provide the current plugin build (Windows R10 but R9.x could be provided) and a C4D scene with a figure being tested. It prints the Euler 'Old' rotation vector 'rotV', the delta from VectorAngle() 'oDV', and the 'New' rotation vector 'nDV = rotV + oDV'. Whenever the bone is rotated past a certain extent, the delta values jump to like 90-delta or 180-delta and the deformations show the results (the bones don't since they are being rotated directly by the HPB in this case). The mesh deformations (and other areas) on the other hand require the correct ordered Euler rotations for deformation weighting - and thus that need being sought so vigorously.
Now, I could simply check for deltas beyond some small value but this sounds naive (and probably is). It is also interesting that the Vector retrieved from VectorAngle() is always three matched values no matter which axis is being rotated (or how many!).
-
THE POST BELOW IS MORE THAN 5 YEARS OLD. RELATED SUPPORT INFORMATION MIGHT BE OUTDATED OR DEPRECATED
On 25/07/2007 at 16:26, xxxxxxxx wrote:
An update to the code. One thing that was missed in transfer from the test area to the live area (IK) was the second Abbas set uses the first Abbas results! I've since broken up this method by rotation order so that certain factors are direct (no more sign variables used on the second Abbas set). Still not 100% stable with HPB anti-gimbal lock, but it's very close with some checks at the end for x,y,z > 360d (while(x > pi2) x -= pi2 type thing) and for differences between x,y,z and ox,oy,oz > 45d (x = ox) to avoid some of the instabilities.
static CHAR* RotOrdStr[7] = { "", "XYZ", "XZY", "YXZ", "YZX", "ZXY", "ZYX" }; // IPPDial.HPBToPoser - Convert HPB to Poser ordered-Euler angles //*---------------------------------------------------------------------------* void IPPDial::HPBToPoser(LONG rotorder, Matrix mat, Vector prot) //*---------------------------------------------------------------------------* { Vector v1; Vector v2; Real sx = 1.0f, sy = 1.0f, sz = 1.0f; Real x, y, z; Real ox = prot.x, oy = prot.y, oz = prot.z; // Take into account Rotation Order and extract Euler angles switch (rotorder) { case ROTORDER_XYZ: v1 = mat.v1; v2 = mat.v2; y = asin(v1.z); if (y < pi05) { if (y > -pi05) { x = Angle(-v2.z, mat.v3.z); z = Angle(-v1.y, v1.x); } else { x = -Angle(v2.x, v2.y); z = 0.0f; } } else { x = Angle(v2.x, v2.y); z = 0.0f; } sy = -sy; break; case ROTORDER_XZY: v1 = mat.v1; v2 = mat.v3; z = asin(-v1.y); if (z < pi05) { if (z > -pi05) { x = Angle(v2.y, mat.v2.y); y = Angle(v1.z, v1.x); } else { x = -Angle(-v2.x, v2.z); y = 0.0f; } } else { x = Angle(-v2.x, v2.z); y = 0.0f; } sz = -sz; break; case ROTORDER_YXZ: v1 = mat.v1; v2 = mat.v2; x = asin(-v2.z); if (x < pi05) { if (x > -pi05) { y = Angle(v1.z, mat.v3.z); z = Angle(v2.x, v2.y); } else { y = -Angle(-v1.y, v1.x); z = 0.0f; } } else { y = Angle(-v1.y, v1.x); z = 0.0f; } sx = -sx; break; case ROTORDER_YZX: v1 = mat.v3; v2 = mat.v2; z = asin(v2.x); if (z < pi05) { if (z > -pi05) { y = Angle(-v1.x, mat.v1.x); x = Angle(-v2.z, v2.y); } else { y = -Angle(v1.y, v1.z); x = 0.0f; } } else { y = Angle(v1.y, v1.z); x = 0.0f; } sz = -sz; break; case ROTORDER_ZXY: v1 = mat.v1; v2 = mat.v3; x = asin(v2.y); if (x < pi05) { if (x > -pi05) { z = Angle(-v1.y, mat.v2.y); y = Angle(-v2.x, v2.z); } else { z = -Angle(v1.z, v1.x); y = 0.0f; } } else { z = Angle(v1.z, v1.x); y = 0.0f; } sx = -sx; break; case ROTORDER_ZYX: v1 = mat.v1; v2 = mat.v3; y = asin(-v2.x); if (y < pi05) { if (y > -pi05) { z = Angle(mat.v2.x, v1.x); x = Angle(v2.y, v2.z); } else { z = -Angle(-v1.y, v1.z); x = 0.0f; } } else { z = Angle(-v1.y, v1.z); x = 0.0f; } sy = -sy; break; } Real d1x, d1y, d1z, d2x, d2y, d2z; Real dc, dp, dn; Real adc, adp, adn; // Find minimum relation in full revolutions (factors of 2*pi or 360d) // r10 = r +/- (n*pi2) // r11 = r +/- (n*pi2) // r12 = r +/- (n*pi2) // - X dp = x + pi2; dn = x - pi2; adc = Abs(x - ox); adp = Abs(dp - ox); adn = Abs(dn - ox); d1x = (adc < adp) ? ((adc < adn) ? x : dn) : ((adp < adn) ? dp : dn); // - Y dp = y + pi2; dn = y - pi2; adc = Abs(y - oy); adp = Abs(dp - oy); adn = Abs(dn - oy); d1y = (adc < adp) ? ((adc < adn) ? y : dn) : ((adp < adn) ? dp : dn); // - Z dp = z + pi2; dn = z - pi2; adc = Abs(z - oz); adp = Abs(dp - oz); adn = Abs(dn - oz); d1z = (adc < adp) ? ((adc < adn) ? z : dn) : ((adp < adn) ? dp : dn); // Make relative to hemi revolution (factor of pi or 180d) // r20 = pi + r10 +/- (n*pi2) // r21 = pi - r11 +/- (n*pi2) // r22 = pi + r12 +/- (n*pi2) // sx,sy,sz handle subtraction on second rotation (pi - r + (n*pi2)) // - X dc = pi + (sx*d1x); dp = pi + (sx*d1x) + pi2; dn = pi + (sx*d1x) - pi2; adc = Abs(dc - ox); adp = Abs(dp - ox); adn = Abs(dn - ox); d2x = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // - Y dc = pi + (sy*d1y); dp = pi + (sy*d1y) + pi2; dn = pi + (sy*d1y) - pi2; adc = Abs(dc - oy); adp = Abs(dp - oy); adn = Abs(dn - oy); d2y = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // - Z dc = pi + (sz*d1z); dp = pi + (sz*d1z) + pi2; dn = pi + (sz*d1z) - pi2; adc = Abs(dc - oz); adp = Abs(dp - oz); adn = Abs(dn - oz); d2z = (adc < adp) ? ((adc < adn) ? dc : dn) : ((adp < adn) ? dp : dn); // Find closest to reference triplet (prot) Real d1a, d2a; d1a = Abs(d1x - ox); d2a = Abs(d2x - ox); x = (d1a < d2a) ? d1x : d2x; d1a = Abs(d1y - oy); d2a = Abs(d2y - oy); y = (d1a < d2a) ? d1y : d2y; d1a = Abs(d1z - oz); d2a = Abs(d2z - oz); z = -((d1a < d2a) ? d1z : d2z); GePrint("Rotation ("+String(RotOrdStr[rotorder])+") = "+VectorToString(Vector(ToDeg(x),ToDeg(y),ToDeg(z)))+" (xyz)"); }