OpenVDB  3.0.0
Mat3.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2014 DreamWorks Animation LLC
4 //
5 // All rights reserved. This software is distributed under the
6 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
7 //
8 // Redistributions of source code must retain the above copyright
9 // and license notice and the following restrictions and disclaimer.
10 //
11 // * Neither the name of DreamWorks Animation nor the names of
12 // its contributors may be used to endorse or promote products derived
13 // from this software without specific prior written permission.
14 //
15 // THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
16 // "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
17 // LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
18 // A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
19 // OWNER OR CONTRIBUTORS BE LIABLE FOR ANY INDIRECT, INCIDENTAL,
20 // SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
21 // LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
22 // DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
23 // THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
24 // (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
25 // OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
26 // IN NO EVENT SHALL THE COPYRIGHT HOLDERS' AND CONTRIBUTORS' AGGREGATE
27 // LIABILITY FOR ALL CLAIMS REGARDLESS OF THEIR BASIS EXCEED US$250.00.
28 //
30 
31 #ifndef OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED
32 #define OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED
33 
34 #include <iomanip>
35 #include <assert.h>
36 #include <math.h>
37 #include <openvdb/Exceptions.h>
38 #include "Vec3.h"
39 #include "Mat.h"
40 
41 
42 namespace openvdb {
44 namespace OPENVDB_VERSION_NAME {
45 namespace math {
46 
47 template<typename T> class Vec3;
48 template<typename T> class Mat4;
49 template<typename T> class Quat;
50 
53 template<typename T>
54 class Mat3: public Mat<3, T>
55 {
56 public:
58  typedef T value_type;
59  typedef T ValueType;
60  typedef Mat<3, T> MyBase;
62  Mat3() {}
63 
66  Mat3(const Quat<T> &q)
67  { setToRotation(q); }
68 
69 
71 
76  template<typename Source>
77  Mat3(Source a, Source b, Source c,
78  Source d, Source e, Source f,
79  Source g, Source h, Source i)
80  {
81  MyBase::mm[0] = static_cast<ValueType>(a);
82  MyBase::mm[1] = static_cast<ValueType>(b);
83  MyBase::mm[2] = static_cast<ValueType>(c);
84  MyBase::mm[3] = static_cast<ValueType>(d);
85  MyBase::mm[4] = static_cast<ValueType>(e);
86  MyBase::mm[5] = static_cast<ValueType>(f);
87  MyBase::mm[6] = static_cast<ValueType>(g);
88  MyBase::mm[7] = static_cast<ValueType>(h);
89  MyBase::mm[8] = static_cast<ValueType>(i);
90  } // constructor1Test
91 
93  template<typename Source>
94  Mat3(const Vec3<Source> &v1, const Vec3<Source> &v2, const Vec3<Source> &v3)
95  { setBasis(v1, v2, v3); }
96 
101  template<typename Source>
102  Mat3(Source *a)
103  {
104  MyBase::mm[0] = a[0];
105  MyBase::mm[1] = a[1];
106  MyBase::mm[2] = a[2];
107  MyBase::mm[3] = a[3];
108  MyBase::mm[4] = a[4];
109  MyBase::mm[5] = a[5];
110  MyBase::mm[6] = a[6];
111  MyBase::mm[7] = a[7];
112  MyBase::mm[8] = a[8];
113  } // constructor1Test
114 
116  Mat3(const Mat<3, T> &m)
117  {
118  for (int i=0; i<3; ++i) {
119  for (int j=0; j<3; ++j) {
120  MyBase::mm[i*3 + j] = m[i][j];
121  }
122  }
123  }
124 
126  template<typename Source>
127  explicit Mat3(const Mat3<Source> &m)
128  {
129  for (int i=0; i<3; ++i) {
130  for (int j=0; j<3; ++j) {
131  MyBase::mm[i*3 + j] = m[i][j];
132  }
133  }
134  }
135 
137  explicit Mat3(const Mat4<T> &m)
138  {
139  for (int i=0; i<3; ++i) {
140  for (int j=0; j<3; ++j) {
141  MyBase::mm[i*3 + j] = m[i][j];
142  }
143  }
144  }
145 
147  static const Mat3<T>& identity() {
148  return sIdentity;
149  }
150 
152  static const Mat3<T>& zero() {
153  return sZero;
154  }
155 
157  void setRow(int i, const Vec3<T> &v)
158  {
159  // assert(i>=0 && i<3);
160  int i3 = i * 3;
161 
162  MyBase::mm[i3+0] = v[0];
163  MyBase::mm[i3+1] = v[1];
164  MyBase::mm[i3+2] = v[2];
165  } // rowColumnTest
166 
168  Vec3<T> row(int i) const
169  {
170  // assert(i>=0 && i<3);
171  return Vec3<T>((*this)(i,0), (*this)(i,1), (*this)(i,2));
172  } // rowColumnTest
173 
175  void setCol(int j, const Vec3<T>& v)
176  {
177  // assert(j>=0 && j<3);
178  MyBase::mm[0+j] = v[0];
179  MyBase::mm[3+j] = v[1];
180  MyBase::mm[6+j] = v[2];
181  } // rowColumnTest
182 
184  Vec3<T> col(int j) const
185  {
186  // assert(j>=0 && j<3);
187  return Vec3<T>((*this)(0,j), (*this)(1,j), (*this)(2,j));
188  } // rowColumnTest
189 
190  // NB: The following two methods were changed to
191  // work around a gccWS5 compiler issue related to strict
192  // aliasing (see FX-475).
193 
195  T* operator[](int i) { return &(MyBase::mm[i*3]); }
198  const T* operator[](int i) const { return &(MyBase::mm[i*3]); }
200 
201  T* asPointer() {return MyBase::mm;}
202  const T* asPointer() const {return MyBase::mm;}
203 
207  T& operator()(int i, int j)
208  {
209  // assert(i>=0 && i<3);
210  // assert(j>=0 && j<3);
211  return MyBase::mm[3*i+j];
212  } // trivial
213 
217  T operator()(int i, int j) const
218  {
219  // assert(i>=0 && i<3);
220  // assert(j>=0 && j<3);
221  return MyBase::mm[3*i+j];
222  } // trivial
223 
225  void setBasis(const Vec3<T> &v1, const Vec3<T> &v2, const Vec3<T> &v3)
226  {
227  MyBase::mm[0] = v1[0];
228  MyBase::mm[1] = v1[1];
229  MyBase::mm[2] = v1[2];
230  MyBase::mm[3] = v2[0];
231  MyBase::mm[4] = v2[1];
232  MyBase::mm[5] = v2[2];
233  MyBase::mm[6] = v3[0];
234  MyBase::mm[7] = v3[1];
235  MyBase::mm[8] = v3[2];
236  } // setBasisTest
237 
239  void setSymmetric(const Vec3<T> &vdiag, const Vec3<T> &vtri)
240  {
241  MyBase::mm[0] = vdiag[0];
242  MyBase::mm[1] = vtri[0];
243  MyBase::mm[2] = vtri[1];
244  MyBase::mm[3] = vtri[0];
245  MyBase::mm[4] = vdiag[1];
246  MyBase::mm[5] = vtri[2];
247  MyBase::mm[6] = vtri[1];
248  MyBase::mm[7] = vtri[2];
249  MyBase::mm[8] = vdiag[2];
250  } // setSymmetricTest
251 
254  static Mat3 symmetric(const Vec3<T> &vdiag, const Vec3<T> &vtri)
255  {
256  return Mat3(
257  vdiag[0], vtri[0], vtri[1],
258  vtri[0], vdiag[1], vtri[2],
259  vtri[1], vtri[2], vdiag[2]
260  );
261  }
262 
264  void setSkew(const Vec3<T> &v)
265  {*this = skew(v);}
266 
270  void setToRotation(const Quat<T> &q)
271  {*this = rotation<Mat3<T> >(q);}
272 
275  void setToRotation(const Vec3<T> &axis, T angle)
276  {*this = rotation<Mat3<T> >(axis, angle);}
277 
279  void setZero()
280  {
281  MyBase::mm[0] = 0;
282  MyBase::mm[1] = 0;
283  MyBase::mm[2] = 0;
284  MyBase::mm[3] = 0;
285  MyBase::mm[4] = 0;
286  MyBase::mm[5] = 0;
287  MyBase::mm[6] = 0;
288  MyBase::mm[7] = 0;
289  MyBase::mm[8] = 0;
290  } // trivial
291 
293  void setIdentity()
294  {
295  MyBase::mm[0] = 1;
296  MyBase::mm[1] = 0;
297  MyBase::mm[2] = 0;
298  MyBase::mm[3] = 0;
299  MyBase::mm[4] = 1;
300  MyBase::mm[5] = 0;
301  MyBase::mm[6] = 0;
302  MyBase::mm[7] = 0;
303  MyBase::mm[8] = 1;
304  } // trivial
305 
307  template<typename Source>
308  const Mat3& operator=(const Mat3<Source> &m)
309  {
310  const Source *src = m.asPointer();
311 
312  // don't suppress type conversion warnings
313  std::copy(src, (src + this->numElements()), MyBase::mm);
314  return *this;
315  } // opEqualToTest
316 
318  bool eq(const Mat3 &m, T eps=1.0e-8) const
319  {
320  return (isApproxEqual(MyBase::mm[0],m.mm[0],eps) &&
321  isApproxEqual(MyBase::mm[1],m.mm[1],eps) &&
322  isApproxEqual(MyBase::mm[2],m.mm[2],eps) &&
323  isApproxEqual(MyBase::mm[3],m.mm[3],eps) &&
324  isApproxEqual(MyBase::mm[4],m.mm[4],eps) &&
325  isApproxEqual(MyBase::mm[5],m.mm[5],eps) &&
326  isApproxEqual(MyBase::mm[6],m.mm[6],eps) &&
327  isApproxEqual(MyBase::mm[7],m.mm[7],eps) &&
328  isApproxEqual(MyBase::mm[8],m.mm[8],eps));
329  } // trivial
330 
333  {
334  return Mat3<T>(
335  -MyBase::mm[0], -MyBase::mm[1], -MyBase::mm[2],
336  -MyBase::mm[3], -MyBase::mm[4], -MyBase::mm[5],
337  -MyBase::mm[6], -MyBase::mm[7], -MyBase::mm[8]
338  );
339  } // trivial
340 
342  // friend Mat3 operator*(T scalar, const Mat3& m) {
343  // return m*scalar;
344  // }
345 
347  template <typename S>
348  const Mat3<T>& operator*=(S scalar)
349  {
350  MyBase::mm[0] *= scalar;
351  MyBase::mm[1] *= scalar;
352  MyBase::mm[2] *= scalar;
353  MyBase::mm[3] *= scalar;
354  MyBase::mm[4] *= scalar;
355  MyBase::mm[5] *= scalar;
356  MyBase::mm[6] *= scalar;
357  MyBase::mm[7] *= scalar;
358  MyBase::mm[8] *= scalar;
359  return *this;
360  }
361 
363  template <typename S>
364  const Mat3<T> &operator+=(const Mat3<S> &m1)
365  {
366  const S *s = m1.asPointer();
367 
368  MyBase::mm[0] += s[0];
369  MyBase::mm[1] += s[1];
370  MyBase::mm[2] += s[2];
371  MyBase::mm[3] += s[3];
372  MyBase::mm[4] += s[4];
373  MyBase::mm[5] += s[5];
374  MyBase::mm[6] += s[6];
375  MyBase::mm[7] += s[7];
376  MyBase::mm[8] += s[8];
377  return *this;
378  }
379 
381  template <typename S>
382  const Mat3<T> &operator-=(const Mat3<S> &m1)
383  {
384  const S *s = m1.asPointer();
385 
386  MyBase::mm[0] -= s[0];
387  MyBase::mm[1] -= s[1];
388  MyBase::mm[2] -= s[2];
389  MyBase::mm[3] -= s[3];
390  MyBase::mm[4] -= s[4];
391  MyBase::mm[5] -= s[5];
392  MyBase::mm[6] -= s[6];
393  MyBase::mm[7] -= s[7];
394  MyBase::mm[8] -= s[8];
395  return *this;
396  }
397 
399  template <typename S>
400  const Mat3<T> &operator*=(const Mat3<S> &m1)
401  {
402  Mat3<T> m0(*this);
403  const T* s0 = m0.asPointer();
404  const S* s1 = m1.asPointer();
405 
406  MyBase::mm[0] = static_cast<T>(s0[0] * s1[0] +
407  s0[1] * s1[3] +
408  s0[2] * s1[6]);
409  MyBase::mm[1] = static_cast<T>(s0[0] * s1[1] +
410  s0[1] * s1[4] +
411  s0[2] * s1[7]);
412  MyBase::mm[2] = static_cast<T>(s0[0] * s1[2] +
413  s0[1] * s1[5] +
414  s0[2] * s1[8]);
415 
416  MyBase::mm[3] = static_cast<T>(s0[3] * s1[0] +
417  s0[4] * s1[3] +
418  s0[5] * s1[6]);
419  MyBase::mm[4] = static_cast<T>(s0[3] * s1[1] +
420  s0[4] * s1[4] +
421  s0[5] * s1[7]);
422  MyBase::mm[5] = static_cast<T>(s0[3] * s1[2] +
423  s0[4] * s1[5] +
424  s0[5] * s1[8]);
425 
426  MyBase::mm[6] = static_cast<T>(s0[6] * s1[0] +
427  s0[7] * s1[3] +
428  s0[8] * s1[6]);
429  MyBase::mm[7] = static_cast<T>(s0[6] * s1[1] +
430  s0[7] * s1[4] +
431  s0[8] * s1[7]);
432  MyBase::mm[8] = static_cast<T>(s0[6] * s1[2] +
433  s0[7] * s1[5] +
434  s0[8] * s1[8]);
435 
436  return *this;
437  }
438 
440  Mat3 adjoint() const
441  {
442  return Mat3<T>(
443  MyBase::mm[4] * MyBase::mm[8] - MyBase::mm[5] * MyBase::mm[7],
444  MyBase::mm[2] * MyBase::mm[7] - MyBase::mm[1] * MyBase::mm[8],
445  MyBase::mm[1] * MyBase::mm[5] - MyBase::mm[2] * MyBase::mm[4],
446  MyBase::mm[5] * MyBase::mm[6] - MyBase::mm[3] * MyBase::mm[8],
447  MyBase::mm[0] * MyBase::mm[8] - MyBase::mm[2] * MyBase::mm[6],
448  MyBase::mm[2] * MyBase::mm[3] - MyBase::mm[0] * MyBase::mm[5],
449  MyBase::mm[3] * MyBase::mm[7] - MyBase::mm[4] * MyBase::mm[6],
450  MyBase::mm[1] * MyBase::mm[6] - MyBase::mm[0] * MyBase::mm[7],
451  MyBase::mm[0] * MyBase::mm[4] - MyBase::mm[1] * MyBase::mm[3]);
452  } // adjointTest
453 
455  Mat3 transpose() const
456  {
457  return Mat3<T>(
458  MyBase::mm[0], MyBase::mm[3], MyBase::mm[6],
459  MyBase::mm[1], MyBase::mm[4], MyBase::mm[7],
460  MyBase::mm[2], MyBase::mm[5], MyBase::mm[8]);
461 
462  } // transposeTest
463 
466  Mat3 inverse(T tolerance = 0) const
467  {
468  Mat3<T> inv(adjoint());
469 
470  T det = inv.mm[0]*MyBase::mm[0] + inv.mm[1]*MyBase::mm[3] + inv.mm[2]*MyBase::mm[6];
471 
472  // If the determinant is 0, m was singular and "this" will contain junk.
473  if (isApproxEqual(det,0.0,tolerance))
474  {
475  OPENVDB_THROW(ArithmeticError, "Inversion of singular 3x3 matrix");
476  }
477  return inv * (T(1)/det);
478  } // invertTest
479 
481  T det() const
482  {
483  T co00 = MyBase::mm[4]*MyBase::mm[8] - MyBase::mm[5]*MyBase::mm[7];
484  T co10 = MyBase::mm[5]*MyBase::mm[6] - MyBase::mm[3]*MyBase::mm[8];
485  T co20 = MyBase::mm[3]*MyBase::mm[7] - MyBase::mm[4]*MyBase::mm[6];
486  T d = MyBase::mm[0]*co00 + MyBase::mm[1]*co10 + MyBase::mm[2]*co20;
487  return d;
488  } // determinantTest
489 
491  T trace() const
492  {
493  return MyBase::mm[0]+MyBase::mm[4]+MyBase::mm[8];
494  }
495 
500  Mat3 snapBasis(Axis axis, const Vec3<T> &direction)
501  {
502  return snapBasis(*this, axis, direction);
503  }
504 
507  template<typename T0>
508  Vec3<T0> transform(const Vec3<T0> &v) const
509  {
510  return static_cast< Vec3<T0> >(v * *this);
511  } // xformVectorTest
512 
515  template<typename T0>
517  {
518  return static_cast< Vec3<T0> >(*this * v);
519  } // xformTVectorTest
520 
525  template<typename T0>
526  Mat3 snappedBasis(Axis axis, const Vec3<T0>& direction) const
527  {
528  return snapBasis(*this, axis, direction);
529  }
530 
531 private:
532  static const Mat3<T> sIdentity;
533  static const Mat3<T> sZero;
534 }; // class Mat3
535 
536 
537 template <typename T>
538 const Mat3<T> Mat3<T>::sIdentity = Mat3<T>(1, 0, 0,
539  0, 1, 0,
540  0, 0, 1);
541 
542 template <typename T>
543 const Mat3<T> Mat3<T>::sZero = Mat3<T>(0, 0, 0,
544  0, 0, 0,
545  0, 0, 0);
546 
549 template <typename T0, typename T1>
550 bool operator==(const Mat3<T0> &m0, const Mat3<T1> &m1)
551 {
552  const T0 *t0 = m0.asPointer();
553  const T1 *t1 = m1.asPointer();
554 
555  for (int i=0; i<9; ++i) {
556  if (!isExactlyEqual(t0[i], t1[i])) return false;
557  }
558  return true;
559 }
560 
563 template <typename T0, typename T1>
564 bool operator!=(const Mat3<T0> &m0, const Mat3<T1> &m1) { return !(m0 == m1); }
565 
568 template <typename S, typename T>
570 { return m*scalar; }
571 
574 template <typename S, typename T>
576 {
578  result *= scalar;
579  return result;
580 }
581 
584 template <typename T0, typename T1>
586 {
588  result += m1;
589  return result;
590 }
591 
594 template <typename T0, typename T1>
596 {
598  result -= m1;
599  return result;
600 }
601 
602 
607 template <typename T0, typename T1>
609 {
611  result *= m1;
612  return result;
613 }
614 
617 template<typename T, typename MT>
618 inline Vec3<typename promote<T, MT>::type>
619 operator*(const Mat3<MT> &_m, const Vec3<T> &_v)
620 {
621  MT const *m = _m.asPointer();
623  _v[0]*m[0] + _v[1]*m[1] + _v[2]*m[2],
624  _v[0]*m[3] + _v[1]*m[4] + _v[2]*m[5],
625  _v[0]*m[6] + _v[1]*m[7] + _v[2]*m[8]);
626 }
627 
630 template<typename T, typename MT>
632 operator*(const Vec3<T> &_v, const Mat3<MT> &_m)
633 {
634  MT const *m = _m.asPointer();
636  _v[0]*m[0] + _v[1]*m[3] + _v[2]*m[6],
637  _v[0]*m[1] + _v[1]*m[4] + _v[2]*m[7],
638  _v[0]*m[2] + _v[1]*m[5] + _v[2]*m[8]);
639 }
640 
643 template<typename T, typename MT>
644 inline Vec3<T> &operator *= (Vec3<T> &_v, const Mat3<MT> &_m)
645 {
646  Vec3<T> mult = _v * _m;
647  _v = mult;
648  return _v;
649 }
650 
653 template <typename T>
654 Mat3<T> outerProduct(const Vec3<T>& v1, const Vec3<T>& v2)
655 {
656  Mat3<T> m;
657 
658  m.setBasis(Vec3<T>(v1[0]*v2[0], v1[1]*v2[0], v1[2]*v2[0]),
659  Vec3<T>(v1[0]*v2[1], v1[1]*v2[1], v1[2]*v2[1]),
660  Vec3<T>(v1[0]*v2[2], v1[1]*v2[2], v1[2]*v2[2]));
661 
662  return m;
663 } // outerproductTest
664 
667 
668 #if DWREAL_IS_DOUBLE == 1
669 typedef Mat3d Mat3f;
670 #else
671 typedef Mat3s Mat3f;
672 #endif // DWREAL_IS_DOUBLE
673 
674 
678 template<typename T, typename T0>
679 Mat3<T> powLerp(const Mat3<T0> &m1, const Mat3<T0> &m2, T t)
680 {
681  Mat3<T> x = m1.inverse() * m2;
682  powSolve(x, x, t);
683  Mat3<T> m = m1 * x;
684  return m;
685 }
686 
687 
688 namespace {
689  template<typename T>
690  void pivot(int i, int j, Mat3<T>& S, Vec3<T>& D, Mat3<T>& Q)
691  {
692  const int& n = Mat3<T>::size; // should be 3
693  T temp;
695  double cotan_of_2_theta;
696  double tan_of_theta;
697  double cosin_of_theta;
698  double sin_of_theta;
699  double z;
700 
701  double Sij = S(i,j);
702 
703  double Sjj_minus_Sii = D[j] - D[i];
704 
705  if (fabs(Sjj_minus_Sii) * (10*math::Tolerance<T>::value()) > fabs(Sij)) {
706  tan_of_theta = Sij / Sjj_minus_Sii;
707  } else {
709  cotan_of_2_theta = 0.5*Sjj_minus_Sii / Sij ;
710 
711  if (cotan_of_2_theta < 0.) {
712  tan_of_theta =
713  -1./(sqrt(1. + cotan_of_2_theta*cotan_of_2_theta) - cotan_of_2_theta);
714  } else {
715  tan_of_theta =
716  1./(sqrt(1. + cotan_of_2_theta*cotan_of_2_theta) + cotan_of_2_theta);
717  }
718  }
719 
720  cosin_of_theta = 1./sqrt( 1. + tan_of_theta * tan_of_theta);
721  sin_of_theta = cosin_of_theta * tan_of_theta;
722  z = tan_of_theta * Sij;
723  S(i,j) = 0;
724  D[i] -= z;
725  D[j] += z;
726  for (int k = 0; k < i; ++k) {
727  temp = S(k,i);
728  S(k,i) = cosin_of_theta * temp - sin_of_theta * S(k,j);
729  S(k,j)= sin_of_theta * temp + cosin_of_theta * S(k,j);
730  }
731  for (int k = i+1; k < j; ++k) {
732  temp = S(i,k);
733  S(i,k) = cosin_of_theta * temp - sin_of_theta * S(k,j);
734  S(k,j) = sin_of_theta * temp + cosin_of_theta * S(k,j);
735  }
736  for (int k = j+1; k < n; ++k) {
737  temp = S(i,k);
738  S(i,k) = cosin_of_theta * temp - sin_of_theta * S(j,k);
739  S(j,k) = sin_of_theta * temp + cosin_of_theta * S(j,k);
740  }
741  for (int k = 0; k < n; ++k)
742  {
743  temp = Q(k,i);
744  Q(k,i) = cosin_of_theta * temp - sin_of_theta*Q(k,j);
745  Q(k,j) = sin_of_theta * temp + cosin_of_theta*Q(k,j);
746  }
747  }
748 }
749 
750 
756 template<typename T>
758  unsigned int MAX_ITERATIONS=250)
759 {
762  Q = Mat3<T>::identity();
763  int n = Mat3<T>::size; // should be 3
764 
766  Mat3<T> S(input);
767 
768  for (int i = 0; i < n; ++i) {
769  D[i] = S(i,i);
770  }
771 
772  unsigned int iterations(0);
775  do {
778  double er = 0;
779  for (int i = 0; i < n; ++i) {
780  for (int j = i+1; j < n; ++j) {
781  er += fabs(S(i,j));
782  }
783  }
784  if (std::abs(er) < math::Tolerance<T>::value()) {
785  return true;
786  }
787  iterations++;
788 
789  T max_element = 0;
790  int ip = 0;
791  int jp = 0;
793  for (int i = 0; i < n; ++i) {
794  for (int j = i+1; j < n; ++j){
795 
796  if ( fabs(D[i]) * (10*math::Tolerance<T>::value()) > fabs(S(i,j))) {
798  S(i,j) = 0;
799  }
800  if (fabs(S(i,j)) > max_element) {
801  max_element = fabs(S(i,j));
802  ip = i;
803  jp = j;
804  }
805  }
806  }
807  pivot(ip, jp, S, D, Q);
808  } while (iterations < MAX_ITERATIONS);
809 
810  return false;
811 }
812 
813 } // namespace math
814 } // namespace OPENVDB_VERSION_NAME
815 } // namespace openvdb
816 
817 #endif // OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED
818 
819 // Copyright (c) 2012-2014 DreamWorks Animation LLC
820 // All rights reserved. This software is distributed under the
821 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
const Mat3< T > & operator-=(const Mat3< S > &m1)
Returns m0, where for .
Definition: Mat3.h:382
Vec3< typename promote< T, MT >::type > operator*(const Vec3< T > &_v, const Mat3< MT > &_m)
Returns v, where for .
Definition: Mat3.h:632
Mat3< typename promote< S, T >::type > operator*(const Mat3< T > &m, S scalar)
Returns M, where for .
Definition: Mat3.h:575
void setCol(int j, const Vec3< T > &v)
Set jth column to vector v.
Definition: Mat3.h:175
void setSkew(const Vec3< T > &v)
Set the matrix as cross product of the given vector.
Definition: Mat3.h:264
Definition: Mat.h:146
const Mat3< T > & operator*=(const Mat3< S > &m1)
Returns m0, where for .
Definition: Mat3.h:400
const Mat3< T > & operator*=(S scalar)
Multiplication operator, e.g. M = scalar * M;.
Definition: Mat3.h:348
void setRow(int i, const Vec3< T > &v)
Set ith row to vector v.
Definition: Mat3.h:157
Mat3< float > Mat3s
Definition: Mat3.h:665
const T * asPointer() const
Definition: Mat3.h:202
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
void setSymmetric(const Vec3< T > &vdiag, const Vec3< T > &vtri)
Set diagonal and symmetric triangular components.
Definition: Mat3.h:239
Mat3< typename promote< T0, T1 >::type > operator*(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Matrix multiplication.
Definition: Mat3.h:608
Mat3()
Trivial constructor, the matrix is NOT initialized.
Definition: Mat3.h:62
4x4 -matrix class.
Definition: Mat3.h:48
Definition: Mat.h:145
void setBasis(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:225
Mat3< typename promote< T0, T1 >::type > operator-(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Returns M, where for .
Definition: Mat3.h:595
const Mat3< T > & operator+=(const Mat3< S > &m1)
Returns m0, where for .
Definition: Mat3.h:364
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:757
Mat3(const Mat< 3, T > &m)
Copy constructor.
Definition: Mat3.h:116
Mat3 adjoint() const
returns adjoint of m
Definition: Mat3.h:440
T & operator()(int i, int j)
Definition: Mat3.h:207
MatType skew(const Vec3< typename MatType::value_type > &skew)
Return a matrix as the cross product of the given vector.
Definition: Mat.h:687
Vec3< T > row(int i) const
Get ith row, e.g. Vec3d v = m.row(1);.
Definition: Mat3.h:168
Vec3< T0 > pretransform(const Vec3< T0 > &v) const
Definition: Mat3.h:516
Mat< 3, T > MyBase
Definition: Mat3.h:60
void setIdentity()
Set "this" matrix to identity.
Definition: Mat3.h:293
Definition: Mat.h:52
#define OPENVDB_VERSION_NAME
Definition: version.h:43
static const Mat3< T > & zero()
Predefined constant for zero matrix.
Definition: Mat3.h:152
3x3 matrix class.
Definition: Mat3.h:54
T trace() const
Trace of matrix.
Definition: Mat3.h:491
T * asPointer()
Definition: Mat3.h:201
Mat3 snapBasis(Axis axis, const Vec3< T > &direction)
Definition: Mat3.h:500
T det() const
Determinant of matrix.
Definition: Mat3.h:481
T operator()(int i, int j) const
Definition: Mat3.h:217
Mat3< T > operator-() const
Negation operator, for e.g. m1 = -m2;.
Definition: Mat3.h:332
const T * operator[](int i) const
Definition: Mat3.h:198
Mat3(Source *a)
Definition: Mat3.h:102
Vec3< T > col(int j) const
Get jth column, e.g. Vec3d v = m.col(0);.
Definition: Mat3.h:184
Mat3(const Mat4< T > &m)
Conversion from Mat4 (copies top left)
Definition: Mat3.h:137
Definition: Exceptions.h:39
Mat3< typename promote< S, T >::type > operator*(S scalar, const Mat3< T > &m)
Returns M, where for .
Definition: Mat3.h:569
void setToRotation(const Vec3< T > &axis, T angle)
Set this matrix to the rotation specified by axis and angle.
Definition: Mat3.h:275
bool isApproxEqual(const Hermite &lhs, const Hermite &rhs)
Definition: Hermite.h:470
static const Mat3< T > & identity()
Predefined constant for identity matrix.
Definition: Mat3.h:147
Vec3< typename promote< T, MT >::type > operator*(const Mat3< MT > &_m, const Vec3< T > &_v)
Returns v, where for .
Definition: Mat3.h:619
T ValueType
Definition: Mat3.h:59
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:455
Mat3< T > outerProduct(const Vec3< T > &v1, const Vec3< T > &v2)
Definition: Mat3.h:654
Mat3< double > Mat3d
Definition: Mat3.h:666
Definition: Exceptions.h:78
Axis
Definition: Math.h:822
Mat3< T > powLerp(const Mat3< T0 > &m1, const Mat3< T0 > &m2, T t)
Definition: Mat3.h:679
void setZero()
Set this matrix to zero.
Definition: Mat3.h:279
Mat3(const Quat< T > &q)
Definition: Mat3.h:66
Mat3< typename promote< T0, T1 >::type > operator+(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Returns M, where for .
Definition: Mat3.h:585
T value_type
Data type held by the matrix.
Definition: Mat3.h:58
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:77
static Mat3 symmetric(const Vec3< T > &vdiag, const Vec3< T > &vtri)
Definition: Mat3.h:254
bool operator==(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Equality operator, does exact floating point comparisons.
Definition: Mat3.h:550
Mat3 inverse(T tolerance=0) const
Definition: Mat3.h:466
void setToRotation(const Quat< T > &q)
Set this matrix to the rotation matrix specified by the quaternion.
Definition: Mat3.h:270
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:391
Mat3(const Vec3< Source > &v1, const Vec3< Source > &v2, const Vec3< Source > &v3)
Construct matrix given basis vectors (columns)
Definition: Mat3.h:94
void powSolve(const MatType &aA, MatType &aB, double aPower, double aTol=0.01)
Definition: Mat.h:776
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Mat3s Mat3f
Definition: Mat3.h:671
Vec3< T0 > transform(const Vec3< T0 > &v) const
Definition: Mat3.h:508
Mat3(const Mat3< Source > &m)
Conversion constructor.
Definition: Mat3.h:127
const Mat3 & operator=(const Mat3< Source > &m)
Assignment operator.
Definition: Mat3.h:308
bool eq(const Mat3 &m, T eps=1.0e-8) const
Test if "this" is equivalent to m with tolerance of eps value.
Definition: Mat3.h:318
T angle(const Vec2< T > &v1, const Vec2< T > &v2)
Definition: Vec2.h:446
Mat3 snappedBasis(Axis axis, const Vec3< T0 > &direction) const
Definition: Mat3.h:526
bool operator!=(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Inequality operator, does exact floating point comparisons.
Definition: Mat3.h:564
Tolerance for floating-point comparison.
Definition: Math.h:116
T mm[SIZE *SIZE]
Definition: Mat.h:141