diff options
Diffstat (limited to 'core/math/matrix3.cpp')
-rw-r--r-- | core/math/matrix3.cpp | 380 |
1 files changed, 263 insertions, 117 deletions
diff --git a/core/math/matrix3.cpp b/core/math/matrix3.cpp index 71e6b62212..1fabfbbd4c 100644 --- a/core/math/matrix3.cpp +++ b/core/math/matrix3.cpp @@ -5,7 +5,7 @@ /* GODOT ENGINE */ /* http://www.godotengine.org */ /*************************************************************************/ -/* Copyright (c) 2007-2016 Juan Linietsky, Ariel Manzur. */ +/* Copyright (c) 2007-2017 Juan Linietsky, Ariel Manzur. */ /* */ /* Permission is hereby granted, free of charge, to any person obtaining */ /* a copy of this software and associated documentation files (the */ @@ -33,7 +33,7 @@ #define cofac(row1,col1, row2, col2)\ (elements[row1][col1] * elements[row2][col2] - elements[row1][col2] * elements[row2][col1]) -void Matrix3::from_z(const Vector3& p_z) { +void Basis::from_z(const Vector3& p_z) { if (Math::abs(p_z.z) > Math_SQRT12 ) { @@ -53,7 +53,7 @@ void Matrix3::from_z(const Vector3& p_z) { elements[2]=p_z; } -void Matrix3::invert() { +void Basis::invert() { real_t co[3]={ @@ -72,7 +72,8 @@ void Matrix3::invert() { } -void Matrix3::orthonormalize() { +void Basis::orthonormalize() { + ERR_FAIL_COND(determinant() == 0); // Gram-Schmidt Process @@ -92,100 +93,230 @@ void Matrix3::orthonormalize() { } -Matrix3 Matrix3::orthonormalized() const { +Basis Basis::orthonormalized() const { - Matrix3 c = *this; + Basis c = *this; c.orthonormalize(); return c; } +bool Basis::is_orthogonal() const { + Basis id; + Basis m = (*this)*transposed(); -Matrix3 Matrix3::inverse() const { + return isequal_approx(id,m); +} + +bool Basis::is_rotation() const { + return Math::isequal_approx(determinant(), 1) && is_orthogonal(); +} + + +bool Basis::is_symmetric() const { + + if (Math::abs(elements[0][1] - elements[1][0]) > CMP_EPSILON) + return false; + if (Math::abs(elements[0][2] - elements[2][0]) > CMP_EPSILON) + return false; + if (Math::abs(elements[1][2] - elements[2][1]) > CMP_EPSILON) + return false; + + return true; +} + + +Basis Basis::diagonalize() { + + //NOTE: only implemented for symmetric matrices + //with the Jacobi iterative method method + + ERR_FAIL_COND_V(!is_symmetric(), Basis()); + + const int ite_max = 1024; + + real_t off_matrix_norm_2 = elements[0][1] * elements[0][1] + elements[0][2] * elements[0][2] + elements[1][2] * elements[1][2]; + + int ite = 0; + Basis acc_rot; + while (off_matrix_norm_2 > CMP_EPSILON2 && ite++ < ite_max ) { + real_t el01_2 = elements[0][1] * elements[0][1]; + real_t el02_2 = elements[0][2] * elements[0][2]; + real_t el12_2 = elements[1][2] * elements[1][2]; + // Find the pivot element + int i, j; + if (el01_2 > el02_2) { + if (el12_2 > el01_2) { + i = 1; + j = 2; + } else { + i = 0; + j = 1; + } + } else { + if (el12_2 > el02_2) { + i = 1; + j = 2; + } else { + i = 0; + j = 2; + } + } + + // Compute the rotation angle + real_t angle; + if (Math::abs(elements[j][j] - elements[i][i]) < CMP_EPSILON) { + angle = Math_PI / 4; + } else { + angle = 0.5 * Math::atan(2 * elements[i][j] / (elements[j][j] - elements[i][i])); + } - Matrix3 inv=*this; + // Compute the rotation matrix + Basis rot; + rot.elements[i][i] = rot.elements[j][j] = Math::cos(angle); + rot.elements[i][j] = - (rot.elements[j][i] = Math::sin(angle)); + + // Update the off matrix norm + off_matrix_norm_2 -= elements[i][j] * elements[i][j]; + + // Apply the rotation + *this = rot * *this * rot.transposed(); + acc_rot = rot * acc_rot; + } + + return acc_rot; +} + +Basis Basis::inverse() const { + + Basis inv=*this; inv.invert(); return inv; } -void Matrix3::transpose() { +void Basis::transpose() { SWAP(elements[0][1],elements[1][0]); SWAP(elements[0][2],elements[2][0]); SWAP(elements[1][2],elements[2][1]); } -Matrix3 Matrix3::transposed() const { +Basis Basis::transposed() const { - Matrix3 tr=*this; + Basis tr=*this; tr.transpose(); return tr; } -void Matrix3::scale(const Vector3& p_scale) { +// Multiplies the matrix from left by the scaling matrix: M -> S.M +// See the comment for Basis::rotated for further explanation. +void Basis::scale(const Vector3& p_scale) { elements[0][0]*=p_scale.x; - elements[1][0]*=p_scale.x; - elements[2][0]*=p_scale.x; - elements[0][1]*=p_scale.y; + elements[0][1]*=p_scale.x; + elements[0][2]*=p_scale.x; + elements[1][0]*=p_scale.y; elements[1][1]*=p_scale.y; - elements[2][1]*=p_scale.y; - elements[0][2]*=p_scale.z; - elements[1][2]*=p_scale.z; + elements[1][2]*=p_scale.y; + elements[2][0]*=p_scale.z; + elements[2][1]*=p_scale.z; elements[2][2]*=p_scale.z; } -Matrix3 Matrix3::scaled( const Vector3& p_scale ) const { +Basis Basis::scaled( const Vector3& p_scale ) const { - Matrix3 m = *this; + Basis m = *this; m.scale(p_scale); return m; } -Vector3 Matrix3::get_scale() const { - - return Vector3( +Vector3 Basis::get_scale() const { + // We are assuming M = R.S, and performing a polar decomposition to extract R and S. + // FIXME: We eventually need a proper polar decomposition. + // As a cheap workaround until then, to ensure that R is a proper rotation matrix with determinant +1 + // (such that it can be represented by a Quat or Euler angles), we absorb the sign flip into the scaling matrix. + // As such, it works in conjuction with get_rotation(). + real_t det_sign = determinant() > 0 ? 1 : -1; + return det_sign*Vector3( Vector3(elements[0][0],elements[1][0],elements[2][0]).length(), Vector3(elements[0][1],elements[1][1],elements[2][1]).length(), Vector3(elements[0][2],elements[1][2],elements[2][2]).length() ); } -void Matrix3::rotate(const Vector3& p_axis, real_t p_phi) { - *this = *this * Matrix3(p_axis, p_phi); +// Multiplies the matrix from left by the rotation matrix: M -> R.M +// Note that this does *not* rotate the matrix itself. +// +// The main use of Basis is as Transform.basis, which is used a the transformation matrix +// of 3D object. Rotate here refers to rotation of the object (which is R * (*this)), +// not the matrix itself (which is R * (*this) * R.transposed()). +Basis Basis::rotated(const Vector3& p_axis, real_t p_phi) const { + return Basis(p_axis, p_phi) * (*this); } -Matrix3 Matrix3::rotated(const Vector3& p_axis, real_t p_phi) const { +void Basis::rotate(const Vector3& p_axis, real_t p_phi) { + *this = rotated(p_axis, p_phi); +} - return *this * Matrix3(p_axis, p_phi); +Basis Basis::rotated(const Vector3& p_euler) const { + return Basis(p_euler) * (*this); +} +void Basis::rotate(const Vector3& p_euler) { + *this = rotated(p_euler); } -Vector3 Matrix3::get_euler() const { +Vector3 Basis::get_rotation() const { + // Assumes that the matrix can be decomposed into a proper rotation and scaling matrix as M = R.S, + // and returns the Euler angles corresponding to the rotation part, complementing get_scale(). + // See the comment in get_scale() for further information. + Basis m = orthonormalized(); + real_t det = m.determinant(); + if (det < 0) { + // Ensure that the determinant is 1, such that result is a proper rotation matrix which can be represented by Euler angles. + m.scale(Vector3(-1,-1,-1)); + } + + return m.get_euler(); +} +// get_euler returns a vector containing the Euler angles in the format +// (a1,a2,a3), where a3 is the angle of the first rotation, and a1 is the last +// (following the convention they are commonly defined in the literature). +// +// The current implementation uses XYZ convention (Z is the first rotation), +// so euler.z is the angle of the (first) rotation around Z axis and so on, +// +// And thus, assuming the matrix is a rotation matrix, this function returns +// the angles in the decomposition R = X(a1).Y(a2).Z(a3) where Z(a) rotates +// around the z-axis by a and so on. +Vector3 Basis::get_euler() const { + + // Euler angles in XYZ convention. + // See https://en.wikipedia.org/wiki/Euler_angles#Rotation_matrix + // // rot = cy*cz -cy*sz sy - // cz*sx*sy+cx*sz cx*cz-sx*sy*sz -cy*sx - // -cx*cz*sy+sx*sz cz*sx+cx*sy*sz cx*cy - - Matrix3 m = *this; - m.orthonormalize(); + // cz*sx*sy+cx*sz cx*cz-sx*sy*sz -cy*sx + // -cx*cz*sy+sx*sz cz*sx+cx*sy*sz cx*cy Vector3 euler; - euler.y = Math::asin(m[0][2]); + ERR_FAIL_COND_V(is_rotation() == false, euler); + + euler.y = Math::asin(elements[0][2]); if ( euler.y < Math_PI*0.5) { if ( euler.y > -Math_PI*0.5) { - euler.x = Math::atan2(-m[1][2],m[2][2]); - euler.z = Math::atan2(-m[0][1],m[0][0]); + euler.x = Math::atan2(-elements[1][2],elements[2][2]); + euler.z = Math::atan2(-elements[0][1],elements[0][0]); } else { - real_t r = Math::atan2(m[1][0],m[1][1]); + real_t r = Math::atan2(elements[1][0],elements[1][1]); euler.z = 0.0; euler.x = euler.z - r; } } else { - real_t r = Math::atan2(m[0][1],m[1][1]); + real_t r = Math::atan2(elements[0][1],elements[1][1]); euler.z = 0; euler.x = r - euler.z; } @@ -195,43 +326,59 @@ Vector3 Matrix3::get_euler() const { } -void Matrix3::set_euler(const Vector3& p_euler) { +// set_euler expects a vector containing the Euler angles in the format +// (c,b,a), where a is the angle of the first rotation, and c is the last. +// The current implementation uses XYZ convention (Z is the first rotation). +void Basis::set_euler(const Vector3& p_euler) { real_t c, s; c = Math::cos(p_euler.x); s = Math::sin(p_euler.x); - Matrix3 xmat(1.0,0.0,0.0,0.0,c,-s,0.0,s,c); + Basis xmat(1.0,0.0,0.0,0.0,c,-s,0.0,s,c); c = Math::cos(p_euler.y); s = Math::sin(p_euler.y); - Matrix3 ymat(c,0.0,s,0.0,1.0,0.0,-s,0.0,c); + Basis ymat(c,0.0,s,0.0,1.0,0.0,-s,0.0,c); c = Math::cos(p_euler.z); s = Math::sin(p_euler.z); - Matrix3 zmat(c,-s,0.0,s,c,0.0,0.0,0.0,1.0); + Basis zmat(c,-s,0.0,s,c,0.0,0.0,0.0,1.0); //optimizer will optimize away all this anyway *this = xmat*(ymat*zmat); } -bool Matrix3::operator==(const Matrix3& p_matrix) const { +bool Basis::isequal_approx(const Basis& a, const Basis& b) const { + + for (int i=0;i<3;i++) { + for (int j=0;j<3;j++) { + if (Math::isequal_approx(a.elements[i][j],b.elements[i][j]) == false) + return false; + } + } + + return true; +} + +bool Basis::operator==(const Basis& p_matrix) const { for (int i=0;i<3;i++) { for (int j=0;j<3;j++) { - if (elements[i][j]!=p_matrix.elements[i][j]) + if (elements[i][j] != p_matrix.elements[i][j]) return false; } } return true; } -bool Matrix3::operator!=(const Matrix3& p_matrix) const { + +bool Basis::operator!=(const Basis& p_matrix) const { return (!(*this==p_matrix)); } -Matrix3::operator String() const { +Basis::operator String() const { String mtx; for (int i=0;i<3;i++) { @@ -248,12 +395,10 @@ Matrix3::operator String() const { return mtx; } -Matrix3::operator Quat() const { +Basis::operator Quat() const { + ERR_FAIL_COND_V(is_rotation() == false, Quat()); - Matrix3 m=*this; - m.orthonormalize(); - - real_t trace = m.elements[0][0] + m.elements[1][1] + m.elements[2][2]; + real_t trace = elements[0][0] + elements[1][1] + elements[2][2]; real_t temp[4]; if (trace > 0.0) @@ -262,66 +407,66 @@ Matrix3::operator Quat() const { temp[3]=(s * 0.5); s = 0.5 / s; - temp[0]=((m.elements[2][1] - m.elements[1][2]) * s); - temp[1]=((m.elements[0][2] - m.elements[2][0]) * s); - temp[2]=((m.elements[1][0] - m.elements[0][1]) * s); + temp[0]=((elements[2][1] - elements[1][2]) * s); + temp[1]=((elements[0][2] - elements[2][0]) * s); + temp[2]=((elements[1][0] - elements[0][1]) * s); } else { - int i = m.elements[0][0] < m.elements[1][1] ? - (m.elements[1][1] < m.elements[2][2] ? 2 : 1) : - (m.elements[0][0] < m.elements[2][2] ? 2 : 0); + int i = elements[0][0] < elements[1][1] ? + (elements[1][1] < elements[2][2] ? 2 : 1) : + (elements[0][0] < elements[2][2] ? 2 : 0); int j = (i + 1) % 3; int k = (i + 2) % 3; - real_t s = Math::sqrt(m.elements[i][i] - m.elements[j][j] - m.elements[k][k] + 1.0); + real_t s = Math::sqrt(elements[i][i] - elements[j][j] - elements[k][k] + 1.0); temp[i] = s * 0.5; s = 0.5 / s; - temp[3] = (m.elements[k][j] - m.elements[j][k]) * s; - temp[j] = (m.elements[j][i] + m.elements[i][j]) * s; - temp[k] = (m.elements[k][i] + m.elements[i][k]) * s; + temp[3] = (elements[k][j] - elements[j][k]) * s; + temp[j] = (elements[j][i] + elements[i][j]) * s; + temp[k] = (elements[k][i] + elements[i][k]) * s; } return Quat(temp[0],temp[1],temp[2],temp[3]); } -static const Matrix3 _ortho_bases[24]={ - Matrix3(1, 0, 0, 0, 1, 0, 0, 0, 1), - Matrix3(0, -1, 0, 1, 0, 0, 0, 0, 1), - Matrix3(-1, 0, 0, 0, -1, 0, 0, 0, 1), - Matrix3(0, 1, 0, -1, 0, 0, 0, 0, 1), - Matrix3(1, 0, 0, 0, 0, -1, 0, 1, 0), - Matrix3(0, 0, 1, 1, 0, 0, 0, 1, 0), - Matrix3(-1, 0, 0, 0, 0, 1, 0, 1, 0), - Matrix3(0, 0, -1, -1, 0, 0, 0, 1, 0), - Matrix3(1, 0, 0, 0, -1, 0, 0, 0, -1), - Matrix3(0, 1, 0, 1, 0, 0, 0, 0, -1), - Matrix3(-1, 0, 0, 0, 1, 0, 0, 0, -1), - Matrix3(0, -1, 0, -1, 0, 0, 0, 0, -1), - Matrix3(1, 0, 0, 0, 0, 1, 0, -1, 0), - Matrix3(0, 0, -1, 1, 0, 0, 0, -1, 0), - Matrix3(-1, 0, 0, 0, 0, -1, 0, -1, 0), - Matrix3(0, 0, 1, -1, 0, 0, 0, -1, 0), - Matrix3(0, 0, 1, 0, 1, 0, -1, 0, 0), - Matrix3(0, -1, 0, 0, 0, 1, -1, 0, 0), - Matrix3(0, 0, -1, 0, -1, 0, -1, 0, 0), - Matrix3(0, 1, 0, 0, 0, -1, -1, 0, 0), - Matrix3(0, 0, 1, 0, -1, 0, 1, 0, 0), - Matrix3(0, 1, 0, 0, 0, 1, 1, 0, 0), - Matrix3(0, 0, -1, 0, 1, 0, 1, 0, 0), - Matrix3(0, -1, 0, 0, 0, -1, 1, 0, 0) +static const Basis _ortho_bases[24]={ + Basis(1, 0, 0, 0, 1, 0, 0, 0, 1), + Basis(0, -1, 0, 1, 0, 0, 0, 0, 1), + Basis(-1, 0, 0, 0, -1, 0, 0, 0, 1), + Basis(0, 1, 0, -1, 0, 0, 0, 0, 1), + Basis(1, 0, 0, 0, 0, -1, 0, 1, 0), + Basis(0, 0, 1, 1, 0, 0, 0, 1, 0), + Basis(-1, 0, 0, 0, 0, 1, 0, 1, 0), + Basis(0, 0, -1, -1, 0, 0, 0, 1, 0), + Basis(1, 0, 0, 0, -1, 0, 0, 0, -1), + Basis(0, 1, 0, 1, 0, 0, 0, 0, -1), + Basis(-1, 0, 0, 0, 1, 0, 0, 0, -1), + Basis(0, -1, 0, -1, 0, 0, 0, 0, -1), + Basis(1, 0, 0, 0, 0, 1, 0, -1, 0), + Basis(0, 0, -1, 1, 0, 0, 0, -1, 0), + Basis(-1, 0, 0, 0, 0, -1, 0, -1, 0), + Basis(0, 0, 1, -1, 0, 0, 0, -1, 0), + Basis(0, 0, 1, 0, 1, 0, -1, 0, 0), + Basis(0, -1, 0, 0, 0, 1, -1, 0, 0), + Basis(0, 0, -1, 0, -1, 0, -1, 0, 0), + Basis(0, 1, 0, 0, 0, -1, -1, 0, 0), + Basis(0, 0, 1, 0, -1, 0, 1, 0, 0), + Basis(0, 1, 0, 0, 0, 1, 1, 0, 0), + Basis(0, 0, -1, 0, 1, 0, 1, 0, 0), + Basis(0, -1, 0, 0, 0, -1, 1, 0, 0) }; -int Matrix3::get_orthogonal_index() const { +int Basis::get_orthogonal_index() const { //could be sped up if i come up with a way - Matrix3 orth=*this; + Basis orth=*this; for(int i=0;i<3;i++) { for(int j=0;j<3;j++) { - float v = orth[i][j]; + real_t v = orth[i][j]; if (v>0.5) v=1.0; else if (v<-0.5) @@ -344,7 +489,7 @@ int Matrix3::get_orthogonal_index() const { return 0; } -void Matrix3::set_orthogonal_index(int p_index){ +void Basis::set_orthogonal_index(int p_index){ //there only exist 24 orthogonal bases in r3 ERR_FAIL_INDEX(p_index,24); @@ -355,12 +500,13 @@ void Matrix3::set_orthogonal_index(int p_index){ } -void Matrix3::get_axis_and_angle(Vector3 &r_axis,real_t& r_angle) const { +void Basis::get_axis_and_angle(Vector3 &r_axis,real_t& r_angle) const { + ERR_FAIL_COND(is_rotation() == false); - double angle,x,y,z; // variables for result - double epsilon = 0.01; // margin to allow for rounding errors - double epsilon2 = 0.1; // margin to distinguish between 0 and 180 degrees + real_t angle,x,y,z; // variables for result + real_t epsilon = 0.01; // margin to allow for rounding errors + real_t epsilon2 = 0.1; // margin to distinguish between 0 and 180 degrees if ( (Math::abs(elements[1][0]-elements[0][1])< epsilon) && (Math::abs(elements[2][0]-elements[0][2])< epsilon) @@ -379,12 +525,12 @@ void Matrix3::get_axis_and_angle(Vector3 &r_axis,real_t& r_angle) const { } // otherwise this singularity is angle = 180 angle = Math_PI; - double xx = (elements[0][0]+1)/2; - double yy = (elements[1][1]+1)/2; - double zz = (elements[2][2]+1)/2; - double xy = (elements[1][0]+elements[0][1])/4; - double xz = (elements[2][0]+elements[0][2])/4; - double yz = (elements[2][1]+elements[1][2])/4; + real_t xx = (elements[0][0]+1)/2; + real_t yy = (elements[1][1]+1)/2; + real_t zz = (elements[2][2]+1)/2; + real_t xy = (elements[1][0]+elements[0][1])/4; + real_t xz = (elements[2][0]+elements[0][2])/4; + real_t yz = (elements[2][1]+elements[1][2])/4; if ((xx > yy) && (xx > zz)) { // elements[0][0] is the largest diagonal term if (xx< epsilon) { x = 0; @@ -421,28 +567,27 @@ void Matrix3::get_axis_and_angle(Vector3 &r_axis,real_t& r_angle) const { return; } // as we have reached here there are no singularities so we can handle normally - double s = Math::sqrt((elements[1][2] - elements[2][1])*(elements[1][2] - elements[2][1]) + real_t s = Math::sqrt((elements[1][2] - elements[2][1])*(elements[1][2] - elements[2][1]) +(elements[2][0] - elements[0][2])*(elements[2][0] - elements[0][2]) - +(elements[0][1] - elements[1][0])*(elements[0][1] - elements[1][0])); // used to normalise - if (Math::abs(s) < 0.001) s=1; - // prevent divide by zero, should not happen if matrix is orthogonal and should be - // caught by singularity test above, but I've left it in just in case + +(elements[0][1] - elements[1][0])*(elements[0][1] - elements[1][0])); // s=|axis||sin(angle)|, used to normalise + angle = Math::acos(( elements[0][0] + elements[1][1] + elements[2][2] - 1)/2); - x = (elements[1][2] - elements[2][1])/s; - y = (elements[2][0] - elements[0][2])/s; - z = (elements[0][1] - elements[1][0])/s; + if (angle < 0) s = -s; + x = (elements[2][1] - elements[1][2])/s; + y = (elements[0][2] - elements[2][0])/s; + z = (elements[1][0] - elements[0][1])/s; r_axis=Vector3(x,y,z); r_angle=angle; } -Matrix3::Matrix3(const Vector3& p_euler) { +Basis::Basis(const Vector3& p_euler) { set_euler( p_euler ); } -Matrix3::Matrix3(const Quat& p_quat) { +Basis::Basis(const Quat& p_quat) { real_t d = p_quat.length_squared(); real_t s = 2.0 / d; @@ -456,7 +601,8 @@ Matrix3::Matrix3(const Quat& p_quat) { } -Matrix3::Matrix3(const Vector3& p_axis, real_t p_phi) { +Basis::Basis(const Vector3& p_axis, real_t p_phi) { + // Rotation matrix from axis and angle, see https://en.wikipedia.org/wiki/Rotation_matrix#Rotation_matrix_from_axis_and_angle Vector3 axis_sq(p_axis.x*p_axis.x,p_axis.y*p_axis.y,p_axis.z*p_axis.z); @@ -464,15 +610,15 @@ Matrix3::Matrix3(const Vector3& p_axis, real_t p_phi) { real_t sine= Math::sin(p_phi); elements[0][0] = axis_sq.x + cosine * ( 1.0 - axis_sq.x ); - elements[0][1] = p_axis.x * p_axis.y * ( 1.0 - cosine ) + p_axis.z * sine; - elements[0][2] = p_axis.z * p_axis.x * ( 1.0 - cosine ) - p_axis.y * sine; + elements[0][1] = p_axis.x * p_axis.y * ( 1.0 - cosine ) - p_axis.z * sine; + elements[0][2] = p_axis.z * p_axis.x * ( 1.0 - cosine ) + p_axis.y * sine; - elements[1][0] = p_axis.x * p_axis.y * ( 1.0 - cosine ) - p_axis.z * sine; + elements[1][0] = p_axis.x * p_axis.y * ( 1.0 - cosine ) + p_axis.z * sine; elements[1][1] = axis_sq.y + cosine * ( 1.0 - axis_sq.y ); - elements[1][2] = p_axis.y * p_axis.z * ( 1.0 - cosine ) + p_axis.x * sine; + elements[1][2] = p_axis.y * p_axis.z * ( 1.0 - cosine ) - p_axis.x * sine; - elements[2][0] = p_axis.z * p_axis.x * ( 1.0 - cosine ) + p_axis.y * sine; - elements[2][1] = p_axis.y * p_axis.z * ( 1.0 - cosine ) - p_axis.x * sine; + elements[2][0] = p_axis.z * p_axis.x * ( 1.0 - cosine ) - p_axis.y * sine; + elements[2][1] = p_axis.y * p_axis.z * ( 1.0 - cosine ) + p_axis.x * sine; elements[2][2] = axis_sq.z + cosine * ( 1.0 - axis_sq.z ); } |