## [cgkit-commits] SF.net SVN: cgkit:[252] cgkit/trunk

 [cgkit-commits] SF.net SVN: cgkit:[252] cgkit/trunk From: - 2008-08-25 15:09:02 ```Revision: 252 http://cgkit.svn.sourceforge.net/cgkit/?rev=252&view=rev Author: mbaas Date: 2008-08-25 15:08:52 +0000 (Mon, 25 Aug 2008) Log Message: ----------- Replaced the broken toEuler*() methods with a more robust version (taken from the comp.graphics.algorithms FAQ wiki). Modified Paths: -------------- cgkit/trunk/cgkit/light/cgtypes/mat3.py cgkit/trunk/changelog.txt cgkit/trunk/supportlib/include/mat3.h cgkit/trunk/supportlib/include/vec3.h cgkit/trunk/unittests/test_mat3.py cgkit/trunk/unittests/test_mat3_light.py Modified: cgkit/trunk/cgkit/light/cgtypes/mat3.py =================================================================== --- cgkit/trunk/cgkit/light/cgtypes/mat3.py 2008-08-25 15:05:45 UTC (rev 251) +++ cgkit/trunk/cgkit/light/cgtypes/mat3.py 2008-08-25 15:08:52 UTC (rev 252) @@ -595,140 +595,32 @@ def toEulerYXZ(self): """Return the Euler angles of a rotation matrix.""" - global _epsilon - - r1 = self.getRow(0) - r2 = self.getRow(1) - r3 = self.getRow(2) - - B = -r2.z - - x = math.asin(B) - - A = math.cos(x) - - if (A>_epsilon): - y = math.acos(r3.z/A) - z = math.acos(r2.y/A) - else: - z = 0.0 - y = math.acos(r1.x) - + y,x,z = self._getRotation(2, True, True, True) return (x,y,z) def toEulerZXY(self): """Return the Euler angles of a rotation matrix.""" - global _epsilon - - r1 = self.getRow(0) - r2 = self.getRow(1) - r3 = self.getRow(2) - - B = r3.y - - x = math.asin(B) - - A = math.cos(x) - - if (A>_epsilon): - y = math.acos(r3.z/A) - z = math.acos(r2.y/A) - else: - z = 0.0 - y = math.acos(r1.x) - + z,x,y = self._getRotation(1, False, True, True) return (x,y,z) def toEulerZYX(self): """Return the Euler angles of a rotation matrix.""" - global _epsilon - - r1 = self.getRow(0) - r2 = self.getRow(1) - r3 = self.getRow(2) - - D = -r3.x - - y = math.asin(D) - - C = math.cos(y) - - if (C>_epsilon): - x = math.acos(r3.z/C) - z = math.acos(r1.x/C) - else: - z = 0.0 - x = math.acos(-r2.y) - + z,y,x = self._getRotation(0, True, True, True) return (x,y,z) def toEulerYZX(self): """Return the Euler angles of a rotation matrix.""" - global _epsilon - - r1 = self.getRow(0) - r2 = self.getRow(1) - r3 = self.getRow(2) - - F = r2.x - - z = math.asin(F) - - E = math.cos(z) - - if (E>_epsilon): - x = math.acos(r2.y/E) - y = math.acos(r1.x/E) - else: - y = 0.0 - x = math.asin(r3.y) - + y,z,x = self._getRotation(0, False, True, True) return (x,y,z) def toEulerXZY(self): """Return the Euler angles of a rotation matrix.""" - global _epsilon - - r1 = self.getRow(0) - r2 = self.getRow(1) - r3 = self.getRow(2) - - F = -r1.y - - z = math.asin(F) - - E = math.cos(z) - - if (E>_epsilon): - x = math.acos(r2.y/E) - y = math.acos(r1.x/E) - else: - y = 0.0 - x = math.acos(r3.z) - + x,z,y = self._getRotation(1, True, True, True) return (x,y,z) def toEulerXYZ(self): """Return the Euler angles of a rotation matrix.""" - global _epsilon - - r1 = self.getRow(0) - r2 = self.getRow(1) - r3 = self.getRow(2) - - D = r1.z - - y = math.asin(D) - - C = math.cos(y) - - if (C>_epsilon): - x = math.acos(r3.z/C) - z = math.acos(r1.x/C) - else: - z = 0.0 - x = math.acos(r2.y) - + x,y,z = self._getRotation(2, False, True, True) return (x,y,z) def fromToRotation(_from, to): @@ -809,6 +701,94 @@ fromToRotation = staticmethod(fromToRotation) + def _getRotation(self, i, neg, alt, rev): + """Get Euler angles in any of the 24 different conventions. + + The first four argument select a particular convention. The last three + output arguments receive the angles. The order of the angles depends + on the convention. + + See http://www.cgafaq.info/wiki/Euler_angles_from_matrix for the + algorithm used. + + i: The index of the first axis (global rotations, s) or last axis (local rotations, r). 0=XZX, 1=YXY, 2=ZYZ + neg: If true, the convention contains an odd permutation of the convention defined by i alone (i.e. the middle axis is replaced. For example, XZX -> XYX) + alt: If true, the first and last axes are different. Local rotations: The first axis changes. For example, XZX -> YZX + rev: If true, the first and last angle are exchanged. This toggles between global/local rotations. In all the concrete getRotation*() functions this is always true because all the functions assume local rotations. + """ + v = [self[0,i], self[1,i], self[2,i]] + + j,k,h = _eulerIndices(i, neg, alt) + + a = v[h] + b = v[k] + c,s,r = _eulerGivens(a, b) + v[h] = r + s1 = c*self[k,j] - s*self[h,j] + c1 = c*self[k,k] - s*self[h,k] + r1 = math.atan2(s1, c1) + r2 = math.atan2(v[j], v[i]) + r3 = math.atan2(s, c) + if alt: + r3 = -r3 + if neg: + r1 = -r1 + r2 = -r2 + r3 = -r3 + if rev: + tmp = r1 + r1 = r3 + r3 = tmp + return r1,r2,r3 + +def _eulerIndices(i, neg, alt): + """Helper function for _getRotation().""" + next = [1, 2, 0, 1] + j = next[i+int(neg)] + k = 3-i-j + h = next[k+(1^int(neg)^int(alt))] + return j,k,h + +def _eulerGivens(a, b): + """Helper function for _getRotation().""" + global _epsilon + + absa = abs(a) + absb = abs(b) + # b=0? + if absb<=_epsilon: + if a>=0: + c = 1.0 + else: + c = -1.0 + return (c, 0.0, absa) + # a=0? + elif absa<=_epsilon: + if b>=0: + s = 1.0 + else: + s = -1.0 + return (0.0, s, absb) + # General case + else: + if absb>absa: + t = a/b + u = math.sqrt(1.0+t*t) + if b<0: + u = -u + s = 1.0/u + c = s*t + r = b*u + else: + t = b/a + u = math.sqrt(1.0+t*t) + if (a<0): + u = -u + c = 1.0/u + s = c*t + r = a*u + return c,s,r + ###################################################################### if __name__=="__main__": Modified: cgkit/trunk/changelog.txt =================================================================== --- cgkit/trunk/changelog.txt 2008-08-25 15:05:45 UTC (rev 251) +++ cgkit/trunk/changelog.txt 2008-08-25 15:08:52 UTC (rev 252) @@ -5,6 +5,15 @@ Bug fixes/enhancements: +- mat3: The toEuler*() methods were broken. Depending on the angles, they + could return invalid values. +- sl/noise: When only cgkit light was installed, the sl module couldn't be + imported because it imports noise which has no pure Python equivalent yet. + There is a dummy noise module now in the light version which allows importing + the module. An error is only triggered when one of the functions gets called. +- rply: The sources can now also be compiled using a C++ compiler. Thanks + to Chris Foster for the patch! +- ribexport: Added RenderMan support for the torus geometry. - ri/cri: Passing unicode strings in a parameter list has caused an exception. ---------------------------------------------------------------------- Modified: cgkit/trunk/supportlib/include/mat3.h =================================================================== --- cgkit/trunk/supportlib/include/mat3.h 2008-08-25 15:05:45 UTC (rev 251) +++ cgkit/trunk/supportlib/include/mat3.h 2008-08-25 15:08:52 UTC (rev 252) @@ -119,6 +119,7 @@ mat3(const mat3& A); T& at(short i, short j); + const T& at(short i, short j) const; // set_ and get_ methods mat3& setIdentity(); @@ -219,6 +220,10 @@ vec3 r2; vec3 r3; + void _getRotation(int i, bool neg, bool alt, bool rev, T& r1, T& r2, T& r3) const; + void _eulerIndices(int i, bool neg, bool alt, int &j, int &k, int&h) const; + void _eulerGivens(T a, T b, T& c, T& s, T& v) const; + // Constructor with three vectors: Init the rows with those vectors. // This constructor is private because its semantics depends // on the actual implementation (that is, the three vectors are @@ -380,7 +385,30 @@ } } +/** + Return element [i,j]. + The indices range from 0 to 2. Other values will cause an EIndexError + exception to be thrown. + + @param i Row + @param j Column + @return Reference to element [i,j] + @exception EIndexError + */ +template +inline const T& mat3::at(short i, short j) const +{ + switch(i) + { + case 0: return r1[j]; + case 1: return r2[j]; + case 2: return r3[j]; + default: throw EIndexError(); + } +} + + /*----------------------------setIdentity------------------------------*//** Set identity matrix. @@ -1584,34 +1612,136 @@ } /** - Get Euler angles (order: yxz). + Helper function for the _getRotation() function. + */ +template +void mat3::_eulerIndices(int i, bool neg, bool alt, int& j, int& k, int& h) const +{ + int next[4] = {1, 2, 0, 1}; + j = next[i+int(neg)]; + k = 3-i-j; + h = next[k+(1^int(neg)^int(alt))]; +} - \param[out] x Angle around x (radians) - \param[out] y Angle around y (radians) - \param[out] z Angle around z (radians) +/** + Helper function for the _getRotation() function. */ template -void mat3::getRotationYXZ(T& x, T& y, T& z) const +void mat3::_eulerGivens(T a, T b, T& c, T& s, T& r) const { - T B = -r2.z; + T absa = xabs(a); + T absb = xabs(b); + // b=0? + if (absb<=vec3::epsilon) + { + if (a>=0) + c = 1.0; + else + c = -1.0; + s = 0.0; + r = absa; + } + // a=0? + else if (absa<=vec3::epsilon) + { + c = 0.0; + if (b>=0) + s = 1.0; + else + s = -1.0; + r = absb; + } + // General case + else + { + if (absb>absa) + { + T t = a/b; + T u = sqrt(1.0+t*t); + if (b<0) + u = -u; + s = 1.0/u; + c = s*t; + r = b*u; + } + else + { + T t = b/a; + T u = sqrt(1.0+t*t); + if (a<0) + u = -u; + c = 1.0/u; + s = c*t; + r = a*u; + } + } +} - x = asin(B); +/** + Get Euler angles in any of the 24 different conventions. - T A = cos(x); + The first four argument select a particular convention. The last three + output arguments receive the angles. The order of the angles depends + on the convention. - if (A>vec3::epsilon) + See http://www.cgafaq.info/wiki/Euler_angles_from_matrix for the + algorithm used. + + \param i The index of the first axis (global rotations, s) or last axis (local rotations, r). 0=XZX, 1=YXY, 2=ZYZ + \param neg If true, the convention contains an odd permutation of the convention defined by i alone (i.e. the middle axis is replaced. For example, XZX -> XYX) + \param alt If true, the first and last axes are different. Local rotations: The first axis changes. For example, XZX -> YZX + \param rev If true, the first and last angle are exchanged. This toggles between global/local rotations. In all the concrete getRotation*() functions this is always true because all the functions assume local rotations. + \param[out] r1 Angle around first axis (radians) + \param[out] r2 Angle around second axis (radians) + \param[out] r3 Angle around third axis (radians) + */ +template +void mat3::_getRotation(int i, bool neg, bool alt, bool rev, T& r1, T& r2, T& r3) const +{ + int j, k, h; + T v[3] = {at(0,i), at(1,i), at(2,i)}; + + _eulerIndices(i, neg, alt, j, k, h); + + T a = v[h]; + T b = v[k]; + T c, s; + _eulerGivens(a, b, c, s, v[h]); + T s1 = c*at(k,j) - s*at(h,j); + T c1 = c*at(k,k) - s*at(h,k); + r1 = atan2(s1, c1); + r2 = atan2(v[j], v[i]); + r3 = atan2(s, c); + if (alt) + r3 = -r3; + if (neg) { - y = acos(r3.z/A); - z = acos(r2.y/A); + r1 = -r1; + r2 = -r2; + r3 = -r3; } - else + if (rev) { - z = 0; - y = acos(r1.x); - } + T tmp = r1; + r1 = r3; + r3 = tmp; + } } /** + Get Euler angles (order: yxz). + + \param[out] x Angle around x (radians) + \param[out] y Angle around y (radians) + \param[out] z Angle around z (radians) + */ +template +void mat3::getRotationYXZ(T& x, T& y, T& z) const +{ + _getRotation(2, 1, 1, 1, y, x, z); +} + +/** Get Euler angles (order: zxy). \param[out] x Angle around x (radians) @@ -1621,22 +1751,7 @@ template void mat3::getRotationZXY(T& x, T& y, T& z) const { - T B = r3.y; - - x = asin(B); - - T A = cos(x); - - if (A>vec3::epsilon) - { - y = acos(r3.z/A); - z = acos(r2.y/A); - } - else - { - z = 0; - y = acos(r1.x); - } + _getRotation(1, 0, 1, 1, z, x, y); } /** @@ -1649,22 +1764,7 @@ template void mat3::getRotationZYX(T& x, T& y, T& z) const { - T D = -r3.x; - - y = asin(D); - - T C = cos(y); - - if (C>vec3::epsilon) - { - x = acos(r3.z/C); - z = acos(r1.x/C); - } - else - { - z = 0; - x = acos(-r2.y); - } + _getRotation(0, 1, 1, 1, z, y, x); } /** @@ -1677,22 +1777,7 @@ template void mat3::getRotationYZX(T& x, T& y, T& z) const { - T F = r2.x; - - z = asin(F); - - T E = cos(z); - - if (E>vec3::epsilon) - { - x = acos(r2.y/E); - y = acos(r1.x/E); - } - else - { - y = 0; - x = asin(r3.y); - } + _getRotation(0, 0, 1, 1, y, z, x); } /** @@ -1705,22 +1790,7 @@ template void mat3::getRotationXZY(T& x, T& y, T& z) const { - T F = -r1.y; - - z = asin(F); - - T E = cos(z); - - if (E>vec3::epsilon) - { - x = acos(r2.y/E); - y = acos(r1.x/E); - } - else - { - y = 0; - x = acos(r3.z); - } + _getRotation(1, 1, 1, 1, x, z, y); } /** @@ -1733,22 +1803,7 @@ template void mat3::getRotationXYZ(T& x, T& y, T& z) const { - T D = r1.z; - - y = asin(D); - - T C = cos(y); - - if (C>vec3::epsilon) - { - x = acos(r3.z/C); - z = acos(r1.x/C); - } - else - { - z = 0; - x = acos(r2.y); - } + _getRotation(2, 0, 1, 1, x, y, z); } /** Modified: cgkit/trunk/supportlib/include/vec3.h =================================================================== --- cgkit/trunk/supportlib/include/vec3.h 2008-08-25 15:05:45 UTC (rev 251) +++ cgkit/trunk/supportlib/include/vec3.h 2008-08-25 15:08:52 UTC (rev 252) @@ -170,6 +170,7 @@ void get_polar(T& r, T& theta, T& phi) const; T& operator[](int); + const T& operator[](int) const; vec3& operator+=(const vec3& v); vec3& operator-=(const vec3& v); @@ -415,6 +416,27 @@ } } +/*-----------------------------operator[]-------------------------------*//** + Index operator. + + Return component i. The index is checked. + + @param index Component index (0,1 or 2). + @return Component i. + @exception EIndexError Index is out of range. +*//*------------------------------------------------------------------------*/ +template +const T& vec3::operator[](int index) const +{ + switch (index) + { + case 0: return x; + case 1: return y; + case 2: return z; + default: throw EIndexError(); + } +} + /*-----------------------------operator+=-------------------------------*//** vector += vector. Modified: cgkit/trunk/unittests/test_mat3.py =================================================================== --- cgkit/trunk/unittests/test_mat3.py 2008-08-25 15:05:45 UTC (rev 251) +++ cgkit/trunk/unittests/test_mat3.py 2008-08-25 15:08:52 UTC (rev 252) @@ -1,10 +1,10 @@ # Test the mat3 class import unittest +#from cgkit.light.cgtypes import * from cgkit.cgtypes import * -#from cgkit.light.cgtypes import * import math, os, pickle, cPickle, sys -from cgkit.sl import radians +from cgkit.sl import degrees, radians class TestMat3(unittest.TestCase): @@ -511,26 +511,34 @@ matrix that is composed by 3 individual rotations. """ - angle = {"X":radians(20), "Y":radians(30), "Z":radians(40)} + angles = [{"X":radians(20), "Y":radians(30), "Z":radians(40)}, + {"X":radians(0), "Y":radians(0), "Z":radians(0)}, + {"X":radians(350), "Y":radians(0), "Z":radians(0)}, + {"X":radians(0), "Y":radians(350), "Z":radians(0)}, + {"X":radians(0), "Y":radians(0), "Z":radians(350)},] axis = {"X":vec3(1,0,0), "Y":vec3(0,1,0), "Z":vec3(0,0,1)} for order in ["XYZ", "YZX", "ZXY", "XZY", "YXZ", "ZYX"]: - R1 = mat3.rotation(angle[order[0]], axis[order[0]]) - R2 = mat3.rotation(angle[order[1]], axis[order[1]]) - R3 = mat3.rotation(angle[order[2]], axis[order[2]]) - # Each rotation is about the *global* axis, so these rotations - # have to be applied just in the opposite order than mentioned - # in the fromEuler*() method name. - C = R1*R2*R3 - exec 'E = mat3.fromEuler%s(angle["X"], angle["Y"], angle["Z"])'%order - self.assertEqual(E, C) + for angle in angles: + R1 = mat3.rotation(angle[order[0]], axis[order[0]]) + R2 = mat3.rotation(angle[order[1]], axis[order[1]]) + R3 = mat3.rotation(angle[order[2]], axis[order[2]]) + # Each rotation is about the *global* axis, so these rotations + # have to be applied just in the opposite order than mentioned + # in the fromEuler*() method name. + C = R1*R2*R3 + exec 'E = mat3.fromEuler%s(angle["X"], angle["Y"], angle["Z"])'%order + self.assertEqual(E, C) + + exec 'x,y,z = E.toEuler%s()'%order + exec 'E2 = mat3.fromEuler%s(x, y, z)'%order + if E2!=E: +# print E +# print E2 + msg = "The matrix E2 generated from the toEuler() angles doesn't match the original matrix E.\n" + msg += "Original angles: (%s, %s, %s), toEuler angles: (%s, %s, %s)"%(degrees(angle["X"]),degrees(angle["Y"]),degrees(angle["Z"]),degrees(x),degrees(y),degrees(z)) + self.fail(msg) - exec 'x,y,z = E.toEuler%s()'%order - self.assertAlmostEqual(x, angle["X"], 8) - self.assertAlmostEqual(y, angle["Y"], 8) - self.assertAlmostEqual(z, angle["Z"], 8) - - ###################################################################### def testFromToRotation(self): R = mat3.fromToRotation((1,0,0), (0,1,0)) Modified: cgkit/trunk/unittests/test_mat3_light.py =================================================================== --- cgkit/trunk/unittests/test_mat3_light.py 2008-08-25 15:05:45 UTC (rev 251) +++ cgkit/trunk/unittests/test_mat3_light.py 2008-08-25 15:08:52 UTC (rev 252) @@ -1,10 +1,9 @@ # Test the mat3 class import unittest -#from cgkit.cgtypes import * from cgkit.light.cgtypes import * import math, os, pickle, cPickle, sys -from cgkit.sl import radians +from cgkit.sl import degrees, radians class TestMat3(unittest.TestCase): @@ -511,26 +510,34 @@ matrix that is composed by 3 individual rotations. """ - angle = {"X":radians(20), "Y":radians(30), "Z":radians(40)} + angles = [{"X":radians(20), "Y":radians(30), "Z":radians(40)}, + {"X":radians(0), "Y":radians(0), "Z":radians(0)}, + {"X":radians(350), "Y":radians(0), "Z":radians(0)}, + {"X":radians(0), "Y":radians(350), "Z":radians(0)}, + {"X":radians(0), "Y":radians(0), "Z":radians(350)},] axis = {"X":vec3(1,0,0), "Y":vec3(0,1,0), "Z":vec3(0,0,1)} for order in ["XYZ", "YZX", "ZXY", "XZY", "YXZ", "ZYX"]: - R1 = mat3.rotation(angle[order[0]], axis[order[0]]) - R2 = mat3.rotation(angle[order[1]], axis[order[1]]) - R3 = mat3.rotation(angle[order[2]], axis[order[2]]) - # Each rotation is about the *global* axis, so these rotations - # have to be applied just in the opposite order than mentioned - # in the fromEuler*() method name. - C = R1*R2*R3 - exec 'E = mat3.fromEuler%s(angle["X"], angle["Y"], angle["Z"])'%order - self.assertEqual(E, C) + for angle in angles: + R1 = mat3.rotation(angle[order[0]], axis[order[0]]) + R2 = mat3.rotation(angle[order[1]], axis[order[1]]) + R3 = mat3.rotation(angle[order[2]], axis[order[2]]) + # Each rotation is about the *global* axis, so these rotations + # have to be applied just in the opposite order than mentioned + # in the fromEuler*() method name. + C = R1*R2*R3 + exec 'E = mat3.fromEuler%s(angle["X"], angle["Y"], angle["Z"])'%order + self.assertEqual(E, C) + + exec 'x,y,z = E.toEuler%s()'%order + exec 'E2 = mat3.fromEuler%s(x, y, z)'%order + if E2!=E: +# print E +# print E2 + msg = "The matrix E2 generated from the toEuler() angles doesn't match the original matrix E.\n" + msg += "Original angles: (%s, %s, %s), toEuler angles: (%s, %s, %s)"%(degrees(angle["X"]),degrees(angle["Y"]),degrees(angle["Z"]),degrees(x),degrees(y),degrees(z)) + self.fail(msg) - exec 'x,y,z = E.toEuler%s()'%order - self.assertAlmostEqual(x, angle["X"], 8) - self.assertAlmostEqual(y, angle["Y"], 8) - self.assertAlmostEqual(z, angle["Z"], 8) - - ###################################################################### def testFromToRotation(self): R = mat3.fromToRotation((1,0,0), (0,1,0)) This was sent by the SourceForge.net collaborative development platform, the world's largest Open Source development site. ```