Blender V2.61 - r43446
|
00001 #include "MT_Optimize.h" 00002 00003 /* 00004 * This is a supposedly faster inverter than the cofactor 00005 * computation. It uses an LU decomposition sort of thing. */ 00006 GEN_INLINE void MT_Matrix4x4::invert() { 00007 /* normalize row 0 */ 00008 00009 int i,j,k; 00010 00011 for (i=1; i < 4; i++) m_el[0][i] /= m_el[0][0]; 00012 for (i=1; i < 4; i++) { 00013 for (j=i; j < 4; j++) { // do a column of L 00014 MT_Scalar sum = 0.0; 00015 for (k = 0; k < i; k++) 00016 sum += m_el[j][k] * m_el[k][i]; 00017 m_el[j][i] -= sum; 00018 } 00019 if (i == 3) continue; 00020 for (j=i+1; j < 4; j++) { // do a row of U 00021 MT_Scalar sum = 0.0; 00022 for (k = 0; k < i; k++) 00023 sum += m_el[i][k]*m_el[k][j]; 00024 m_el[i][j] = 00025 (m_el[i][j]-sum) / m_el[i][i]; 00026 } 00027 } 00028 for (i = 0; i < 4; i++ ) // invert L 00029 for (j = i; j < 4; j++ ) { 00030 MT_Scalar x = 1.0; 00031 if ( i != j ) { 00032 x = 0.0; 00033 for (k = i; k < j; k++ ) 00034 x -= m_el[j][k]*m_el[k][i]; 00035 } 00036 m_el[j][i] = x / m_el[j][j]; 00037 } 00038 for (i = 0; i < 4; i++ ) // invert U 00039 for (j = i; j < 4; j++ ) { 00040 if ( i == j ) continue; 00041 MT_Scalar sum = 0.0; 00042 for (k = i; k < j; k++ ) 00043 sum += m_el[k][j]*( (i==k) ? 1.0 : m_el[i][k] ); 00044 m_el[i][j] = -sum; 00045 } 00046 for (i = 0; i < 4; i++ ) // final inversion 00047 for (j = 0; j < 4; j++ ) { 00048 MT_Scalar sum = 0.0; 00049 for (k = ((i>j)?i:j); k < 4; k++ ) 00050 sum += ((j==k)?1.0:m_el[j][k])*m_el[k][i]; 00051 m_el[j][i] = sum; 00052 } 00053 } 00054 00055 GEN_INLINE MT_Matrix4x4 MT_Matrix4x4::inverse() const 00056 { 00057 MT_Matrix4x4 invmat = *this; 00058 00059 invmat.invert(); 00060 00061 return invmat; 00062 } 00063 00064 GEN_INLINE MT_Matrix4x4& MT_Matrix4x4::operator*=(const MT_Matrix4x4& m) 00065 { 00066 setValue(m.tdot(0, m_el[0]), m.tdot(1, m_el[0]), m.tdot(2, m_el[0]), m.tdot(3, m_el[0]), 00067 m.tdot(0, m_el[1]), m.tdot(1, m_el[1]), m.tdot(2, m_el[1]), m.tdot(3, m_el[1]), 00068 m.tdot(0, m_el[2]), m.tdot(1, m_el[2]), m.tdot(2, m_el[2]), m.tdot(3, m_el[2]), 00069 m.tdot(0, m_el[3]), m.tdot(1, m_el[3]), m.tdot(2, m_el[3]), m.tdot(3, m_el[3])); 00070 return *this; 00071 00072 } 00073 00074 GEN_INLINE MT_Vector4 operator*(const MT_Matrix4x4& m, const MT_Vector4& v) { 00075 return MT_Vector4(MT_dot(m[0], v), MT_dot(m[1], v), MT_dot(m[2], v), MT_dot(m[3], v)); 00076 } 00077 00078 GEN_INLINE MT_Vector4 operator*(const MT_Vector4& v, const MT_Matrix4x4& m) { 00079 return MT_Vector4(m.tdot(0, v), m.tdot(1, v), m.tdot(2, v), m.tdot(3, v)); 00080 } 00081 00082 GEN_INLINE MT_Matrix4x4 operator*(const MT_Matrix4x4& m1, const MT_Matrix4x4& m2) { 00083 return 00084 MT_Matrix4x4(m2.tdot(0, m1[0]), m2.tdot(1, m1[0]), m2.tdot(2, m1[0]), m2.tdot(3, m1[0]), 00085 m2.tdot(0, m1[1]), m2.tdot(1, m1[1]), m2.tdot(2, m1[1]), m2.tdot(3, m1[1]), 00086 m2.tdot(0, m1[2]), m2.tdot(1, m1[2]), m2.tdot(2, m1[2]), m2.tdot(3, m1[2]), 00087 m2.tdot(0, m1[3]), m2.tdot(1, m1[3]), m2.tdot(2, m1[3]), m2.tdot(3, m1[3])); 00088 } 00089 00090 00091 GEN_INLINE MT_Matrix4x4 MT_Matrix4x4::transposed() const { 00092 return MT_Matrix4x4(m_el[0][0], m_el[1][0], m_el[2][0], m_el[3][0], 00093 m_el[0][1], m_el[1][1], m_el[2][1], m_el[3][1], 00094 m_el[0][2], m_el[1][2], m_el[2][2], m_el[3][2], 00095 m_el[0][3], m_el[1][3], m_el[2][3], m_el[3][3]); 00096 } 00097 00098 GEN_INLINE void MT_Matrix4x4::transpose() { 00099 *this = transposed(); 00100 } 00101 00102 GEN_INLINE MT_Matrix4x4 MT_Matrix4x4::absolute() const { 00103 return 00104 MT_Matrix4x4(MT_abs(m_el[0][0]), MT_abs(m_el[0][1]), MT_abs(m_el[0][2]), MT_abs(m_el[0][3]), 00105 MT_abs(m_el[1][0]), MT_abs(m_el[1][1]), MT_abs(m_el[1][2]), MT_abs(m_el[1][3]), 00106 MT_abs(m_el[2][0]), MT_abs(m_el[2][1]), MT_abs(m_el[2][2]), MT_abs(m_el[2][3]), 00107 MT_abs(m_el[3][0]), MT_abs(m_el[3][1]), MT_abs(m_el[3][2]), MT_abs(m_el[3][3])); 00108 }