HDK
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
Mat3.h
Go to the documentation of this file.
1 // Copyright Contributors to the OpenVDB Project
2 // SPDX-License-Identifier: MPL-2.0
3 
4 #ifndef OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED
5 #define OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED
6 
7 #include <openvdb/Exceptions.h>
8 #include "Vec3.h"
9 #include "Mat.h"
10 #include <algorithm> // for std::copy()
11 #include <cassert>
12 #include <cmath>
13 #include <iomanip>
14 
15 
16 namespace openvdb {
18 namespace OPENVDB_VERSION_NAME {
19 namespace math {
20 
21 template<typename T> class Vec3;
22 template<typename T> class Mat4;
23 template<typename T> class Quat;
24 
25 /// @class Mat3 Mat3.h
26 /// @brief 3x3 matrix class.
27 template<typename T>
28 class Mat3: public Mat<3, T>
29 {
30 public:
31  /// Data type held by the matrix.
32  using value_type = T;
33  using ValueType = T;
34  using MyBase = Mat<3, T>;
35  /// Trivial constructor, the matrix is NOT initialized
36  Mat3() {}
37 
38  /// Constructor given the quaternion rotation, e.g. Mat3f m(q);
39  /// The quaternion is normalized and used to construct the matrix
40  Mat3(const Quat<T> &q)
41  { setToRotation(q); }
42 
43 
44  /// Constructor given array of elements, the ordering is in row major form:
45  /** @verbatim
46  a b c
47  d e f
48  g h i
49  @endverbatim */
50  template<typename Source>
51  Mat3(Source a, Source b, Source c,
52  Source d, Source e, Source f,
53  Source g, Source h, Source i)
54  {
55  MyBase::mm[0] = static_cast<T>(a);
56  MyBase::mm[1] = static_cast<T>(b);
57  MyBase::mm[2] = static_cast<T>(c);
58  MyBase::mm[3] = static_cast<T>(d);
59  MyBase::mm[4] = static_cast<T>(e);
60  MyBase::mm[5] = static_cast<T>(f);
61  MyBase::mm[6] = static_cast<T>(g);
62  MyBase::mm[7] = static_cast<T>(h);
63  MyBase::mm[8] = static_cast<T>(i);
64  } // constructor1Test
65 
66  /// Construct matrix from rows or columns vectors (defaults to rows
67  /// for historical reasons)
68  template<typename Source>
69  Mat3(const Vec3<Source> &v1, const Vec3<Source> &v2, const Vec3<Source> &v3, bool rows = true)
70  {
71  if (rows) {
72  this->setRows(v1, v2, v3);
73  } else {
74  this->setColumns(v1, v2, v3);
75  }
76  }
77 
78  /// Constructor given array of elements, the ordering is in row major form:\n
79  /// a[0] a[1] a[2]\n
80  /// a[3] a[4] a[5]\n
81  /// a[6] a[7] a[8]\n
82  template<typename Source>
83  Mat3(Source *a)
84  {
85  MyBase::mm[0] = static_cast<T>(a[0]);
86  MyBase::mm[1] = static_cast<T>(a[1]);
87  MyBase::mm[2] = static_cast<T>(a[2]);
88  MyBase::mm[3] = static_cast<T>(a[3]);
89  MyBase::mm[4] = static_cast<T>(a[4]);
90  MyBase::mm[5] = static_cast<T>(a[5]);
91  MyBase::mm[6] = static_cast<T>(a[6]);
92  MyBase::mm[7] = static_cast<T>(a[7]);
93  MyBase::mm[8] = static_cast<T>(a[8]);
94  } // constructor1Test
95 
96  /// Copy constructor
97  Mat3(const Mat<3, T> &m)
98  {
99  for (int i=0; i<3; ++i) {
100  for (int j=0; j<3; ++j) {
101  MyBase::mm[i*3 + j] = m[i][j];
102  }
103  }
104  }
105 
106  /// Conversion constructor
107  template<typename Source>
108  explicit Mat3(const Mat3<Source> &m)
109  {
110  for (int i=0; i<3; ++i) {
111  for (int j=0; j<3; ++j) {
112  MyBase::mm[i*3 + j] = static_cast<T>(m[i][j]);
113  }
114  }
115  }
116 
117  /// Conversion from Mat4 (copies top left)
118  explicit Mat3(const Mat4<T> &m)
119  {
120  for (int i=0; i<3; ++i) {
121  for (int j=0; j<3; ++j) {
122  MyBase::mm[i*3 + j] = m[i][j];
123  }
124  }
125  }
126 
127  /// Predefined constant for identity matrix
128  static const Mat3<T>& identity() {
129  static const Mat3<T> sIdentity = Mat3<T>(
130  1, 0, 0,
131  0, 1, 0,
132  0, 0, 1
133  );
134  return sIdentity;
135  }
136 
137  /// Predefined constant for zero matrix
138  static const Mat3<T>& zero() {
139  static const Mat3<T> sZero = Mat3<T>(
140  0, 0, 0,
141  0, 0, 0,
142  0, 0, 0
143  );
144  return sZero;
145  }
146 
147  /// Set ith row to vector v
148  void setRow(int i, const Vec3<T> &v)
149  {
150  // assert(i>=0 && i<3);
151  int i3 = i * 3;
152 
153  MyBase::mm[i3+0] = v[0];
154  MyBase::mm[i3+1] = v[1];
155  MyBase::mm[i3+2] = v[2];
156  } // rowColumnTest
157 
158  /// Get ith row, e.g. Vec3d v = m.row(1);
159  Vec3<T> row(int i) const
160  {
161  // assert(i>=0 && i<3);
162  return Vec3<T>((*this)(i,0), (*this)(i,1), (*this)(i,2));
163  } // rowColumnTest
164 
165  /// Set jth column to vector v
166  void setCol(int j, const Vec3<T>& v)
167  {
168  // assert(j>=0 && j<3);
169  MyBase::mm[0+j] = v[0];
170  MyBase::mm[3+j] = v[1];
171  MyBase::mm[6+j] = v[2];
172  } // rowColumnTest
173 
174  /// Get jth column, e.g. Vec3d v = m.col(0);
175  Vec3<T> col(int j) const
176  {
177  // assert(j>=0 && j<3);
178  return Vec3<T>((*this)(0,j), (*this)(1,j), (*this)(2,j));
179  } // rowColumnTest
180 
181  /// Alternative indexed reference to the elements
182  /// Note that the indices are row first and column second.
183  /// e.g. m(0,0) = 1;
184  T& operator()(int i, int j)
185  {
186  // assert(i>=0 && i<3);
187  // assert(j>=0 && j<3);
188  return MyBase::mm[3*i+j];
189  } // trivial
190 
191  /// Alternative indexed constant reference to the elements,
192  /// Note that the indices are row first and column second.
193  /// e.g. float f = m(1,0);
194  T operator()(int i, int j) const
195  {
196  // assert(i>=0 && i<3);
197  // assert(j>=0 && j<3);
198  return MyBase::mm[3*i+j];
199  } // trivial
200 
201  /// Set the rows of this matrix to the vectors v1, v2, v3
202  void setRows(const Vec3<T> &v1, const Vec3<T> &v2, const Vec3<T> &v3)
203  {
204  MyBase::mm[0] = v1[0];
205  MyBase::mm[1] = v1[1];
206  MyBase::mm[2] = v1[2];
207  MyBase::mm[3] = v2[0];
208  MyBase::mm[4] = v2[1];
209  MyBase::mm[5] = v2[2];
210  MyBase::mm[6] = v3[0];
211  MyBase::mm[7] = v3[1];
212  MyBase::mm[8] = v3[2];
213  } // setRows
214 
215  /// Set the columns of this matrix to the vectors v1, v2, v3
216  void setColumns(const Vec3<T> &v1, const Vec3<T> &v2, const Vec3<T> &v3)
217  {
218  MyBase::mm[0] = v1[0];
219  MyBase::mm[1] = v2[0];
220  MyBase::mm[2] = v3[0];
221  MyBase::mm[3] = v1[1];
222  MyBase::mm[4] = v2[1];
223  MyBase::mm[5] = v3[1];
224  MyBase::mm[6] = v1[2];
225  MyBase::mm[7] = v2[2];
226  MyBase::mm[8] = v3[2];
227  } // setColumns
228 
229  /// Set diagonal and symmetric triangular components
230  void setSymmetric(const Vec3<T> &vdiag, const Vec3<T> &vtri)
231  {
232  MyBase::mm[0] = vdiag[0];
233  MyBase::mm[1] = vtri[0];
234  MyBase::mm[2] = vtri[1];
235  MyBase::mm[3] = vtri[0];
236  MyBase::mm[4] = vdiag[1];
237  MyBase::mm[5] = vtri[2];
238  MyBase::mm[6] = vtri[1];
239  MyBase::mm[7] = vtri[2];
240  MyBase::mm[8] = vdiag[2];
241  } // setSymmetricTest
242 
243  /// Return a matrix with the prescribed diagonal and symmetric triangular components.
244  static Mat3 symmetric(const Vec3<T> &vdiag, const Vec3<T> &vtri)
245  {
246  return Mat3(
247  vdiag[0], vtri[0], vtri[1],
248  vtri[0], vdiag[1], vtri[2],
249  vtri[1], vtri[2], vdiag[2]
250  );
251  }
252 
253  /// Set the matrix as cross product of the given vector
254  void setSkew(const Vec3<T> &v)
255  {*this = skew(v);}
256 
257  /// @brief Set this matrix to the rotation matrix specified by the quaternion
258  /// @details The quaternion is normalized and used to construct the matrix.
259  /// Note that the matrix is transposed to match post-multiplication semantics.
260  void setToRotation(const Quat<T> &q)
261  {*this = rotation<Mat3<T> >(q);}
262 
263  /// @brief Set this matrix to the rotation specified by @a axis and @a angle
264  /// @details The axis must be unit vector
265  void setToRotation(const Vec3<T> &axis, T angle)
266  {*this = rotation<Mat3<T> >(axis, angle);}
267 
268  /// Set this matrix to zero
269  void setZero()
270  {
271  MyBase::mm[0] = 0;
272  MyBase::mm[1] = 0;
273  MyBase::mm[2] = 0;
274  MyBase::mm[3] = 0;
275  MyBase::mm[4] = 0;
276  MyBase::mm[5] = 0;
277  MyBase::mm[6] = 0;
278  MyBase::mm[7] = 0;
279  MyBase::mm[8] = 0;
280  } // trivial
281 
282  /// Set this matrix to identity
283  void setIdentity()
284  {
285  MyBase::mm[0] = 1;
286  MyBase::mm[1] = 0;
287  MyBase::mm[2] = 0;
288  MyBase::mm[3] = 0;
289  MyBase::mm[4] = 1;
290  MyBase::mm[5] = 0;
291  MyBase::mm[6] = 0;
292  MyBase::mm[7] = 0;
293  MyBase::mm[8] = 1;
294  } // trivial
295 
296  /// Assignment operator
297  template<typename Source>
298  const Mat3& operator=(const Mat3<Source> &m)
299  {
300  const Source *src = m.asPointer();
301 
302  // don't suppress type conversion warnings
303  std::copy(src, (src + this->numElements()), MyBase::mm);
304  return *this;
305  } // opEqualToTest
306 
307  /// Return @c true if this matrix is equivalent to @a m within a tolerance of @a eps.
308  bool eq(const Mat3 &m, T eps=1.0e-8) const
309  {
310  return (isApproxEqual(MyBase::mm[0],m.mm[0],eps) &&
311  isApproxEqual(MyBase::mm[1],m.mm[1],eps) &&
312  isApproxEqual(MyBase::mm[2],m.mm[2],eps) &&
313  isApproxEqual(MyBase::mm[3],m.mm[3],eps) &&
314  isApproxEqual(MyBase::mm[4],m.mm[4],eps) &&
315  isApproxEqual(MyBase::mm[5],m.mm[5],eps) &&
316  isApproxEqual(MyBase::mm[6],m.mm[6],eps) &&
317  isApproxEqual(MyBase::mm[7],m.mm[7],eps) &&
318  isApproxEqual(MyBase::mm[8],m.mm[8],eps));
319  } // trivial
320 
321  /// Negation operator, for e.g. m1 = -m2;
323  {
324  return Mat3<T>(
325  -MyBase::mm[0], -MyBase::mm[1], -MyBase::mm[2],
326  -MyBase::mm[3], -MyBase::mm[4], -MyBase::mm[5],
327  -MyBase::mm[6], -MyBase::mm[7], -MyBase::mm[8]
328  );
329  } // trivial
330 
331  /// Multiplication operator, e.g. M = scalar * M;
332  // friend Mat3 operator*(T scalar, const Mat3& m) {
333  // return m*scalar;
334  // }
335 
336  /// Multiply each element of this matrix by @a scalar.
337  template <typename S>
338  const Mat3<T>& operator*=(S scalar)
339  {
340  MyBase::mm[0] *= scalar;
341  MyBase::mm[1] *= scalar;
342  MyBase::mm[2] *= scalar;
343  MyBase::mm[3] *= scalar;
344  MyBase::mm[4] *= scalar;
345  MyBase::mm[5] *= scalar;
346  MyBase::mm[6] *= scalar;
347  MyBase::mm[7] *= scalar;
348  MyBase::mm[8] *= scalar;
349  return *this;
350  }
351 
352  /// Add each element of the given matrix to the corresponding element of this matrix.
353  template <typename S>
354  const Mat3<T> &operator+=(const Mat3<S> &m1)
355  {
356  const S *s = m1.asPointer();
357 
358  MyBase::mm[0] += s[0];
359  MyBase::mm[1] += s[1];
360  MyBase::mm[2] += s[2];
361  MyBase::mm[3] += s[3];
362  MyBase::mm[4] += s[4];
363  MyBase::mm[5] += s[5];
364  MyBase::mm[6] += s[6];
365  MyBase::mm[7] += s[7];
366  MyBase::mm[8] += s[8];
367  return *this;
368  }
369 
370  /// Subtract each element of the given matrix from the corresponding element of this matrix.
371  template <typename S>
372  const Mat3<T> &operator-=(const Mat3<S> &m1)
373  {
374  const S *s = m1.asPointer();
375 
376  MyBase::mm[0] -= s[0];
377  MyBase::mm[1] -= s[1];
378  MyBase::mm[2] -= s[2];
379  MyBase::mm[3] -= s[3];
380  MyBase::mm[4] -= s[4];
381  MyBase::mm[5] -= s[5];
382  MyBase::mm[6] -= s[6];
383  MyBase::mm[7] -= s[7];
384  MyBase::mm[8] -= s[8];
385  return *this;
386  }
387 
388  /// Multiply this matrix by the given matrix.
389  template <typename S>
390  const Mat3<T> &operator*=(const Mat3<S> &m1)
391  {
392  Mat3<T> m0(*this);
393  const T* s0 = m0.asPointer();
394  const S* s1 = m1.asPointer();
395 
396  MyBase::mm[0] = static_cast<T>(s0[0] * s1[0] +
397  s0[1] * s1[3] +
398  s0[2] * s1[6]);
399  MyBase::mm[1] = static_cast<T>(s0[0] * s1[1] +
400  s0[1] * s1[4] +
401  s0[2] * s1[7]);
402  MyBase::mm[2] = static_cast<T>(s0[0] * s1[2] +
403  s0[1] * s1[5] +
404  s0[2] * s1[8]);
405 
406  MyBase::mm[3] = static_cast<T>(s0[3] * s1[0] +
407  s0[4] * s1[3] +
408  s0[5] * s1[6]);
409  MyBase::mm[4] = static_cast<T>(s0[3] * s1[1] +
410  s0[4] * s1[4] +
411  s0[5] * s1[7]);
412  MyBase::mm[5] = static_cast<T>(s0[3] * s1[2] +
413  s0[4] * s1[5] +
414  s0[5] * s1[8]);
415 
416  MyBase::mm[6] = static_cast<T>(s0[6] * s1[0] +
417  s0[7] * s1[3] +
418  s0[8] * s1[6]);
419  MyBase::mm[7] = static_cast<T>(s0[6] * s1[1] +
420  s0[7] * s1[4] +
421  s0[8] * s1[7]);
422  MyBase::mm[8] = static_cast<T>(s0[6] * s1[2] +
423  s0[7] * s1[5] +
424  s0[8] * s1[8]);
425 
426  return *this;
427  }
428 
429  /// @brief Return the cofactor matrix of this matrix.
430  Mat3 cofactor() const
431  {
432  return Mat3<T>(
433  MyBase::mm[4] * MyBase::mm[8] - MyBase::mm[5] * MyBase::mm[7],
434  MyBase::mm[5] * MyBase::mm[6] - MyBase::mm[3] * MyBase::mm[8],
435  MyBase::mm[3] * MyBase::mm[7] - MyBase::mm[4] * MyBase::mm[6],
436  MyBase::mm[2] * MyBase::mm[7] - MyBase::mm[1] * MyBase::mm[8],
437  MyBase::mm[0] * MyBase::mm[8] - MyBase::mm[2] * MyBase::mm[6],
438  MyBase::mm[1] * MyBase::mm[6] - MyBase::mm[0] * MyBase::mm[7],
439  MyBase::mm[1] * MyBase::mm[5] - MyBase::mm[2] * MyBase::mm[4],
440  MyBase::mm[2] * MyBase::mm[3] - MyBase::mm[0] * MyBase::mm[5],
441  MyBase::mm[0] * MyBase::mm[4] - MyBase::mm[1] * MyBase::mm[3]);
442  }
443 
444  /// Return the adjoint of this matrix, i.e., the transpose of its cofactor.
445  Mat3 adjoint() const
446  {
447  return Mat3<T>(
448  MyBase::mm[4] * MyBase::mm[8] - MyBase::mm[5] * MyBase::mm[7],
449  MyBase::mm[2] * MyBase::mm[7] - MyBase::mm[1] * MyBase::mm[8],
450  MyBase::mm[1] * MyBase::mm[5] - MyBase::mm[2] * MyBase::mm[4],
451  MyBase::mm[5] * MyBase::mm[6] - MyBase::mm[3] * MyBase::mm[8],
452  MyBase::mm[0] * MyBase::mm[8] - MyBase::mm[2] * MyBase::mm[6],
453  MyBase::mm[2] * MyBase::mm[3] - MyBase::mm[0] * MyBase::mm[5],
454  MyBase::mm[3] * MyBase::mm[7] - MyBase::mm[4] * MyBase::mm[6],
455  MyBase::mm[1] * MyBase::mm[6] - MyBase::mm[0] * MyBase::mm[7],
456  MyBase::mm[0] * MyBase::mm[4] - MyBase::mm[1] * MyBase::mm[3]);
457 
458  } // adjointTest
459 
460  /// returns transpose of this
461  Mat3 transpose() const
462  {
463  return Mat3<T>(
464  MyBase::mm[0], MyBase::mm[3], MyBase::mm[6],
465  MyBase::mm[1], MyBase::mm[4], MyBase::mm[7],
466  MyBase::mm[2], MyBase::mm[5], MyBase::mm[8]);
467 
468  } // transposeTest
469 
470  /// returns inverse of this
471  /// @throws ArithmeticError if singular
472  Mat3 inverse(T tolerance = 0) const
473  {
474  Mat3<T> inv(this->adjoint());
475 
476  const T det = inv.mm[0]*MyBase::mm[0] + inv.mm[1]*MyBase::mm[3] + inv.mm[2]*MyBase::mm[6];
477 
478  // If the determinant is 0, m was singular and the result will be invalid.
479  if (isApproxEqual(det,T(0.0),tolerance)) {
480  OPENVDB_THROW(ArithmeticError, "Inversion of singular 3x3 matrix");
481  }
482  return inv * (T(1)/det);
483  } // invertTest
484 
485  /// Determinant of matrix
486  T det() const
487  {
488  const T co00 = MyBase::mm[4]*MyBase::mm[8] - MyBase::mm[5]*MyBase::mm[7];
489  const T co10 = MyBase::mm[5]*MyBase::mm[6] - MyBase::mm[3]*MyBase::mm[8];
490  const T co20 = MyBase::mm[3]*MyBase::mm[7] - MyBase::mm[4]*MyBase::mm[6];
491  return MyBase::mm[0]*co00 + MyBase::mm[1]*co10 + MyBase::mm[2]*co20;
492  } // determinantTest
493 
494  /// Trace of matrix
495  T trace() const
496  {
497  return MyBase::mm[0]+MyBase::mm[4]+MyBase::mm[8];
498  }
499 
500  /// This function snaps a specific axis to a specific direction,
501  /// preserving scaling. It does this using minimum energy, thus
502  /// posing a unique solution if basis & direction arent parralel.
503  /// Direction need not be unit.
505  {
506  return snapMatBasis(*this, axis, direction);
507  }
508 
509  /// Return the transformed vector by this matrix.
510  /// This function is equivalent to post-multiplying the matrix.
511  template<typename T0>
512  Vec3<T0> transform(const Vec3<T0> &v) const
513  {
514  return static_cast< Vec3<T0> >(v * *this);
515  } // xformVectorTest
516 
517  /// Return the transformed vector by transpose of this matrix.
518  /// This function is equivalent to pre-multiplying the matrix.
519  template<typename T0>
521  {
522  return static_cast< Vec3<T0> >(*this * v);
523  } // xformTVectorTest
524 
525 
526  /// @brief Treat @a diag as a diagonal matrix and return the product
527  /// of this matrix with @a diag (from the right).
528  Mat3 timesDiagonal(const Vec3<T>& diag) const
529  {
530  Mat3 ret(*this);
531 
532  ret.mm[0] *= diag(0);
533  ret.mm[1] *= diag(1);
534  ret.mm[2] *= diag(2);
535  ret.mm[3] *= diag(0);
536  ret.mm[4] *= diag(1);
537  ret.mm[5] *= diag(2);
538  ret.mm[6] *= diag(0);
539  ret.mm[7] *= diag(1);
540  ret.mm[8] *= diag(2);
541  return ret;
542  }
543 }; // class Mat3
544 
545 
546 /// @relates Mat3
547 /// @brief Equality operator, does exact floating point comparisons
548 template <typename T0, typename T1>
549 bool operator==(const Mat3<T0> &m0, const Mat3<T1> &m1)
550 {
551  const T0 *t0 = m0.asPointer();
552  const T1 *t1 = m1.asPointer();
553 
554  for (int i=0; i<9; ++i) {
555  if (!isExactlyEqual(t0[i], t1[i])) return false;
556  }
557  return true;
558 }
559 
560 /// @relates Mat3
561 /// @brief Inequality operator, does exact floating point comparisons
562 template <typename T0, typename T1>
563 bool operator!=(const Mat3<T0> &m0, const Mat3<T1> &m1) { return !(m0 == m1); }
564 
565 /// @relates Mat3
566 /// @brief Multiply each element of the given matrix by @a scalar and return the result.
567 template <typename S, typename T>
569 { return m*scalar; }
570 
571 /// @relates Mat3
572 /// @brief Multiply each element of the given matrix by @a scalar and return the result.
573 template <typename S, typename T>
575 {
577  result *= scalar;
578  return result;
579 }
580 
581 /// @relates Mat3
582 /// @brief Add corresponding elements of @a m0 and @a m1 and return the result.
583 template <typename T0, typename T1>
585 {
587  result += m1;
588  return result;
589 }
590 
591 /// @relates Mat3
592 /// @brief Subtract corresponding elements of @a m0 and @a m1 and return the result.
593 template <typename T0, typename T1>
595 {
597  result -= m1;
598  return result;
599 }
600 
601 
602 /// @brief Multiply @a m0 by @a m1 and return the resulting matrix.
603 template <typename T0, typename T1>
605 {
607  result *= m1;
608  return result;
609 }
610 
611 /// @relates Mat3
612 /// @brief Multiply @a _m by @a _v and return the resulting vector.
613 template<typename T, typename MT>
615 operator*(const Mat3<MT> &_m, const Vec3<T> &_v)
616 {
617  MT const *m = _m.asPointer();
619  _v[0]*m[0] + _v[1]*m[1] + _v[2]*m[2],
620  _v[0]*m[3] + _v[1]*m[4] + _v[2]*m[5],
621  _v[0]*m[6] + _v[1]*m[7] + _v[2]*m[8]);
622 }
623 
624 /// @relates Mat3
625 /// @brief Multiply @a _v by @a _m and return the resulting vector.
626 template<typename T, typename MT>
628 operator*(const Vec3<T> &_v, const Mat3<MT> &_m)
629 {
630  MT const *m = _m.asPointer();
632  _v[0]*m[0] + _v[1]*m[3] + _v[2]*m[6],
633  _v[0]*m[1] + _v[1]*m[4] + _v[2]*m[7],
634  _v[0]*m[2] + _v[1]*m[5] + _v[2]*m[8]);
635 }
636 
637 /// @relates Mat3
638 /// @brief Multiply @a _v by @a _m and replace @a _v with the resulting vector.
639 template<typename T, typename MT>
640 inline Vec3<T> &operator *= (Vec3<T> &_v, const Mat3<MT> &_m)
641 {
642  Vec3<T> mult = _v * _m;
643  _v = mult;
644  return _v;
645 }
646 
647 /// Returns outer product of v1, v2, i.e. v1 v2^T if v1 and v2 are
648 /// column vectors, e.g. M = Mat3f::outerproduct(v1,v2);
649 template <typename T>
651 {
652  return Mat3<T>(v1[0]*v2[0], v1[0]*v2[1], v1[0]*v2[2],
653  v1[1]*v2[0], v1[1]*v2[1], v1[1]*v2[2],
654  v1[2]*v2[0], v1[2]*v2[1], v1[2]*v2[2]);
655 }// outerProduct
656 
657 
658 /// Interpolate the rotation between m1 and m2 using Mat::powSolve.
659 /// Unlike slerp, translation is not treated independently.
660 /// This results in smoother animation results.
661 template<typename T, typename T0>
662 Mat3<T> powLerp(const Mat3<T0> &m1, const Mat3<T0> &m2, T t)
663 {
664  Mat3<T> x = m1.inverse() * m2;
665  powSolve(x, x, t);
666  Mat3<T> m = m1 * x;
667  return m;
668 }
669 
670 
671 namespace mat3_internal {
672 
673 template<typename T>
674 inline void
675 pivot(int i, int j, Mat3<T>& S, Vec3<T>& D, Mat3<T>& Q)
676 {
677  const int& n = Mat3<T>::size; // should be 3
678  T temp;
679  /// scratch variables used in pivoting
680  double cotan_of_2_theta;
681  double tan_of_theta;
682  double cosin_of_theta;
683  double sin_of_theta;
684  double z;
685 
686  double Sij = S(i,j);
687 
688  double Sjj_minus_Sii = D[j] - D[i];
689 
690  if (fabs(Sjj_minus_Sii) * (10*math::Tolerance<T>::value()) > fabs(Sij)) {
691  tan_of_theta = Sij / Sjj_minus_Sii;
692  } else {
693  /// pivot on Sij
694  cotan_of_2_theta = 0.5*Sjj_minus_Sii / Sij ;
695 
696  if (cotan_of_2_theta < 0.) {
697  tan_of_theta =
698  -1./(sqrt(1. + cotan_of_2_theta*cotan_of_2_theta) - cotan_of_2_theta);
699  } else {
700  tan_of_theta =
701  1./(sqrt(1. + cotan_of_2_theta*cotan_of_2_theta) + cotan_of_2_theta);
702  }
703  }
704 
705  cosin_of_theta = 1./sqrt( 1. + tan_of_theta * tan_of_theta);
706  sin_of_theta = cosin_of_theta * tan_of_theta;
707  z = tan_of_theta * Sij;
708  S(i,j) = 0;
709  D[i] -= z;
710  D[j] += z;
711  for (int k = 0; k < i; ++k) {
712  temp = S(k,i);
713  S(k,i) = cosin_of_theta * temp - sin_of_theta * S(k,j);
714  S(k,j)= sin_of_theta * temp + cosin_of_theta * S(k,j);
715  }
716  for (int k = i+1; k < j; ++k) {
717  temp = S(i,k);
718  S(i,k) = cosin_of_theta * temp - sin_of_theta * S(k,j);
719  S(k,j) = sin_of_theta * temp + cosin_of_theta * S(k,j);
720  }
721  for (int k = j+1; k < n; ++k) {
722  temp = S(i,k);
723  S(i,k) = cosin_of_theta * temp - sin_of_theta * S(j,k);
724  S(j,k) = sin_of_theta * temp + cosin_of_theta * S(j,k);
725  }
726  for (int k = 0; k < n; ++k)
727  {
728  temp = Q(k,i);
729  Q(k,i) = cosin_of_theta * temp - sin_of_theta*Q(k,j);
730  Q(k,j) = sin_of_theta * temp + cosin_of_theta*Q(k,j);
731  }
732 }
733 
734 } // namespace mat3_internal
735 
736 
737 /// @brief Use Jacobi iterations to decompose a symmetric 3x3 matrix
738 /// (diagonalize and compute eigenvectors)
739 /// @details This is based on the "Efficient numerical diagonalization of Hermitian 3x3 matrices"
740 /// Joachim Kopp. arXiv.org preprint: physics/0610206
741 /// with the addition of largest pivot
742 template<typename T>
743 inline bool
745  unsigned int MAX_ITERATIONS=250)
746 {
747  /// use Givens rotation matrix to eliminate off-diagonal entries.
748  /// initialize the rotation matrix as idenity
749  Q = Mat3<T>::identity();
750  int n = Mat3<T>::size; // should be 3
751 
752  /// temp matrix. Assumed to be symmetric
753  Mat3<T> S(input);
754 
755  for (int i = 0; i < n; ++i) {
756  D[i] = S(i,i);
757  }
758 
759  unsigned int iterations(0);
760  /// Just iterate over all the non-diagonal enteries
761  /// using the largest as a pivot.
762  do {
763  /// check for absolute convergence
764  /// are symmetric off diagonals all zero
765  double er = 0;
766  for (int i = 0; i < n; ++i) {
767  for (int j = i+1; j < n; ++j) {
768  er += fabs(S(i,j));
769  }
770  }
771  if (std::abs(er) < math::Tolerance<T>::value()) {
772  return true;
773  }
774  iterations++;
775 
776  T max_element = 0;
777  int ip = 0;
778  int jp = 0;
779  /// loop over all the off-diagonals above the diagonal
780  for (int i = 0; i < n; ++i) {
781  for (int j = i+1; j < n; ++j){
782 
783  if ( fabs(D[i]) * (10*math::Tolerance<T>::value()) > fabs(S(i,j))) {
784  /// value too small to pivot on
785  S(i,j) = 0;
786  }
787  if (fabs(S(i,j)) > max_element) {
788  max_element = fabs(S(i,j));
789  ip = i;
790  jp = j;
791  }
792  }
793  }
794  mat3_internal::pivot(ip, jp, S, D, Q);
795  } while (iterations < MAX_ITERATIONS);
796 
797  return false;
798 }
799 
800 template<typename T>
801 inline Mat3<T>
802 Abs(const Mat3<T>& m)
803 {
804  Mat3<T> out;
805  const T* ip = m.asPointer();
806  T* op = out.asPointer();
807  for (unsigned i = 0; i < 9; ++i, ++op, ++ip) *op = math::Abs(*ip);
808  return out;
809 }
810 
811 template<typename Type1, typename Type2>
812 inline Mat3<Type1>
813 cwiseAdd(const Mat3<Type1>& m, const Type2 s)
814 {
815  Mat3<Type1> out;
816  const Type1* ip = m.asPointer();
817  Type1* op = out.asPointer();
818  for (unsigned i = 0; i < 9; ++i, ++op, ++ip) {
820  *op = *ip + s;
822  }
823  return out;
824 }
825 
826 template<typename T>
827 inline bool
828 cwiseLessThan(const Mat3<T>& m0, const Mat3<T>& m1)
829 {
830  return cwiseLessThan<3, T>(m0, m1);
831 }
832 
833 template<typename T>
834 inline bool
835 cwiseGreaterThan(const Mat3<T>& m0, const Mat3<T>& m1)
836 {
837  return cwiseGreaterThan<3, T>(m0, m1);
838 }
839 
842 using Mat3f = Mat3d;
843 
844 } // namespace math
845 
846 
847 template<> inline math::Mat3s zeroVal<math::Mat3s>() { return math::Mat3s::zero(); }
848 template<> inline math::Mat3d zeroVal<math::Mat3d>() { return math::Mat3d::zero(); }
849 
850 } // namespace OPENVDB_VERSION_NAME
851 } // namespace openvdb
852 
853 #endif // OPENVDB_MATH_MAT3_H_HAS_BEEN_INCLUDED
GLdouble s
Definition: glew.h:1390
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_END
Definition: Platform.h:188
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:568
GLenum GLenum GLenum input
Definition: glew.h:13879
void setToRotation(const Quat< T > &q)
Set this matrix to the rotation matrix specified by the quaternion.
Definition: Mat3.h:260
GLenum src
Definition: glew.h:2410
bool cwiseGreaterThan(const Mat< SIZE, T > &m0, const Mat< SIZE, T > &m1)
Definition: Mat.h:1044
void pivot(int i, int j, Mat3< T > &S, Vec3< T > &D, Mat3< T > &Q)
Definition: Mat3.h:675
Mat3< T > outerProduct(const Vec3< T > &v1, const Vec3< T > &v2)
Definition: Mat3.h:650
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:628
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:434
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:51
Mat3(const Mat4< T > &m)
Conversion from Mat4 (copies top left)
Definition: Mat3.h:118
Mat3 snapBasis(Axis axis, const Vec3< T > &direction)
Definition: Mat3.h:504
IMF_EXPORT IMATH_NAMESPACE::V3f direction(const IMATH_NAMESPACE::Box2i &dataWindow, const IMATH_NAMESPACE::V2f &pixelPosition)
void setToRotation(const Vec3< T > &axis, T angle)
Set this matrix to the rotation specified by axis and angle.
Definition: Mat3.h:265
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:584
Vec3< T > col(int j) const
Get jth column, e.g. Vec3d v = m.col(0);.
Definition: Mat3.h:175
Mat3< Type1 > cwiseAdd(const Mat3< Type1 > &m, const Type2 s)
Definition: Mat3.h:813
GLboolean GLboolean GLboolean GLboolean a
Definition: glew.h:9477
void setZero()
Set this matrix to zero.
Definition: Mat3.h:269
vfloat4 sqrt(const vfloat4 &a)
Definition: simd.h:7231
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:604
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:167
const GLdouble * m
Definition: glew.h:9124
const GLdouble * v
Definition: glew.h:1391
Mat3< T > operator-() const
Negation operator, for e.g. m1 = -m2;.
Definition: Mat3.h:322
GLdouble angle
Definition: glew.h:9135
Tolerance for floating-point comparison.
Definition: Math.h:137
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:202
Mat3()
Trivial constructor, the matrix is NOT initialized.
Definition: Mat3.h:36
GLdouble GLdouble z
Definition: glew.h:1559
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:216
Mat3(const Vec3< Source > &v1, const Vec3< Source > &v2, const Vec3< Source > &v3, bool rows=true)
Definition: Mat3.h:69
Vec3< T0 > transform(const Vec3< T0 > &v) const
Definition: Mat3.h:512
Mat3 transpose() const
returns transpose of this
Definition: Mat3.h:461
GLdouble GLdouble GLdouble GLdouble q
Definition: glew.h:1414
GLfloat GLfloat GLfloat v2
Definition: glew.h:1856
static const Mat3< T > & identity()
Predefined constant for identity matrix.
Definition: Mat3.h:128
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:244
bool isApproxEqual(const Type &a, const Type &b, const Type &tolerance)
Return true if a is equal to b to within the given tolerance.
Definition: Math.h:397
IMATH_INTERNAL_NAMESPACE_HEADER_ENTER T abs(T a)
Definition: ImathFun.h:55
GLclampf f
Definition: glew.h:3499
GLint GLint GLint GLint GLint x
Definition: glew.h:1252
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:744
Coord Abs(const Coord &xyz)
Definition: Coord.h:515
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:594
Mat3 inverse(T tolerance=0) const
Definition: Mat3.h:472
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:574
void setCol(int j, const Vec3< T > &v)
Set jth column to vector v.
Definition: Mat3.h:166
Vec3< T > row(int i) const
Get ith row, e.g. Vec3d v = m.row(1);.
Definition: Mat3.h:159
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:308
GLsizei n
Definition: glew.h:4040
const GLfloat * c
Definition: glew.h:16296
Mat3(const Mat3< Source > &m)
Conversion constructor.
Definition: Mat3.h:108
bool operator!=(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Inequality operator, does exact floating point comparisons.
Definition: Mat3.h:563
T angle(const Vec2< T > &v1, const Vec2< T > &v2)
Definition: Vec2.h:445
bool cwiseLessThan(const Mat< SIZE, T > &m0, const Mat< SIZE, T > &m1)
Definition: Mat.h:1030
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:354
T det() const
Determinant of matrix.
Definition: Mat3.h:486
GLfloat GLfloat GLfloat GLfloat h
Definition: glew.h:8011
Mat3 adjoint() const
Return the adjoint of this matrix, i.e., the transpose of its cofactor.
Definition: Mat3.h:445
const Vec2< S > & operator*=(Vec2< S > &v, const Matrix33< T > &m)
Definition: ImathMatrix.h:3322
static const Mat3< T > & zero()
Predefined constant for zero matrix.
Definition: Mat3.h:138
const Mat3< T > & operator*=(S scalar)
Multiplication operator, e.g. M = scalar * M;.
Definition: Mat3.h:338
Mat3(const Mat< 3, T > &m)
Copy constructor.
Definition: Mat3.h:97
void setSkew(const Vec3< T > &v)
Set the matrix as cross product of the given vector.
Definition: Mat3.h:254
GLdouble GLdouble GLdouble b
Definition: glew.h:9122
void setRow(int i, const Vec3< T > &v)
Set ith row to vector v.
Definition: Mat3.h:148
void setSymmetric(const Vec3< T > &vdiag, const Vec3< T > &vtri)
Set diagonal and symmetric triangular components.
Definition: Mat3.h:230
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:528
T * asPointer()
Direct access to the internal data.
Definition: Mat.h:116
#define OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN
Bracket code with OPENVDB_NO_TYPE_CONVERSION_WARNING_BEGIN/_END, to inhibit warnings about type conve...
Definition: Platform.h:187
T trace() const
Trace of matrix.
Definition: Mat3.h:495
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s1
Definition: glew.h:12681
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t1
Definition: glew.h:12681
OIIO_API bool copy(string_view from, string_view to, std::string &err)
const Mat3 & operator=(const Mat3< Source > &m)
Assignment operator.
Definition: Mat3.h:298
void setIdentity()
Set this matrix to identity.
Definition: Mat3.h:283
GLuint64EXT * result
Definition: glew.h:14007
Vec3< T0 > pretransform(const Vec3< T0 > &v) const
Definition: Mat3.h:520
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat s0
Definition: glew.h:12681
GLuint GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat GLfloat t0
Definition: glew.h:12681
Mat3< T > powLerp(const Mat3< T0 > &m1, const Mat3< T0 > &m2, T t)
Definition: Mat3.h:662
const Mat3< T > & operator*=(const Mat3< S > &m1)
Multiply this matrix by the given matrix.
Definition: Mat3.h:390
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:766
#define OPENVDB_VERSION_NAME
The version namespace name for this library version.
Definition: version.h:113
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:615
GLfloat GLfloat GLfloat GLfloat v3
Definition: glew.h:1860
GLdouble GLdouble t
Definition: glew.h:1398
MatType skew(const Vec3< typename MatType::value_type > &skew)
Return a matrix as the cross product of the given vector.
Definition: Mat.h:723
GLfloat GLfloat v1
Definition: glew.h:1852
bool operator==(const Mat3< T0 > &m0, const Mat3< T1 > &m1)
Equality operator, does exact floating point comparisons.
Definition: Mat3.h:549
T operator()(int i, int j) const
Definition: Mat3.h:194
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:82
Mat3 cofactor() const
Return the cofactor matrix of this matrix.
Definition: Mat3.h:430
void powSolve(const MatType &aA, MatType &aB, double aPower, double aTol=0.01)
Definition: Mat.h:837
GLboolean GLboolean g
Definition: glew.h:9477
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:372