OpenVDB  2.0.0
Mat3.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2013 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] = a;
82  MyBase::mm[1] = b;
83  MyBase::mm[2] = c;
84  MyBase::mm[3] = d;
85  MyBase::mm[4] = e;
86  MyBase::mm[5] = f;
87  MyBase::mm[6] = g;
88  MyBase::mm[7] = h;
89  MyBase::mm[8] = 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-2013 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/ )
void setCol(int j, const Vec3< T > &v)
Set jth column to vector v.
Definition: Mat3.h:175
Mat3< typename promote< S, T >::type > operator*(const Mat3< T > &m, S scalar)
Returns M, where for .
Definition: Mat3.h:575
const T * asPointer() const
Definition: Mat3.h:202
const Mat3< T > & operator-=(const Mat3< S > &m1)
Returns m0, where for .
Definition: Mat3.h:382
const Mat3< T > & operator*=(const Mat3< S > &m1)
Returns m0, where for .
Definition: Mat3.h:400
void setRow(int i, const Vec3< T > &v)
Set ith row to vector v.
Definition: Mat3.h:157
const Mat3< T > & operator*=(S scalar)
Multiplication operator, e.g. M = scalar * M;.
Definition: Mat3.h:348
Mat3s Mat3f
Definition: Mat3.h:671
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
Mat3< double > Mat3d
Definition: Mat3.h:666
void setBasis(const Vec3< T > &v1, const Vec3< T > &v2, const Vec3< T > &v3)
Set the columns of &quot;this&quot; matrix to the vectors v1, v2, v3.
Definition: Mat3.h:225
Definition: Exceptions.h:78
MatType skew(const Vec3< typename MatType::value_type > &skew)
Definition: Mat.h:689
void setSymmetric(const Vec3< T > &vdiag, const Vec3< T > &vtri)
Set diagonal and symmetric triangular components.
Definition: Mat3.h:239
T angle(const Vec2< T > &v1, const Vec2< T > &v2)
Definition: Vec2.h:431
Tolerance for floating-point comparison.
Definition: Math.h:116
Vec3< T > row(int i) const
Get ith row, e.g. Vec3d v = m.row(1);.
Definition: Mat3.h:168
Mat3()
Trivial constructor, the matrix is NOT initialized.
Definition: Mat3.h:62
const Mat3< T > & operator+=(const Mat3< S > &m1)
Returns m0, where for .
Definition: Mat3.h:364
void setIdentity()
Set &quot;this&quot; matrix to identity.
Definition: Mat3.h:293
Mat3(const Mat< 3, T > &m)
Copy constructor.
Definition: Mat3.h:116
Mat3< typename promote< T0, T1 >::type > operator-(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Returns M, where for .
Definition: Mat3.h:595
4x4 -matrix class.
Definition: Mat3.h:48
T * asPointer()
Definition: Mat3.h:201
Mat3 snapBasis(Axis axis, const Vec3< T > &direction)
Definition: Mat3.h:500
Mat3 adjoint() const
returns adjoint of m
Definition: Mat3.h:440
T det() const
Determinant of matrix.
Definition: Mat3.h:481
T & operator()(int i, int j)
Definition: Mat3.h:207
Mat< 3, T > MyBase
Definition: Mat3.h:60
Vec3< T0 > pretransform(const Vec3< T0 > &v) const
Definition: Mat3.h:516
#define OPENVDB_VERSION_NAME
Definition: version.h:45
const T * operator[](int i) const
Definition: Mat3.h:198
T mm[SIZE *SIZE]
Definition: Mat.h:145
Mat3< T > operator-() const
Negation operator, for e.g. m1 = -m2;.
Definition: Mat3.h:332
Vec3< T > col(int j) const
Get jth column, e.g. Vec3d v = m.col(0);.
Definition: Mat3.h:184
T operator()(int i, int j) const
Definition: Mat3.h:217
T trace() const
Trace of matrix.
Definition: Mat3.h:491
Mat3(const Mat4< T > &m)
Conversion from Mat4 (copies top left)
Definition: Mat3.h:137
Mat3(Source *a)
Definition: Mat3.h:102
void setToRotation(const Vec3< T > &axis, T angle)
Set this matrix to the rotation specified by axis and angle.
Definition: Mat3.h:275
static const Mat3< T > & zero()
Predefined constant for zero matrix.
Definition: Mat3.h:152
Mat3< typename promote< S, T >::type > operator*(S scalar, const Mat3< T > &m)
Returns M, where for .
Definition: Mat3.h:569
Vec3< typename promote< T, MT >::type > operator*(const Mat3< MT > &_m, const Vec3< T > &_v)
Returns v, where for .
Definition: Mat3.h:619
static const Mat3< T > & identity()
Predefined constant for identity matrix.
Definition: Mat3.h:147
3x3 matrix class.
Definition: Mat3.h:54
Mat3< T > powLerp(const Mat3< T0 > &m1, const Mat3< T0 > &m2, T t)
Definition: Mat3.h:679
void powSolve(const MatType &aA, MatType &aB, double aPower, double aTol=0.01)
Definition: Mat.h:779
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:353
Definition: Mat.h:150
Mat3(const Quat< T > &q)
Definition: Mat3.h:66
T ValueType
Definition: Mat3.h:59
T value_type
Data type held by the matrix.
Definition: Mat3.h:58
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:455
Mat3 inverse(T tolerance=0) const
Definition: Mat3.h:466
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
static Mat3 symmetric(const Vec3< T > &vdiag, const Vec3< T > &vtri)
Definition: Mat3.h:254
Mat3< typename promote< T0, T1 >::type > operator+(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Returns M, where for .
Definition: Mat3.h:585
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
Mat3(const Vec3< Source > &v1, const Vec3< Source > &v2, const Vec3< Source > &v3)
Construct matrix given basis vectors (columns)
Definition: Mat3.h:94
void setZero()
Set this matrix to zero.
Definition: Mat3.h:279
Mat3(const Mat3< Source > &m)
Conversion constructor.
Definition: Mat3.h:127
Definition: Mat.h:149
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:56
Axis
Definition: Math.h:762
void setToRotation(const Quat< T > &q)
Set this matrix to the rotation matrix specified by the quaternion.
Definition: Mat3.h:270
bool isApproxEqual(const Hermite &lhs, const Hermite &rhs)
Definition: Hermite.h:470
bool operator==(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Equality operator, does exact floating point comparisons.
Definition: Mat3.h:550
const Mat3 & operator=(const Mat3< Source > &m)
Assignment operator.
Definition: Mat3.h:308
Mat3< float > Mat3s
Definition: Mat3.h:665
Mat3< T > outerProduct(const Vec3< T > &v1, const Vec3< T > &v2)
Definition: Mat3.h:654
bool eq(const Mat3 &m, T eps=1.0e-8) const
Test if &quot;this&quot; is equivalent to m with tolerance of eps value.
Definition: Mat3.h:318
Mat3< typename promote< T0, T1 >::type > operator*(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Matrix multiplication.
Definition: Mat3.h:608
Definition: Mat.h:56
Vec3< T0 > transform(const Vec3< T0 > &v) const
Definition: Mat3.h:508
Mat3 snappedBasis(Axis axis, const Vec3< T0 > &direction) const
Definition: Mat3.h:526
Vec3< typename promote< T, MT >::type > operator*(const Vec3< T > &_v, const Mat3< MT > &_m)
Returns v, where for .
Definition: Mat3.h:632
bool operator!=(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Inequality operator, does exact floating point comparisons.
Definition: Mat3.h:564
void setSkew(const Vec3< T > &v)
Set the matrix as cross product of the given vector.
Definition: Mat3.h:264