31 #ifndef OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED 32 #define OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED 48 template<
typename T>
class Vec3;
49 template<
typename T>
class Mat4;
50 template<
typename T>
class Quat;
77 template<
typename Source>
78 Mat3(Source a, Source b, Source c,
79 Source d, Source e, Source f,
80 Source g, Source h, Source i)
82 MyBase::mm[0] =
static_cast<T
>(a);
83 MyBase::mm[1] =
static_cast<T
>(b);
84 MyBase::mm[2] =
static_cast<T
>(c);
85 MyBase::mm[3] =
static_cast<T
>(d);
86 MyBase::mm[4] =
static_cast<T
>(e);
87 MyBase::mm[5] =
static_cast<T
>(f);
88 MyBase::mm[6] =
static_cast<T
>(g);
89 MyBase::mm[7] =
static_cast<T
>(h);
90 MyBase::mm[8] =
static_cast<T
>(i);
95 template<
typename Source>
99 this->setRows(v1, v2, v3);
101 this->setColumns(v1, v2, v3);
109 template<
typename Source>
112 MyBase::mm[0] =
static_cast<T
>(a[0]);
113 MyBase::mm[1] =
static_cast<T
>(a[1]);
114 MyBase::mm[2] =
static_cast<T
>(a[2]);
115 MyBase::mm[3] =
static_cast<T
>(a[3]);
116 MyBase::mm[4] =
static_cast<T
>(a[4]);
117 MyBase::mm[5] =
static_cast<T
>(a[5]);
118 MyBase::mm[6] =
static_cast<T
>(a[6]);
119 MyBase::mm[7] =
static_cast<T
>(a[7]);
120 MyBase::mm[8] =
static_cast<T
>(a[8]);
126 for (
int i=0; i<3; ++i) {
127 for (
int j=0; j<3; ++j) {
128 MyBase::mm[i*3 + j] = m[i][j];
134 template<
typename Source>
137 for (
int i=0; i<3; ++i) {
138 for (
int j=0; j<3; ++j) {
139 MyBase::mm[i*3 + j] =
static_cast<T
>(m[i][j]);
147 for (
int i=0; i<3; ++i) {
148 for (
int j=0; j<3; ++j) {
149 MyBase::mm[i*3 + j] = m[i][j];
180 MyBase::mm[i3+0] = v[0];
181 MyBase::mm[i3+1] = v[1];
182 MyBase::mm[i3+2] = v[2];
189 return Vec3<T>((*this)(i,0), (*
this)(i,1), (*
this)(i,2));
196 MyBase::mm[0+j] = v[0];
197 MyBase::mm[3+j] = v[1];
198 MyBase::mm[6+j] = v[2];
205 return Vec3<T>((*this)(0,j), (*
this)(1,j), (*
this)(2,j));
213 T* operator[](
int i) {
return &(MyBase::mm[i*3]); }
216 const T*
operator[](
int i)
const {
return &(MyBase::mm[i*3]); }
229 return MyBase::mm[3*i+j];
239 return MyBase::mm[3*i+j];
245 MyBase::mm[0] = v1[0];
246 MyBase::mm[1] = v1[1];
247 MyBase::mm[2] = v1[2];
248 MyBase::mm[3] = v2[0];
249 MyBase::mm[4] = v2[1];
250 MyBase::mm[5] = v2[2];
251 MyBase::mm[6] = v3[0];
252 MyBase::mm[7] = v3[1];
253 MyBase::mm[8] = v3[2];
259 MyBase::mm[0] = v1[0];
260 MyBase::mm[1] = v2[0];
261 MyBase::mm[2] = v3[0];
262 MyBase::mm[3] = v1[1];
263 MyBase::mm[4] = v2[1];
264 MyBase::mm[5] = v3[1];
265 MyBase::mm[6] = v1[2];
266 MyBase::mm[7] = v2[2];
267 MyBase::mm[8] = v3[2];
273 MyBase::mm[0] = vdiag[0];
274 MyBase::mm[1] = vtri[0];
275 MyBase::mm[2] = vtri[1];
276 MyBase::mm[3] = vtri[0];
277 MyBase::mm[4] = vdiag[1];
278 MyBase::mm[5] = vtri[2];
279 MyBase::mm[6] = vtri[1];
280 MyBase::mm[7] = vtri[2];
281 MyBase::mm[8] = vdiag[2];
288 vdiag[0], vtri[0], vtri[1],
289 vtri[0], vdiag[1], vtri[2],
290 vtri[1], vtri[2], vdiag[2]
302 {*
this = rotation<Mat3<T> >(q);}
307 {*
this = rotation<Mat3<T> >(axis,
angle);}
338 template<
typename Source>
344 std::copy(src, (src + this->numElements()), MyBase::mm);
349 bool eq(
const Mat3 &m, T eps=1.0e-8)
const 366 -MyBase::mm[0], -MyBase::mm[1], -MyBase::mm[2],
367 -MyBase::mm[3], -MyBase::mm[4], -MyBase::mm[5],
368 -MyBase::mm[6], -MyBase::mm[7], -MyBase::mm[8]
378 template <
typename S>
381 MyBase::mm[0] *= scalar;
382 MyBase::mm[1] *= scalar;
383 MyBase::mm[2] *= scalar;
384 MyBase::mm[3] *= scalar;
385 MyBase::mm[4] *= scalar;
386 MyBase::mm[5] *= scalar;
387 MyBase::mm[6] *= scalar;
388 MyBase::mm[7] *= scalar;
389 MyBase::mm[8] *= scalar;
394 template <
typename S>
399 MyBase::mm[0] += s[0];
400 MyBase::mm[1] += s[1];
401 MyBase::mm[2] += s[2];
402 MyBase::mm[3] += s[3];
403 MyBase::mm[4] += s[4];
404 MyBase::mm[5] += s[5];
405 MyBase::mm[6] += s[6];
406 MyBase::mm[7] += s[7];
407 MyBase::mm[8] += s[8];
412 template <
typename S>
417 MyBase::mm[0] -= s[0];
418 MyBase::mm[1] -= s[1];
419 MyBase::mm[2] -= s[2];
420 MyBase::mm[3] -= s[3];
421 MyBase::mm[4] -= s[4];
422 MyBase::mm[5] -= s[5];
423 MyBase::mm[6] -= s[6];
424 MyBase::mm[7] -= s[7];
425 MyBase::mm[8] -= s[8];
430 template <
typename S>
437 MyBase::mm[0] =
static_cast<T
>(s0[0] * s1[0] +
440 MyBase::mm[1] =
static_cast<T
>(s0[0] * s1[1] +
443 MyBase::mm[2] =
static_cast<T
>(s0[0] * s1[2] +
447 MyBase::mm[3] =
static_cast<T
>(s0[3] * s1[0] +
450 MyBase::mm[4] =
static_cast<T
>(s0[3] * s1[1] +
453 MyBase::mm[5] =
static_cast<T
>(s0[3] * s1[2] +
457 MyBase::mm[6] =
static_cast<T
>(s0[6] * s1[0] +
460 MyBase::mm[7] =
static_cast<T
>(s0[6] * s1[1] +
463 MyBase::mm[8] =
static_cast<T
>(s0[6] * s1[2] +
474 MyBase::mm[4] * MyBase::mm[8] - MyBase::mm[5] * MyBase::mm[7],
475 MyBase::mm[5] * MyBase::mm[6] - MyBase::mm[3] * MyBase::mm[8],
476 MyBase::mm[3] * MyBase::mm[7] - MyBase::mm[4] * MyBase::mm[6],
477 MyBase::mm[2] * MyBase::mm[7] - MyBase::mm[1] * MyBase::mm[8],
478 MyBase::mm[0] * MyBase::mm[8] - MyBase::mm[2] * MyBase::mm[6],
479 MyBase::mm[1] * MyBase::mm[6] - MyBase::mm[0] * MyBase::mm[7],
480 MyBase::mm[1] * MyBase::mm[5] - MyBase::mm[2] * MyBase::mm[4],
481 MyBase::mm[2] * MyBase::mm[3] - MyBase::mm[0] * MyBase::mm[5],
482 MyBase::mm[0] * MyBase::mm[4] - MyBase::mm[1] * MyBase::mm[3]);
489 MyBase::mm[4] * MyBase::mm[8] - MyBase::mm[5] * MyBase::mm[7],
490 MyBase::mm[2] * MyBase::mm[7] - MyBase::mm[1] * MyBase::mm[8],
491 MyBase::mm[1] * MyBase::mm[5] - MyBase::mm[2] * MyBase::mm[4],
492 MyBase::mm[5] * MyBase::mm[6] - MyBase::mm[3] * MyBase::mm[8],
493 MyBase::mm[0] * MyBase::mm[8] - MyBase::mm[2] * MyBase::mm[6],
494 MyBase::mm[2] * MyBase::mm[3] - MyBase::mm[0] * MyBase::mm[5],
495 MyBase::mm[3] * MyBase::mm[7] - MyBase::mm[4] * MyBase::mm[6],
496 MyBase::mm[1] * MyBase::mm[6] - MyBase::mm[0] * MyBase::mm[7],
497 MyBase::mm[0] * MyBase::mm[4] - MyBase::mm[1] * MyBase::mm[3]);
505 MyBase::mm[0], MyBase::mm[3], MyBase::mm[6],
506 MyBase::mm[1], MyBase::mm[4], MyBase::mm[7],
507 MyBase::mm[2], MyBase::mm[5], MyBase::mm[8]);
517 const T det = inv.
mm[0]*MyBase::mm[0] + inv.
mm[1]*MyBase::mm[3] + inv.
mm[2]*MyBase::mm[6];
523 return inv * (T(1)/det);
529 const T co00 = MyBase::mm[4]*MyBase::mm[8] - MyBase::mm[5]*MyBase::mm[7];
530 const T co10 = MyBase::mm[5]*MyBase::mm[6] - MyBase::mm[3]*MyBase::mm[8];
531 const T co20 = MyBase::mm[3]*MyBase::mm[7] - MyBase::mm[4]*MyBase::mm[6];
532 return MyBase::mm[0]*co00 + MyBase::mm[1]*co10 + MyBase::mm[2]*co20;
538 return MyBase::mm[0]+MyBase::mm[4]+MyBase::mm[8];
552 template<
typename T0>
555 return static_cast< Vec3<T0> >(v * *
this);
560 template<
typename T0>
563 return static_cast< Vec3<T0> >(*
this * v);
573 ret.
mm[0] *= diag(0);
574 ret.
mm[1] *= diag(1);
575 ret.
mm[2] *= diag(2);
576 ret.
mm[3] *= diag(0);
577 ret.
mm[4] *= diag(1);
578 ret.
mm[5] *= diag(2);
579 ret.
mm[6] *= diag(0);
580 ret.
mm[7] *= diag(1);
581 ret.
mm[8] *= diag(2);
589 template <
typename T0,
typename T1>
595 for (
int i=0; i<9; ++i) {
603 template <
typename T0,
typename T1>
608 template <
typename S,
typename T>
614 template <
typename S,
typename T>
624 template <
typename T0,
typename T1>
634 template <
typename T0,
typename T1>
644 template <
typename T0,
typename T1>
654 template<
typename T,
typename MT>
660 _v[0]*m[0] + _v[1]*m[1] + _v[2]*m[2],
661 _v[0]*m[3] + _v[1]*m[4] + _v[2]*m[5],
662 _v[0]*m[6] + _v[1]*m[7] + _v[2]*m[8]);
667 template<
typename T,
typename MT>
673 _v[0]*m[0] + _v[1]*m[3] + _v[2]*m[6],
674 _v[0]*m[1] + _v[1]*m[4] + _v[2]*m[7],
675 _v[0]*m[2] + _v[1]*m[5] + _v[2]*m[8]);
680 template<
typename T,
typename MT>
690 template <
typename T>
693 return Mat3<T>(v1[0]*v2[0], v1[0]*v2[1], v1[0]*v2[2],
694 v1[1]*v2[0], v1[1]*v2[1], v1[1]*v2[2],
695 v1[2]*v2[0], v1[2]*v2[1], v1[2]*v2[2]);
702 template<
typename T,
typename T0>
712 namespace mat3_internal {
721 double cotan_of_2_theta;
723 double cosin_of_theta;
729 double Sjj_minus_Sii = D[j] - D[i];
732 tan_of_theta = Sij / Sjj_minus_Sii;
735 cotan_of_2_theta = 0.5*Sjj_minus_Sii / Sij ;
737 if (cotan_of_2_theta < 0.) {
739 -1./(sqrt(1. + cotan_of_2_theta*cotan_of_2_theta) - cotan_of_2_theta);
742 1./(sqrt(1. + cotan_of_2_theta*cotan_of_2_theta) + cotan_of_2_theta);
746 cosin_of_theta = 1./sqrt( 1. + tan_of_theta * tan_of_theta);
747 sin_of_theta = cosin_of_theta * tan_of_theta;
748 z = tan_of_theta * Sij;
752 for (
int k = 0; k < i; ++k) {
754 S(k,i) = cosin_of_theta * temp - sin_of_theta * S(k,j);
755 S(k,j)= sin_of_theta * temp + cosin_of_theta * S(k,j);
757 for (
int k = i+1; k < j; ++k) {
759 S(i,k) = cosin_of_theta * temp - sin_of_theta * S(k,j);
760 S(k,j) = sin_of_theta * temp + cosin_of_theta * S(k,j);
762 for (
int k = j+1; k < n; ++k) {
764 S(i,k) = cosin_of_theta * temp - sin_of_theta * S(j,k);
765 S(j,k) = sin_of_theta * temp + cosin_of_theta * S(j,k);
767 for (
int k = 0; k < n; ++k)
770 Q(k,i) = cosin_of_theta * temp - sin_of_theta*Q(k,j);
771 Q(k,j) = sin_of_theta * temp + cosin_of_theta*Q(k,j);
786 unsigned int MAX_ITERATIONS=250)
796 for (
int i = 0; i < n; ++i) {
800 unsigned int iterations(0);
807 for (
int i = 0; i < n; ++i) {
808 for (
int j = i+1; j < n; ++j) {
821 for (
int i = 0; i < n; ++i) {
822 for (
int j = i+1; j < n; ++j){
828 if (fabs(S(i,j)) > max_element) {
829 max_element = fabs(S(i,j));
836 }
while (iterations < MAX_ITERATIONS);
849 template<>
inline math::Mat3s zeroVal<math::Mat3s>() {
return math::Mat3s::zero(); }
850 template<>
inline math::Mat3d zeroVal<math::Mat3d>() {
return math::Mat3d::zero(); }
855 #endif // OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:415
void pivot(int i, int j, Mat3< T > &S, Vec3< T > &D, Mat3< T > &Q)
Definition: Mat3.h:716
Mat3(Source a, Source b, Source c, Source d, Source e, Source f, Source g, Source h, Source i)
Constructor given array of elements, the ordering is in row major form:
Definition: Mat3.h:78
3x3 matrix class.
Definition: Mat3.h:55
const Mat3< T > & operator-=(const Mat3< S > &m1)
Subtract each element of the given matrix from the corresponding element of this matrix.
Definition: Mat3.h:413
const Mat3< T > & operator*=(S scalar)
Multiplication operator, e.g. M = scalar * M;.
Definition: Mat3.h:379
bool operator==(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Equality operator, does exact floating point comparisons.
Definition: Mat3.h:590
Mat3 cofactor() const
Return the cofactor matrix of this matrix.
Definition: Mat3.h:471
void setIdentity()
Set this matrix to identity.
Definition: Mat3.h:324
const Mat3 & operator=(const Mat3< Source > &m)
Assignment operator.
Definition: Mat3.h:339
void setToRotation(const Quat< T > &q)
Set this matrix to the rotation matrix specified by the quaternion.
Definition: Mat3.h:301
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:109
Mat3< typename promote< S, T >::type > operator*(S scalar, const Mat3< T > &m)
Multiply each element of the given matrix by scalar and return the result.
Definition: Mat3.h:609
Vec3< T > row(int i) const
Get ith row, e.g. Vec3d v = m.row(1);.
Definition: Mat3.h:186
void setToRotation(const Vec3< T > &axis, T angle)
Set this matrix to the rotation specified by axis and angle.
Definition: Mat3.h:306
Vec3< typename promote< T, MT >::type > operator*(const Mat3< MT > &_m, const Vec3< T > &_v)
Multiply _m by _v and return the resulting vector.
Definition: Mat3.h:656
Mat3< typename promote< T0, T1 >::type > operator+(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Add corresponding elements of m0 and m1 and return the result.
Definition: Mat3.h:625
Definition: Exceptions.h:83
void setRow(int i, const Vec3< T > &v)
Set ith row to vector v.
Definition: Mat3.h:175
void setCol(int j, const Vec3< T > &v)
Set jth column to vector v.
Definition: Mat3.h:193
T value_type
Definition: Mat.h:56
MatType snapMatBasis(const MatType &source, Axis axis, const Vec3< typename MatType::value_type > &direction)
This function snaps a specific axis to a specific direction, preserving scaling.
Definition: Mat.h:783
static const Mat3< T > & zero()
Predefined constant for zero matrix.
Definition: Mat3.h:165
T operator()(int i, int j) const
Definition: Mat3.h:235
const Mat3< T > & operator+=(const Mat3< S > &m1)
Add each element of the given matrix to the corresponding element of this matrix. ...
Definition: Mat3.h:395
MatType skew(const Vec3< typename MatType::value_type > &skew)
Return a matrix as the cross product of the given vector.
Definition: Mat.h:740
T & operator()(int i, int j)
Definition: Mat3.h:225
Mat3 snapBasis(Axis axis, const Vec3< T > &direction)
Definition: Mat3.h:545
bool eq(const Mat3 &m, T eps=1.0e-8) const
Return true if this matrix is equivalent to m within a tolerance of eps.
Definition: Mat3.h:349
T mm[SIZE *SIZE]
Definition: Mat.h:192
Vec3< T0 > transform(const Vec3< T0 > &v) const
Definition: Mat3.h:553
const T * operator[](int i) const
Definition: Mat3.h:216
T angle(const Vec2< T > &v1, const Vec2< T > &v2)
Definition: Vec2.h:472
bool operator!=(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Inequality operator, does exact floating point comparisons.
Definition: Mat3.h:604
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:128
Mat3(const Mat3< Source > &m)
Conversion constructor.
Definition: Mat3.h:135
Mat3< typename promote< T0, T1 >::type > operator*(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Multiply m0 by m1 and return the resulting matrix.
Definition: Mat3.h:645
Mat3< typename promote< S, T >::type > operator*(const Mat3< T > &m, S scalar)
Multiply each element of the given matrix by scalar and return the result.
Definition: Mat3.h:615
bool isApproxEqual(const Type &a, const Type &b)
Return true if a is equal to b to within the default floating-point comparison tolerance.
Definition: Math.h:378
Definition: Exceptions.h:40
void setColumns(const Vec3< T > &v1, const Vec3< T > &v2, const Vec3< T > &v3)
Set the columns of this matrix to the vectors v1, v2, v3.
Definition: Mat3.h:257
Vec3< typename promote< T, MT >::type > operator*(const Vec3< T > &_v, const Mat3< MT > &_m)
Multiply _v by _m and return the resulting vector.
Definition: Mat3.h:669
Vec3< T > col(int j) const
Get jth column, e.g. Vec3d v = m.col(0);.
Definition: Mat3.h:202
T ValueType
Definition: Mat.h:57
T trace() const
Trace of matrix.
Definition: Mat3.h:536
void setSymmetric(const Vec3< T > &vdiag, const Vec3< T > &vtri)
Set diagonal and symmetric triangular components.
Definition: Mat3.h:271
Tolerance for floating-point comparison.
Definition: Math.h:117
const T * asPointer() const
Definition: Mat3.h:220
Axis
Definition: Math.h:876
void setRows(const Vec3< T > &v1, const Vec3< T > &v2, const Vec3< T > &v3)
Set the rows of this matrix to the vectors v1, v2, v3.
Definition: Mat3.h:243
T * asPointer()
Definition: Mat3.h:219
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:502
Vec3< T0 > pretransform(const Vec3< T0 > &v) const
Definition: Mat3.h:561
Mat3< T > powLerp(const Mat3< T0 > &m1, const Mat3< T0 > &m2, T t)
Definition: Mat3.h:703
Mat3(const Vec3< Source > &v1, const Vec3< Source > &v2, const Vec3< Source > &v3, bool rows=true)
Definition: Mat3.h:96
Mat3< typename promote< T0, T1 >::type > operator-(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Subtract corresponding elements of m0 and m1 and return the result.
Definition: Mat3.h:635
Mat3(const Mat< 3, T > &m)
Copy constructor.
Definition: Mat3.h:124
const Mat3< T > & operator*=(const Mat3< S > &m1)
Multiply this matrix by the given matrix.
Definition: Mat3.h:431
Mat3()
Trivial constructor, the matrix is NOT initialized.
Definition: Mat3.h:63
Mat3(Source *a)
Definition: Mat3.h:110
T det() const
Determinant of matrix.
Definition: Mat3.h:527
Mat3 timesDiagonal(const Vec3< T > &diag) const
Treat diag as a diagonal matrix and return the product of this matrix with diag (from the right)...
Definition: Mat3.h:569
Mat3< T > outerProduct(const Vec3< T > &v1, const Vec3< T > &v2)
Definition: Mat3.h:691
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:180
static Mat3 symmetric(const Vec3< T > &vdiag, const Vec3< T > &vtri)
Return a matrix with the prescribed diagonal and symmetric triangular components. ...
Definition: Mat3.h:285
4x4 -matrix class.
Definition: Mat3.h:49
Mat3 inverse(T tolerance=0) const
Definition: Mat3.h:513
bool diagonalizeSymmetricMatrix(const Mat3< T > &input, Mat3< T > &Q, Vec3< T > &D, unsigned int MAX_ITERATIONS=250)
Use Jacobi iterations to decompose a symmetric 3x3 matrix (diagonalize and compute eigenvectors) ...
Definition: Mat3.h:785
void setZero()
Set this matrix to zero.
Definition: Mat3.h:310
void powSolve(const MatType &aA, MatType &aB, double aPower, double aTol=0.01)
Definition: Mat.h:854
Mat3 adjoint() const
Return the adjoint of this matrix, i.e., the transpose of its cofactor.
Definition: Mat3.h:486
static const Mat3< T > & identity()
Predefined constant for identity matrix.
Definition: Mat3.h:155
Mat3(const Quat< T > &q)
Definition: Mat3.h:67
Mat3(const Mat4< T > &m)
Conversion from Mat4 (copies top left)
Definition: Mat3.h:145
void setSkew(const Vec3< T > &v)
Set the matrix as cross product of the given vector.
Definition: Mat3.h:295
Mat3< double > Mat3d
Definition: Mat3.h:843
Mat3< T > operator-() const
Negation operator, for e.g. m1 = -m2;.
Definition: Mat3.h:363