OpenVDB  3.2.0
Mat4.h
Go to the documentation of this file.
1 //
3 // Copyright (c) 2012-2016 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_MAT4_H_HAS_BEEN_INCLUDED
32 #define OPENVDB_MATH_MAT4_H_HAS_BEEN_INCLUDED
33 
34 #include <openvdb/Exceptions.h>
35 #include <openvdb/Platform.h>
36 #include <iomanip>
37 #include <assert.h>
38 #include <math.h>
39 #include <algorithm>
40 #include "Math.h"
41 #include "Mat3.h"
42 #include "Vec3.h"
43 #include "Vec4.h"
44 
45 
46 namespace openvdb {
48 namespace OPENVDB_VERSION_NAME {
49 namespace math {
50 
51 template<typename T> class Vec4;
52 
53 
56 template<typename T>
57 class Mat4: public Mat<4, T>
58 {
59 public:
61  typedef T value_type;
62  typedef T ValueType;
63  typedef Mat<4, T> MyBase;
64 
66  Mat4() {}
67 
69 
75  template<typename Source>
76  Mat4(Source *a)
77  {
78  for (int i = 0; i < 16; i++) {
79  MyBase::mm[i] = a[i];
80  }
81  }
82 
84 
90  template<typename Source>
91  Mat4(Source a, Source b, Source c, Source d,
92  Source e, Source f, Source g, Source h,
93  Source i, Source j, Source k, Source l,
94  Source m, Source n, Source o, Source p)
95  {
96  MyBase::mm[ 0] = T(a);
97  MyBase::mm[ 1] = T(b);
98  MyBase::mm[ 2] = T(c);
99  MyBase::mm[ 3] = T(d);
100 
101  MyBase::mm[ 4] = T(e);
102  MyBase::mm[ 5] = T(f);
103  MyBase::mm[ 6] = T(g);
104  MyBase::mm[ 7] = T(h);
105 
106  MyBase::mm[ 8] = T(i);
107  MyBase::mm[ 9] = T(j);
108  MyBase::mm[10] = T(k);
109  MyBase::mm[11] = T(l);
110 
111  MyBase::mm[12] = T(m);
112  MyBase::mm[13] = T(n);
113  MyBase::mm[14] = T(o);
114  MyBase::mm[15] = T(p);
115  }
116 
119  template<typename Source>
120  Mat4(const Vec4<Source> &v1, const Vec4<Source> &v2,
121  const Vec4<Source> &v3, const Vec4<Source> &v4, bool rows = true)
122  {
123  if (rows) {
124  this->setRows(v1, v2, v3, v4);
125  } else {
126  this->setColumns(v1, v2, v3, v4);
127  }
128  }
129 
131  Mat4(const Mat<4, T> &m)
132  {
133  for (int i = 0; i < 4; ++i) {
134  for (int j = 0; j < 4; ++j) {
135  MyBase::mm[i*4 + j] = m[i][j];
136  }
137  }
138  }
139 
141  template<typename Source>
142  explicit Mat4(const Mat4<Source> &m)
143  {
144  const Source *src = m.asPointer();
145 
146  for (int i=0; i<16; ++i) {
147  MyBase::mm[i] = static_cast<T>(src[i]);
148  }
149  }
150 
152  static const Mat4<T>& identity() {
153  return sIdentity;
154  }
155 
157  static const Mat4<T>& zero() {
158  return sZero;
159  }
160 
162  void setRow(int i, const Vec4<T> &v)
163  {
164  // assert(i>=0 && i<4);
165  int i4 = i * 4;
166  MyBase::mm[i4+0] = v[0];
167  MyBase::mm[i4+1] = v[1];
168  MyBase::mm[i4+2] = v[2];
169  MyBase::mm[i4+3] = v[3];
170  }
171 
173  Vec4<T> row(int i) const
174  {
175  // assert(i>=0 && i<3);
176  return Vec4<T>((*this)(i,0), (*this)(i,1), (*this)(i,2), (*this)(i,3));
177  }
178 
180  void setCol(int j, const Vec4<T>& v)
181  {
182  // assert(j>=0 && j<4);
183  MyBase::mm[ 0+j] = v[0];
184  MyBase::mm[ 4+j] = v[1];
185  MyBase::mm[ 8+j] = v[2];
186  MyBase::mm[12+j] = v[3];
187  }
188 
190  Vec4<T> col(int j) const
191  {
192  // assert(j>=0 && j<4);
193  return Vec4<T>((*this)(0,j), (*this)(1,j), (*this)(2,j), (*this)(3,j));
194  }
195 
197  T* operator[](int i) { return &(MyBase::mm[i<<2]); }
200  const T* operator[](int i) const { return &(MyBase::mm[i<<2]); }
202 
204  T* asPointer() {return MyBase::mm;}
205  const T* asPointer() const {return MyBase::mm;}
206 
210  T& operator()(int i, int j)
211  {
212  // assert(i>=0 && i<4);
213  // assert(j>=0 && j<4);
214  return MyBase::mm[4*i+j];
215  }
216 
220  T operator()(int i, int j) const
221  {
222  // assert(i>=0 && i<4);
223  // assert(j>=0 && j<4);
224  return MyBase::mm[4*i+j];
225  }
226 
228  void setRows(const Vec4<T> &v1, const Vec4<T> &v2,
229  const Vec4<T> &v3, const Vec4<T> &v4)
230  {
231  MyBase::mm[ 0] = v1[0];
232  MyBase::mm[ 1] = v1[1];
233  MyBase::mm[ 2] = v1[2];
234  MyBase::mm[ 3] = v1[3];
235 
236  MyBase::mm[ 4] = v2[0];
237  MyBase::mm[ 5] = v2[1];
238  MyBase::mm[ 6] = v2[2];
239  MyBase::mm[ 7] = v2[3];
240 
241  MyBase::mm[ 8] = v3[0];
242  MyBase::mm[ 9] = v3[1];
243  MyBase::mm[10] = v3[2];
244  MyBase::mm[11] = v3[3];
245 
246  MyBase::mm[12] = v4[0];
247  MyBase::mm[13] = v4[1];
248  MyBase::mm[14] = v4[2];
249  MyBase::mm[15] = v4[3];
250  }
251 
253  void setColumns(const Vec4<T> &v1, const Vec4<T> &v2,
254  const Vec4<T> &v3, const Vec4<T> &v4)
255  {
256  MyBase::mm[ 0] = v1[0];
257  MyBase::mm[ 1] = v2[0];
258  MyBase::mm[ 2] = v3[0];
259  MyBase::mm[ 3] = v4[0];
260 
261  MyBase::mm[ 4] = v1[1];
262  MyBase::mm[ 5] = v2[1];
263  MyBase::mm[ 6] = v3[1];
264  MyBase::mm[ 7] = v4[1];
265 
266  MyBase::mm[ 8] = v1[2];
267  MyBase::mm[ 9] = v2[2];
268  MyBase::mm[10] = v3[2];
269  MyBase::mm[11] = v4[2];
270 
271  MyBase::mm[12] = v1[3];
272  MyBase::mm[13] = v2[3];
273  MyBase::mm[14] = v3[3];
274  MyBase::mm[15] = v4[3];
275  }
276 
278  OPENVDB_DEPRECATED void setBasis(const Vec4<T> &v1, const Vec4<T> &v2,
279  const Vec4<T> &v3, const Vec4<T> &v4)
280  {
281  this->setRows(v1, v2, v3, v4);
282  }
283 
284 
285  // Set "this" matrix to zero
286  void setZero()
287  {
288  MyBase::mm[ 0] = 0;
289  MyBase::mm[ 1] = 0;
290  MyBase::mm[ 2] = 0;
291  MyBase::mm[ 3] = 0;
292  MyBase::mm[ 4] = 0;
293  MyBase::mm[ 5] = 0;
294  MyBase::mm[ 6] = 0;
295  MyBase::mm[ 7] = 0;
296  MyBase::mm[ 8] = 0;
297  MyBase::mm[ 9] = 0;
298  MyBase::mm[10] = 0;
299  MyBase::mm[11] = 0;
300  MyBase::mm[12] = 0;
301  MyBase::mm[13] = 0;
302  MyBase::mm[14] = 0;
303  MyBase::mm[15] = 0;
304  }
305 
307  void setIdentity()
308  {
309  MyBase::mm[ 0] = 1;
310  MyBase::mm[ 1] = 0;
311  MyBase::mm[ 2] = 0;
312  MyBase::mm[ 3] = 0;
313 
314  MyBase::mm[ 4] = 0;
315  MyBase::mm[ 5] = 1;
316  MyBase::mm[ 6] = 0;
317  MyBase::mm[ 7] = 0;
318 
319  MyBase::mm[ 8] = 0;
320  MyBase::mm[ 9] = 0;
321  MyBase::mm[10] = 1;
322  MyBase::mm[11] = 0;
323 
324  MyBase::mm[12] = 0;
325  MyBase::mm[13] = 0;
326  MyBase::mm[14] = 0;
327  MyBase::mm[15] = 1;
328  }
329 
330 
332  void setMat3(const Mat3<T> &m)
333  {
334  for (int i = 0; i < 3; i++)
335  for (int j=0; j < 3; j++)
336  MyBase::mm[i*4+j] = m[i][j];
337  }
338 
339  Mat3<T> getMat3() const
340  {
341  Mat3<T> m;
342 
343  for (int i = 0; i < 3; i++)
344  for (int j = 0; j < 3; j++)
345  m[i][j] = MyBase::mm[i*4+j];
346 
347  return m;
348  }
349 
352  {
353  return Vec3<T>(MyBase::mm[12], MyBase::mm[13], MyBase::mm[14]);
354  }
355 
356  void setTranslation(const Vec3<T> &t)
357  {
358  MyBase::mm[12] = t[0];
359  MyBase::mm[13] = t[1];
360  MyBase::mm[14] = t[2];
361  }
362 
364  template<typename Source>
365  const Mat4& operator=(const Mat4<Source> &m)
366  {
367  const Source *src = m.asPointer();
368 
369  // don't suppress warnings when assigning from different numerical types
370  std::copy(src, (src + this->numElements()), MyBase::mm);
371  return *this;
372  }
373 
375  bool eq(const Mat4 &m, T eps=1.0e-8) const
376  {
377  for (int i = 0; i < 16; i++) {
378  if (!isApproxEqual(MyBase::mm[i], m.mm[i], eps))
379  return false;
380  }
381  return true;
382  }
383 
386  {
387  return Mat4<T>(
388  -MyBase::mm[ 0], -MyBase::mm[ 1], -MyBase::mm[ 2], -MyBase::mm[ 3],
389  -MyBase::mm[ 4], -MyBase::mm[ 5], -MyBase::mm[ 6], -MyBase::mm[ 7],
390  -MyBase::mm[ 8], -MyBase::mm[ 9], -MyBase::mm[10], -MyBase::mm[11],
391  -MyBase::mm[12], -MyBase::mm[13], -MyBase::mm[14], -MyBase::mm[15]
392  );
393  } // trivial
394 
396  template <typename S>
397  const Mat4<T>& operator*=(S scalar)
398  {
399  MyBase::mm[ 0] *= scalar;
400  MyBase::mm[ 1] *= scalar;
401  MyBase::mm[ 2] *= scalar;
402  MyBase::mm[ 3] *= scalar;
403 
404  MyBase::mm[ 4] *= scalar;
405  MyBase::mm[ 5] *= scalar;
406  MyBase::mm[ 6] *= scalar;
407  MyBase::mm[ 7] *= scalar;
408 
409  MyBase::mm[ 8] *= scalar;
410  MyBase::mm[ 9] *= scalar;
411  MyBase::mm[10] *= scalar;
412  MyBase::mm[11] *= scalar;
413 
414  MyBase::mm[12] *= scalar;
415  MyBase::mm[13] *= scalar;
416  MyBase::mm[14] *= scalar;
417  MyBase::mm[15] *= scalar;
418  return *this;
419  }
420 
422  template <typename S>
423  const Mat4<T> &operator+=(const Mat4<S> &m1)
424  {
425  const S* s = m1.asPointer();
426 
427  MyBase::mm[ 0] += s[ 0];
428  MyBase::mm[ 1] += s[ 1];
429  MyBase::mm[ 2] += s[ 2];
430  MyBase::mm[ 3] += s[ 3];
431 
432  MyBase::mm[ 4] += s[ 4];
433  MyBase::mm[ 5] += s[ 5];
434  MyBase::mm[ 6] += s[ 6];
435  MyBase::mm[ 7] += s[ 7];
436 
437  MyBase::mm[ 8] += s[ 8];
438  MyBase::mm[ 9] += s[ 9];
439  MyBase::mm[10] += s[10];
440  MyBase::mm[11] += s[11];
441 
442  MyBase::mm[12] += s[12];
443  MyBase::mm[13] += s[13];
444  MyBase::mm[14] += s[14];
445  MyBase::mm[15] += s[15];
446 
447  return *this;
448  }
449 
451  template <typename S>
452  const Mat4<T> &operator-=(const Mat4<S> &m1)
453  {
454  const S* s = m1.asPointer();
455 
456  MyBase::mm[ 0] -= s[ 0];
457  MyBase::mm[ 1] -= s[ 1];
458  MyBase::mm[ 2] -= s[ 2];
459  MyBase::mm[ 3] -= s[ 3];
460 
461  MyBase::mm[ 4] -= s[ 4];
462  MyBase::mm[ 5] -= s[ 5];
463  MyBase::mm[ 6] -= s[ 6];
464  MyBase::mm[ 7] -= s[ 7];
465 
466  MyBase::mm[ 8] -= s[ 8];
467  MyBase::mm[ 9] -= s[ 9];
468  MyBase::mm[10] -= s[10];
469  MyBase::mm[11] -= s[11];
470 
471  MyBase::mm[12] -= s[12];
472  MyBase::mm[13] -= s[13];
473  MyBase::mm[14] -= s[14];
474  MyBase::mm[15] -= s[15];
475 
476  return *this;
477  }
478 
480  template <typename S>
481  const Mat4<T> &operator*=(const Mat4<S> &m1)
482  {
483  Mat4<T> m0(*this);
484 
485  const T* s0 = m0.asPointer();
486  const S* s1 = m1.asPointer();
487 
488  for (int i = 0; i < 4; i++) {
489  int i4 = 4 * i;
490  MyBase::mm[i4+0] = static_cast<T>(s0[i4+0] * s1[ 0] +
491  s0[i4+1] * s1[ 4] +
492  s0[i4+2] * s1[ 8] +
493  s0[i4+3] * s1[12]);
494 
495  MyBase::mm[i4+1] = static_cast<T>(s0[i4+0] * s1[ 1] +
496  s0[i4+1] * s1[ 5] +
497  s0[i4+2] * s1[ 9] +
498  s0[i4+3] * s1[13]);
499 
500  MyBase::mm[i4+2] = static_cast<T>(s0[i4+0] * s1[ 2] +
501  s0[i4+1] * s1[ 6] +
502  s0[i4+2] * s1[10] +
503  s0[i4+3] * s1[14]);
504 
505  MyBase::mm[i4+3] = static_cast<T>(s0[i4+0] * s1[ 3] +
506  s0[i4+1] * s1[ 7] +
507  s0[i4+2] * s1[11] +
508  s0[i4+3] * s1[15]);
509  }
510  return *this;
511  }
512 
514  Mat4 transpose() const
515  {
516  return Mat4<T>(
517  MyBase::mm[ 0], MyBase::mm[ 4], MyBase::mm[ 8], MyBase::mm[12],
518  MyBase::mm[ 1], MyBase::mm[ 5], MyBase::mm[ 9], MyBase::mm[13],
519  MyBase::mm[ 2], MyBase::mm[ 6], MyBase::mm[10], MyBase::mm[14],
520  MyBase::mm[ 3], MyBase::mm[ 7], MyBase::mm[11], MyBase::mm[15]
521  );
522  }
523 
524 
527  Mat4 inverse(T tolerance = 0) const
528  {
529  //
530  // inv [ A | b ] = [ E | f ] A: 3x3, b: 3x1, c': 1x3 d: 1x1
531  // [ c' | d ] [ g' | h ]
532  //
533  // If A is invertible use
534  //
535  // E = A^-1 + p*h*r
536  // p = A^-1 * b
537  // f = -p * h
538  // g' = -h * c'
539  // h = 1 / (d - c'*p)
540  // r' = c'*A^-1
541  //
542  // Otherwise use gauss-jordan elimination
543  //
544 
545  //
546  // We create this alias to ourself so we can easily use own subscript
547  // operator.
548  const Mat4<T>& m(*this);
549 
550  T m0011 = m[0][0] * m[1][1];
551  T m0012 = m[0][0] * m[1][2];
552  T m0110 = m[0][1] * m[1][0];
553  T m0210 = m[0][2] * m[1][0];
554  T m0120 = m[0][1] * m[2][0];
555  T m0220 = m[0][2] * m[2][0];
556 
557  T detA = m0011 * m[2][2] - m0012 * m[2][1] - m0110 * m[2][2]
558  + m0210 * m[2][1] + m0120 * m[1][2] - m0220 * m[1][1];
559 
560  bool hasPerspective =
561  (!isExactlyEqual(m[0][3], T(0.0)) ||
562  !isExactlyEqual(m[1][3], T(0.0)) ||
563  !isExactlyEqual(m[2][3], T(0.0)) ||
564  !isExactlyEqual(m[3][3], T(1.0)));
565 
566  T det;
567  if (hasPerspective) {
568  det = m[0][3] * det3(m, 1,2,3, 0,2,1)
569  + m[1][3] * det3(m, 2,0,3, 0,2,1)
570  + m[2][3] * det3(m, 3,0,1, 0,2,1)
571  + m[3][3] * detA;
572  } else {
573  det = detA * m[3][3];
574  }
575 
576  Mat4<T> inv;
577  bool invertible;
578 
579  if (isApproxEqual(det,T(0.0),tolerance)) {
580  invertible = false;
581 
582  } else if (isApproxEqual(detA,T(0.0),T(1e-8))) {
583  // det is too small to rely on inversion by subblocks
584  invertible = m.invert(inv, tolerance);
585 
586  } else {
587  invertible = true;
588  detA = 1.0 / detA;
589 
590  //
591  // Calculate A^-1
592  //
593  inv[0][0] = detA * ( m[1][1] * m[2][2] - m[1][2] * m[2][1]);
594  inv[0][1] = detA * (-m[0][1] * m[2][2] + m[0][2] * m[2][1]);
595  inv[0][2] = detA * ( m[0][1] * m[1][2] - m[0][2] * m[1][1]);
596 
597  inv[1][0] = detA * (-m[1][0] * m[2][2] + m[1][2] * m[2][0]);
598  inv[1][1] = detA * ( m[0][0] * m[2][2] - m0220);
599  inv[1][2] = detA * ( m0210 - m0012);
600 
601  inv[2][0] = detA * ( m[1][0] * m[2][1] - m[1][1] * m[2][0]);
602  inv[2][1] = detA * ( m0120 - m[0][0] * m[2][1]);
603  inv[2][2] = detA * ( m0011 - m0110);
604 
605  if (hasPerspective) {
606  //
607  // Calculate r, p, and h
608  //
609  Vec3<T> r;
610  r[0] = m[3][0] * inv[0][0] + m[3][1] * inv[1][0]
611  + m[3][2] * inv[2][0];
612  r[1] = m[3][0] * inv[0][1] + m[3][1] * inv[1][1]
613  + m[3][2] * inv[2][1];
614  r[2] = m[3][0] * inv[0][2] + m[3][1] * inv[1][2]
615  + m[3][2] * inv[2][2];
616 
617  Vec3<T> p;
618  p[0] = inv[0][0] * m[0][3] + inv[0][1] * m[1][3]
619  + inv[0][2] * m[2][3];
620  p[1] = inv[1][0] * m[0][3] + inv[1][1] * m[1][3]
621  + inv[1][2] * m[2][3];
622  p[2] = inv[2][0] * m[0][3] + inv[2][1] * m[1][3]
623  + inv[2][2] * m[2][3];
624 
625  T h = m[3][3] - p.dot(Vec3<T>(m[3][0],m[3][1],m[3][2]));
626  if (isApproxEqual(h,T(0.0),tolerance)) {
627  invertible = false;
628 
629  } else {
630  h = 1.0 / h;
631 
632  //
633  // Calculate h, g, and f
634  //
635  inv[3][3] = h;
636  inv[3][0] = -h * r[0];
637  inv[3][1] = -h * r[1];
638  inv[3][2] = -h * r[2];
639 
640  inv[0][3] = -h * p[0];
641  inv[1][3] = -h * p[1];
642  inv[2][3] = -h * p[2];
643 
644  //
645  // Calculate E
646  //
647  p *= h;
648  inv[0][0] += p[0] * r[0];
649  inv[0][1] += p[0] * r[1];
650  inv[0][2] += p[0] * r[2];
651  inv[1][0] += p[1] * r[0];
652  inv[1][1] += p[1] * r[1];
653  inv[1][2] += p[1] * r[2];
654  inv[2][0] += p[2] * r[0];
655  inv[2][1] += p[2] * r[1];
656  inv[2][2] += p[2] * r[2];
657  }
658  } else {
659  // Equations are much simpler in the non-perspective case
660  inv[3][0] = - (m[3][0] * inv[0][0] + m[3][1] * inv[1][0]
661  + m[3][2] * inv[2][0]);
662  inv[3][1] = - (m[3][0] * inv[0][1] + m[3][1] * inv[1][1]
663  + m[3][2] * inv[2][1]);
664  inv[3][2] = - (m[3][0] * inv[0][2] + m[3][1] * inv[1][2]
665  + m[3][2] * inv[2][2]);
666  inv[0][3] = 0.0;
667  inv[1][3] = 0.0;
668  inv[2][3] = 0.0;
669  inv[3][3] = 1.0;
670  }
671  }
672 
673  if (!invertible) OPENVDB_THROW(ArithmeticError, "Inversion of singular 4x4 matrix");
674  return inv;
675  }
676 
677 
679  T det() const
680  {
681  const T *ap;
682  Mat3<T> submat;
683  T det;
684  T *sp;
685  int i, j, k, sign;
686 
687  det = 0;
688  sign = 1;
689  for (i = 0; i < 4; i++) {
690  ap = &MyBase::mm[ 0];
691  sp = submat.asPointer();
692  for (j = 0; j < 4; j++) {
693  for (k = 0; k < 4; k++) {
694  if ((k != i) && (j != 0)) {
695  *sp++ = *ap;
696  }
697  ap++;
698  }
699  }
700 
701  det += sign * MyBase::mm[i] * submat.det();
702  sign = -sign;
703  }
704 
705  return det;
706  }
707 
709  static Mat4 translation(const Vec3d& v)
710  {
711  return Mat4(
712  T(1), T(0), T(0), T(0),
713  T(0), T(1), T(0), T(0),
714  T(0), T(0), T(1), T(0),
715  T(v.x()), T(v.y()),T(v.z()), T(1));
716  }
717 
719  template <typename T0>
720  void setToTranslation(const Vec3<T0>& v)
721  {
722  MyBase::mm[ 0] = 1;
723  MyBase::mm[ 1] = 0;
724  MyBase::mm[ 2] = 0;
725  MyBase::mm[ 3] = 0;
726 
727  MyBase::mm[ 4] = 0;
728  MyBase::mm[ 5] = 1;
729  MyBase::mm[ 6] = 0;
730  MyBase::mm[ 7] = 0;
731 
732  MyBase::mm[ 8] = 0;
733  MyBase::mm[ 9] = 0;
734  MyBase::mm[10] = 1;
735  MyBase::mm[11] = 0;
736 
737  MyBase::mm[12] = v.x();
738  MyBase::mm[13] = v.y();
739  MyBase::mm[14] = v.z();
740  MyBase::mm[15] = 1;
741  }
742 
744  template <typename T0>
745  void preTranslate(const Vec3<T0>& tr)
746  {
747  Vec3<T> tmp(tr.x(), tr.y(), tr.z());
748  Mat4<T> Tr = Mat4<T>::translation(tmp);
749 
750  *this = Tr * (*this);
751 
752  }
753 
755  template <typename T0>
756  void postTranslate(const Vec3<T0>& tr)
757  {
758  Vec3<T> tmp(tr.x(), tr.y(), tr.z());
759  Mat4<T> Tr = Mat4<T>::translation(tmp);
760 
761  *this = (*this) * Tr;
762 
763  }
764 
765 
767  template <typename T0>
768  void setToScale(const Vec3<T0>& v)
769  {
770  this->setIdentity();
771  MyBase::mm[ 0] = v.x();
772  MyBase::mm[ 5] = v.y();
773  MyBase::mm[10] = v.z();
774  }
775 
776  // Left multiples by the specified scale matrix, i.e. Sc * (*this)
777  template <typename T0>
778  void preScale(const Vec3<T0>& v)
779  {
780  MyBase::mm[ 0] *= v.x();
781  MyBase::mm[ 1] *= v.x();
782  MyBase::mm[ 2] *= v.x();
783  MyBase::mm[ 3] *= v.x();
784 
785  MyBase::mm[ 4] *= v.y();
786  MyBase::mm[ 5] *= v.y();
787  MyBase::mm[ 6] *= v.y();
788  MyBase::mm[ 7] *= v.y();
789 
790  MyBase::mm[ 8] *= v.z();
791  MyBase::mm[ 9] *= v.z();
792  MyBase::mm[10] *= v.z();
793  MyBase::mm[11] *= v.z();
794  }
795 
796 
797 
798  // Right multiples by the specified scale matrix, i.e. (*this) * Sc
799  template <typename T0>
800  void postScale(const Vec3<T0>& v)
801  {
802 
803  MyBase::mm[ 0] *= v.x();
804  MyBase::mm[ 1] *= v.y();
805  MyBase::mm[ 2] *= v.z();
806 
807  MyBase::mm[ 4] *= v.x();
808  MyBase::mm[ 5] *= v.y();
809  MyBase::mm[ 6] *= v.z();
810 
811  MyBase::mm[ 8] *= v.x();
812  MyBase::mm[ 9] *= v.y();
813  MyBase::mm[10] *= v.z();
814 
815  MyBase::mm[12] *= v.x();
816  MyBase::mm[13] *= v.y();
817  MyBase::mm[14] *= v.z();
818 
819  }
820 
821 
825  void setToRotation(Axis axis, T angle) {*this = rotation<Mat4<T> >(axis, angle);}
826 
830  void setToRotation(const Vec3<T>& axis, T angle) {*this = rotation<Mat4<T> >(axis, angle);}
831 
834  void setToRotation(const Vec3<T>& v1, const Vec3<T>& v2) {*this = rotation<Mat4<T> >(v1, v2);}
835 
836 
840  void preRotate(Axis axis, T angle)
841  {
842  T c = static_cast<T>(cos(angle));
843  T s = -static_cast<T>(sin(angle)); // the "-" makes it clockwise
844 
845  switch (axis) {
846  case X_AXIS:
847  {
848  T a4, a5, a6, a7;
849 
850  a4 = c * MyBase::mm[ 4] - s * MyBase::mm[ 8];
851  a5 = c * MyBase::mm[ 5] - s * MyBase::mm[ 9];
852  a6 = c * MyBase::mm[ 6] - s * MyBase::mm[10];
853  a7 = c * MyBase::mm[ 7] - s * MyBase::mm[11];
854 
855 
856  MyBase::mm[ 8] = s * MyBase::mm[ 4] + c * MyBase::mm[ 8];
857  MyBase::mm[ 9] = s * MyBase::mm[ 5] + c * MyBase::mm[ 9];
858  MyBase::mm[10] = s * MyBase::mm[ 6] + c * MyBase::mm[10];
859  MyBase::mm[11] = s * MyBase::mm[ 7] + c * MyBase::mm[11];
860 
861  MyBase::mm[ 4] = a4;
862  MyBase::mm[ 5] = a5;
863  MyBase::mm[ 6] = a6;
864  MyBase::mm[ 7] = a7;
865  }
866  break;
867 
868  case Y_AXIS:
869  {
870  T a0, a1, a2, a3;
871 
872  a0 = c * MyBase::mm[ 0] + s * MyBase::mm[ 8];
873  a1 = c * MyBase::mm[ 1] + s * MyBase::mm[ 9];
874  a2 = c * MyBase::mm[ 2] + s * MyBase::mm[10];
875  a3 = c * MyBase::mm[ 3] + s * MyBase::mm[11];
876 
877  MyBase::mm[ 8] = -s * MyBase::mm[ 0] + c * MyBase::mm[ 8];
878  MyBase::mm[ 9] = -s * MyBase::mm[ 1] + c * MyBase::mm[ 9];
879  MyBase::mm[10] = -s * MyBase::mm[ 2] + c * MyBase::mm[10];
880  MyBase::mm[11] = -s * MyBase::mm[ 3] + c * MyBase::mm[11];
881 
882 
883  MyBase::mm[ 0] = a0;
884  MyBase::mm[ 1] = a1;
885  MyBase::mm[ 2] = a2;
886  MyBase::mm[ 3] = a3;
887  }
888  break;
889 
890  case Z_AXIS:
891  {
892  T a0, a1, a2, a3;
893 
894  a0 = c * MyBase::mm[ 0] - s * MyBase::mm[ 4];
895  a1 = c * MyBase::mm[ 1] - s * MyBase::mm[ 5];
896  a2 = c * MyBase::mm[ 2] - s * MyBase::mm[ 6];
897  a3 = c * MyBase::mm[ 3] - s * MyBase::mm[ 7];
898 
899  MyBase::mm[ 4] = s * MyBase::mm[ 0] + c * MyBase::mm[ 4];
900  MyBase::mm[ 5] = s * MyBase::mm[ 1] + c * MyBase::mm[ 5];
901  MyBase::mm[ 6] = s * MyBase::mm[ 2] + c * MyBase::mm[ 6];
902  MyBase::mm[ 7] = s * MyBase::mm[ 3] + c * MyBase::mm[ 7];
903 
904  MyBase::mm[ 0] = a0;
905  MyBase::mm[ 1] = a1;
906  MyBase::mm[ 2] = a2;
907  MyBase::mm[ 3] = a3;
908  }
909  break;
910 
911  default:
912  assert(axis==X_AXIS || axis==Y_AXIS || axis==Z_AXIS);
913  }
914  }
915 
916 
920  void postRotate(Axis axis, T angle)
921  {
922  T c = static_cast<T>(cos(angle));
923  T s = -static_cast<T>(sin(angle)); // the "-" makes it clockwise
924 
925 
926 
927  switch (axis) {
928  case X_AXIS:
929  {
930  T a2, a6, a10, a14;
931 
932  a2 = c * MyBase::mm[ 2] - s * MyBase::mm[ 1];
933  a6 = c * MyBase::mm[ 6] - s * MyBase::mm[ 5];
934  a10 = c * MyBase::mm[10] - s * MyBase::mm[ 9];
935  a14 = c * MyBase::mm[14] - s * MyBase::mm[13];
936 
937 
938  MyBase::mm[ 1] = c * MyBase::mm[ 1] + s * MyBase::mm[ 2];
939  MyBase::mm[ 5] = c * MyBase::mm[ 5] + s * MyBase::mm[ 6];
940  MyBase::mm[ 9] = c * MyBase::mm[ 9] + s * MyBase::mm[10];
941  MyBase::mm[13] = c * MyBase::mm[13] + s * MyBase::mm[14];
942 
943  MyBase::mm[ 2] = a2;
944  MyBase::mm[ 6] = a6;
945  MyBase::mm[10] = a10;
946  MyBase::mm[14] = a14;
947  }
948  break;
949 
950  case Y_AXIS:
951  {
952  T a2, a6, a10, a14;
953 
954  a2 = c * MyBase::mm[ 2] + s * MyBase::mm[ 0];
955  a6 = c * MyBase::mm[ 6] + s * MyBase::mm[ 4];
956  a10 = c * MyBase::mm[10] + s * MyBase::mm[ 8];
957  a14 = c * MyBase::mm[14] + s * MyBase::mm[12];
958 
959  MyBase::mm[ 0] = c * MyBase::mm[ 0] - s * MyBase::mm[ 2];
960  MyBase::mm[ 4] = c * MyBase::mm[ 4] - s * MyBase::mm[ 6];
961  MyBase::mm[ 8] = c * MyBase::mm[ 8] - s * MyBase::mm[10];
962  MyBase::mm[12] = c * MyBase::mm[12] - s * MyBase::mm[14];
963 
964  MyBase::mm[ 2] = a2;
965  MyBase::mm[ 6] = a6;
966  MyBase::mm[10] = a10;
967  MyBase::mm[14] = a14;
968  }
969  break;
970 
971  case Z_AXIS:
972  {
973  T a1, a5, a9, a13;
974 
975  a1 = c * MyBase::mm[ 1] - s * MyBase::mm[ 0];
976  a5 = c * MyBase::mm[ 5] - s * MyBase::mm[ 4];
977  a9 = c * MyBase::mm[ 9] - s * MyBase::mm[ 8];
978  a13 = c * MyBase::mm[13] - s * MyBase::mm[12];
979 
980  MyBase::mm[ 0] = c * MyBase::mm[ 0] + s * MyBase::mm[ 1];
981  MyBase::mm[ 4] = c * MyBase::mm[ 4] + s * MyBase::mm[ 5];
982  MyBase::mm[ 8] = c * MyBase::mm[ 8] + s * MyBase::mm[ 9];
983  MyBase::mm[12] = c * MyBase::mm[12] + s * MyBase::mm[13];
984 
985  MyBase::mm[ 1] = a1;
986  MyBase::mm[ 5] = a5;
987  MyBase::mm[ 9] = a9;
988  MyBase::mm[13] = a13;
989 
990  }
991  break;
992 
993  default:
994  assert(axis==X_AXIS || axis==Y_AXIS || axis==Z_AXIS);
995  }
996  }
997 
1002  void setToShear(Axis axis0, Axis axis1, T shearby)
1003  {
1004  *this = shear<Mat4<T> >(axis0, axis1, shearby);
1005  }
1006 
1007 
1010  void preShear(Axis axis0, Axis axis1, T shear)
1011  {
1012  int index0 = static_cast<int>(axis0);
1013  int index1 = static_cast<int>(axis1);
1014 
1015  // to row "index1" add a multiple of the index0 row
1016  MyBase::mm[index1 * 4 + 0] += shear * MyBase::mm[index0 * 4 + 0];
1017  MyBase::mm[index1 * 4 + 1] += shear * MyBase::mm[index0 * 4 + 1];
1018  MyBase::mm[index1 * 4 + 2] += shear * MyBase::mm[index0 * 4 + 2];
1019  MyBase::mm[index1 * 4 + 3] += shear * MyBase::mm[index0 * 4 + 3];
1020  }
1021 
1022 
1025  void postShear(Axis axis0, Axis axis1, T shear)
1026  {
1027  int index0 = static_cast<int>(axis0);
1028  int index1 = static_cast<int>(axis1);
1029 
1030  // to collumn "index0" add a multiple of the index1 row
1031  MyBase::mm[index0 + 0] += shear * MyBase::mm[index1 + 0];
1032  MyBase::mm[index0 + 4] += shear * MyBase::mm[index1 + 4];
1033  MyBase::mm[index0 + 8] += shear * MyBase::mm[index1 + 8];
1034  MyBase::mm[index0 + 12] += shear * MyBase::mm[index1 + 12];
1035 
1036  }
1037 
1039  template<typename T0>
1040  Vec4<T0> transform(const Vec4<T0> &v) const
1041  {
1042  return static_cast< Vec4<T0> >(v * *this);
1043  }
1044 
1046  template<typename T0>
1047  Vec3<T0> transform(const Vec3<T0> &v) const
1048  {
1049  return static_cast< Vec3<T0> >(v * *this);
1050  }
1051 
1053  template<typename T0>
1055  {
1056  return static_cast< Vec4<T0> >(*this * v);
1057  }
1058 
1060  template<typename T0>
1062  {
1063  return static_cast< Vec3<T0> >(*this * v);
1064  }
1065 
1067  template<typename T0>
1068  Vec3<T0> transformH(const Vec3<T0> &p) const
1069  {
1070  T0 w;
1071 
1072  // w = p * (*this).col(3);
1073  w = static_cast<T0>(p[0] * MyBase::mm[ 3] + p[1] * MyBase::mm[ 7]
1074  + p[2] * MyBase::mm[11] + MyBase::mm[15]);
1075 
1076  if ( !isExactlyEqual(w , 0.0) ) {
1077  return Vec3<T0>(static_cast<T0>((p[0] * MyBase::mm[ 0] + p[1] * MyBase::mm[ 4] +
1078  p[2] * MyBase::mm[ 8] + MyBase::mm[12]) / w),
1079  static_cast<T0>((p[0] * MyBase::mm[ 1] + p[1] * MyBase::mm[ 5] +
1080  p[2] * MyBase::mm[ 9] + MyBase::mm[13]) / w),
1081  static_cast<T0>((p[0] * MyBase::mm[ 2] + p[1] * MyBase::mm[ 6] +
1082  p[2] * MyBase::mm[10] + MyBase::mm[14]) / w));
1083  }
1084 
1085  return Vec3<T0>(0, 0, 0);
1086  }
1087 
1089  template<typename T0>
1091  {
1092  T0 w;
1093 
1094  // w = p * (*this).col(3);
1095  w = p[0] * MyBase::mm[12] + p[1] * MyBase::mm[13] + p[2] * MyBase::mm[14] + MyBase::mm[15];
1096 
1097  if ( !isExactlyEqual(w , 0.0) ) {
1098  return Vec3<T0>(static_cast<T0>((p[0] * MyBase::mm[ 0] + p[1] * MyBase::mm[ 1] +
1099  p[2] * MyBase::mm[ 2] + MyBase::mm[ 3]) / w),
1100  static_cast<T0>((p[0] * MyBase::mm[ 4] + p[1] * MyBase::mm[ 5] +
1101  p[2] * MyBase::mm[ 6] + MyBase::mm[ 7]) / w),
1102  static_cast<T0>((p[0] * MyBase::mm[ 8] + p[1] * MyBase::mm[ 9] +
1103  p[2] * MyBase::mm[10] + MyBase::mm[11]) / w));
1104  }
1105 
1106  return Vec3<T0>(0, 0, 0);
1107  }
1108 
1110  template<typename T0>
1112  {
1113  return Vec3<T0>(
1114  static_cast<T0>(v[0] * MyBase::mm[ 0] + v[1] * MyBase::mm[ 4] + v[2] * MyBase::mm[ 8]),
1115  static_cast<T0>(v[0] * MyBase::mm[ 1] + v[1] * MyBase::mm[ 5] + v[2] * MyBase::mm[ 9]),
1116  static_cast<T0>(v[0] * MyBase::mm[ 2] + v[1] * MyBase::mm[ 6] + v[2] * MyBase::mm[10]));
1117  }
1118 
1119 
1120 private:
1121  bool invert(Mat4<T> &inverse, T tolerance) const;
1122 
1123  T det2(const Mat4<T> &a, int i0, int i1, int j0, int j1) const {
1124  int i0row = i0 * 4;
1125  int i1row = i1 * 4;
1126  return a.mm[i0row+j0]*a.mm[i1row+j1] - a.mm[i0row+j1]*a.mm[i1row+j0];
1127  }
1128 
1129  T det3(const Mat4<T> &a, int i0, int i1, int i2,
1130  int j0, int j1, int j2) const {
1131  int i0row = i0 * 4;
1132  return a.mm[i0row+j0]*det2(a, i1,i2, j1,j2) +
1133  a.mm[i0row+j1]*det2(a, i1,i2, j2,j0) +
1134  a.mm[i0row+j2]*det2(a, i1,i2, j0,j1);
1135  }
1136 
1137  static const Mat4<T> sIdentity;
1138  static const Mat4<T> sZero;
1139 }; // class Mat4
1140 
1141 
1142 template <typename T>
1143 const Mat4<T> Mat4<T>::sIdentity = Mat4<T>(1, 0, 0, 0,
1144  0, 1, 0, 0,
1145  0, 0, 1, 0,
1146  0, 0, 0, 1);
1147 
1148 template <typename T>
1149 const Mat4<T> Mat4<T>::sZero = Mat4<T>(0, 0, 0, 0,
1150  0, 0, 0, 0,
1151  0, 0, 0, 0,
1152  0, 0, 0, 0);
1153 
1156 template <typename T0, typename T1>
1157 bool operator==(const Mat4<T0> &m0, const Mat4<T1> &m1)
1158 {
1159  const T0 *t0 = m0.asPointer();
1160  const T1 *t1 = m1.asPointer();
1161 
1162  for (int i=0; i<16; ++i) if (!isExactlyEqual(t0[i], t1[i])) return false;
1163  return true;
1164 }
1165 
1168 template <typename T0, typename T1>
1169 bool operator!=(const Mat4<T0> &m0, const Mat4<T1> &m1) { return !(m0 == m1); }
1170 
1173 template <typename S, typename T>
1175 {
1176  return m*scalar;
1177 }
1178 
1181 template <typename S, typename T>
1183 {
1185  result *= scalar;
1186  return result;
1187 }
1188 
1191 template<typename T, typename MT>
1194  const Vec4<T> &_v)
1195 {
1196  MT const *m = _m.asPointer();
1198  _v[0]*m[0] + _v[1]*m[1] + _v[2]*m[2] + _v[3]*m[3],
1199  _v[0]*m[4] + _v[1]*m[5] + _v[2]*m[6] + _v[3]*m[7],
1200  _v[0]*m[8] + _v[1]*m[9] + _v[2]*m[10] + _v[3]*m[11],
1201  _v[0]*m[12] + _v[1]*m[13] + _v[2]*m[14] + _v[3]*m[15]);
1202 }
1203 
1206 template<typename T, typename MT>
1208 operator*(const Vec4<T> &_v,
1209  const Mat4<MT> &_m)
1210 {
1211  MT const *m = _m.asPointer();
1213  _v[0]*m[0] + _v[1]*m[4] + _v[2]*m[8] + _v[3]*m[12],
1214  _v[0]*m[1] + _v[1]*m[5] + _v[2]*m[9] + _v[3]*m[13],
1215  _v[0]*m[2] + _v[1]*m[6] + _v[2]*m[10] + _v[3]*m[14],
1216  _v[0]*m[3] + _v[1]*m[7] + _v[2]*m[11] + _v[3]*m[15]);
1217 }
1218 
1222 template<typename T, typename MT>
1225  const Vec3<T> &_v)
1226 {
1227  MT const *m = _m.asPointer();
1229  _v[0]*m[0] + _v[1]*m[1] + _v[2]*m[2] + m[3],
1230  _v[0]*m[4] + _v[1]*m[5] + _v[2]*m[6] + m[7],
1231  _v[0]*m[8] + _v[1]*m[9] + _v[2]*m[10] + m[11]);
1232 }
1233 
1237 template<typename T, typename MT>
1239 operator*(const Vec3<T> &_v,
1240  const Mat4<MT> &_m)
1241 {
1242  MT const *m = _m.asPointer();
1244  _v[0]*m[0] + _v[1]*m[4] + _v[2]*m[8] + m[12],
1245  _v[0]*m[1] + _v[1]*m[5] + _v[2]*m[9] + m[13],
1246  _v[0]*m[2] + _v[1]*m[6] + _v[2]*m[10] + m[14]);
1247 }
1248 
1251 template <typename T0, typename T1>
1253 operator+(const Mat4<T0> &m0, const Mat4<T1> &m1)
1254 {
1256  result += m1;
1257  return result;
1258 }
1259 
1262 template <typename T0, typename T1>
1264 operator-(const Mat4<T0> &m0, const Mat4<T1> &m1)
1265 {
1267  result -= m1;
1268  return result;
1269 }
1270 
1274 template <typename T0, typename T1>
1276 operator*(const Mat4<T0> &m0, const Mat4<T1> &m1)
1277 {
1279  result *= m1;
1280  return result;
1281 }
1282 
1283 
1287 template<typename T0, typename T1>
1289 {
1290  return Vec3<T1>(
1291  static_cast<T1>(m[0][0]*n[0] + m[0][1]*n[1] + m[0][2]*n[2]),
1292  static_cast<T1>(m[1][0]*n[0] + m[1][1]*n[1] + m[1][2]*n[2]),
1293  static_cast<T1>(m[2][0]*n[0] + m[2][1]*n[1] + m[2][2]*n[2]));
1294 }
1295 
1296 
1298 template<typename T>
1299 bool Mat4<T>::invert(Mat4<T> &inverse, T tolerance) const
1300 {
1301  Mat4<T> temp(*this);
1302  inverse.setIdentity();
1303 
1304  // Forward elimination step
1305  double det = 1.0;
1306  for (int i = 0; i < 4; ++i) {
1307  int row = i;
1308  double max = fabs(temp[i][i]);
1309 
1310  for (int k = i+1; k < 4; ++k) {
1311  if (fabs(temp[k][i]) > max) {
1312  row = k;
1313  max = fabs(temp[k][i]);
1314  }
1315  }
1316 
1317  if (isExactlyEqual(max, 0.0)) return false;
1318 
1319  // must move pivot to row i
1320  if (row != i) {
1321  det = -det;
1322  for (int k = 0; k < 4; ++k) {
1323  std::swap(temp[row][k], temp[i][k]);
1324  std::swap(inverse[row][k], inverse[i][k]);
1325  }
1326  }
1327 
1328  double pivot = temp[i][i];
1329  det *= pivot;
1330 
1331  // scale row i
1332  for (int k = 0; k < 4; ++k) {
1333  temp[i][k] /= pivot;
1334  inverse[i][k] /= pivot;
1335  }
1336 
1337  // eliminate in rows below i
1338  for (int j = i+1; j < 4; ++j) {
1339  double t = temp[j][i];
1340  if (!isExactlyEqual(t, 0.0)) {
1341  // subtract scaled row i from row j
1342  for (int k = 0; k < 4; ++k) {
1343  temp[j][k] -= temp[i][k] * t;
1344  inverse[j][k] -= inverse[i][k] * t;
1345  }
1346  }
1347  }
1348  }
1349 
1350  // Backward elimination step
1351  for (int i = 3; i > 0; --i) {
1352  for (int j = 0; j < i; ++j) {
1353  double t = temp[j][i];
1354 
1355  if (!isExactlyEqual(t, 0.0)) {
1356  for (int k = 0; k < 4; ++k) {
1357  inverse[j][k] -= inverse[i][k]*t;
1358  }
1359  }
1360  }
1361  }
1362  return det*det >= tolerance*tolerance;
1363 }
1364 
1365 template <typename T>
1366 inline bool isAffine(const Mat4<T>& m) {
1367  return (m.col(3) == Vec4<T>(0, 0, 0, 1));
1368 }
1369 
1370 template <typename T>
1371 inline bool hasTranslation(const Mat4<T>& m) {
1372  return (m.row(3) != Vec4<T>(0, 0, 0, 1));
1373 }
1374 
1375 
1378 
1379 #if DWREAL_IS_DOUBLE == 1
1380 typedef Mat4d Mat4f;
1381 #else
1382 typedef Mat4s Mat4f;
1383 #endif // DWREAL_IS_DOUBLE
1384 
1385 } // namespace math
1386 
1387 
1388 template<> inline math::Mat4s zeroVal<math::Mat4s>() { return math::Mat4s::identity(); }
1389 template<> inline math::Mat4d zeroVal<math::Mat4d>() { return math::Mat4d::identity(); }
1390 
1391 } // namespace OPENVDB_VERSION_NAME
1392 } // namespace openvdb
1393 
1394 #endif // OPENVDB_UTIL_MAT4_H_HAS_BEEN_INCLUDED
1395 
1396 // Copyright (c) 2012-2016 DreamWorks Animation LLC
1397 // All rights reserved. This software is distributed under the
1398 // Mozilla Public License 2.0 ( http://www.mozilla.org/MPL/2.0/ )
void setZero()
Definition: Mat4.h:286
void postShear(Axis axis0, Axis axis1, T shear)
Right multiplies a shearing transformation into the matrix.
Definition: Mat4.h:1025
void preShear(Axis axis0, Axis axis1, T shear)
Left multiplies a shearing transformation into the matrix.
Definition: Mat4.h:1010
Mat4< typename promote< T0, T1 >::type > operator-(const Mat4< T0 > &m0, const Mat4< T1 > &m1)
Returns M, where for .
Definition: Mat4.h:1264
T dot(const Vec3< T > &v) const
Dot product.
Definition: Vec3.h:203
const Mat4 & operator=(const Mat4< Source > &m)
Assignment operator.
Definition: Mat4.h:365
Mat4(const Vec4< Source > &v1, const Vec4< Source > &v2, const Vec4< Source > &v3, const Vec4< Source > &v4, bool rows=true)
Definition: Mat4.h:120
#define OPENVDB_DEPRECATED
Definition: Platform.h:49
T & z()
Definition: Vec3.h:99
bool isExactlyEqual(const T0 &a, const T1 &b)
Return true if a is exactly equal to b.
Definition: Math.h:407
void postScale(const Vec3< T0 > &v)
Definition: Mat4.h:800
void setToRotation(Axis axis, T angle)
Sets the matrix to a rotation about the given axis.
Definition: Mat4.h:825
void setRow(int i, const Vec4< T > &v)
Set ith row to vector v.
Definition: Mat4.h:162
General-purpose arithmetic and comparison routines, most of which accept arbitrary value types (or at...
Vec3< T0 > pretransformH(const Vec3< T0 > &p) const
Transform a Vec3 by pre-multiplication, doing homogenous division.
Definition: Mat4.h:1090
Definition: Math.h:857
const T * operator[](int i) const
Definition: Mat4.h:200
const Mat4< T > & operator*=(const Mat4< S > &m1)
Return m, where for .
Definition: Mat4.h:481
static const Mat4< T > & identity()
Predefined constant for identity matrix.
Definition: Mat4.h:152
Definition: Mat.h:146
Mat4(const Mat< 4, T > &m)
Copy constructor.
Definition: Mat4.h:131
#define OPENVDB_THROW(exception, message)
Definition: Exceptions.h:97
T mm[SIZE *SIZE]
Definition: Mat.h:141
T operator()(int i, int j) const
Definition: Mat4.h:220
bool eq(const Mat4 &m, T eps=1.0e-8) const
Test if "this" is equivalent to m with tolerance of eps value.
Definition: Mat4.h:375
Definition: Exceptions.h:78
T * asPointer()
Direct access to the internal data.
Definition: Mat4.h:204
void setIdentity()
Set "this" matrix to identity.
Definition: Mat4.h:307
Vec3< typename promote< T, MT >::type > operator*(const Vec3< T > &_v, const Mat4< MT > &_m)
Returns v, where for .
Definition: Mat4.h:1239
T det() const
Determinant of matrix.
Definition: Mat4.h:679
Mat4< typename promote< S, T >::type > operator*(S scalar, const Mat4< T > &m)
Returns M, where for .
Definition: Mat4.h:1174
Definition: Math.h:858
MatType shear(Axis axis0, Axis axis1, typename MatType::value_type shear)
Set the matrix to a shear along axis0 by a fraction of axis1.
Definition: Mat.h:667
void setColumns(const Vec4< T > &v1, const Vec4< T > &v2, const Vec4< T > &v3, const Vec4< T > &v4)
Set the columns of "this" matrix to the vectors v1, v2, v3, v4.
Definition: Mat4.h:253
T ValueType
Definition: Mat4.h:62
OPENVDB_DEPRECATED void setBasis(const Vec4< T > &v1, const Vec4< T > &v2, const Vec4< T > &v3, const Vec4< T > &v4)
Set the rows of "this" matrix to the vectors v1, v2, v3, v4.
Definition: Mat4.h:278
const boost::disable_if_c< VecTraits< T >::IsVec, T >::type & max(const T &a, const T &b)
Definition: Composite.h:132
Vec3< T0 > transformH(const Vec3< T0 > &p) const
Transform a Vec3 by post-multiplication, doing homogenous divison.
Definition: Mat4.h:1068
Mat4s Mat4f
Definition: Mat4.h:1382
void preRotate(Axis axis, T angle)
Left multiplies by a rotation clock-wiseabout the given axis into this matrix.
Definition: Mat4.h:840
4x4 -matrix class.
Definition: Mat3.h:48
void preTranslate(const Vec3< T0 > &tr)
Left multiples by the specified translation, i.e. Trans * (*this)
Definition: Mat4.h:745
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:370
Vec4< typename promote< T, MT >::type > operator*(const Vec4< T > &_v, const Mat4< MT > &_m)
Returns v, where for .
Definition: Mat4.h:1208
const Mat4< T > & operator+=(const Mat4< S > &m1)
Returns m0, where for .
Definition: Mat4.h:423
Mat4(Source *a)
Constructor given array of elements, the ordering is in row major form:
Definition: Mat4.h:76
Vec3< T0 > pretransform(const Vec3< T0 > &v) const
Transform a Vec3 by pre-multiplication, without homogenous division.
Definition: Mat4.h:1061
void setToScale(const Vec3< T0 > &v)
Sets the matrix to a matrix that scales by v.
Definition: Mat4.h:768
Mat4 inverse(T tolerance=0) const
Definition: Mat4.h:527
#define OPENVDB_VERSION_NAME
Definition: version.h:43
Mat4 transpose() const
Definition: Mat4.h:514
void setMat3(const Mat3< T > &m)
Set upper left to a Mat3.
Definition: Mat4.h:332
bool operator==(const Mat4< T0 > &m0, const Mat4< T1 > &m1)
Equality operator, does exact floating point comparisons.
Definition: Mat4.h:1157
void postRotate(Axis axis, T angle)
Right multiplies by a rotation clock-wiseabout the given axis into this matrix.
Definition: Mat4.h:920
Mat4(Source a, Source b, Source c, Source d, Source e, Source f, Source g, Source h, Source i, Source j, Source k, Source l, Source m, Source n, Source o, Source p)
Constructor given array of elements, the ordering is in row major form:
Definition: Mat4.h:91
T & x()
Reference to the component, e.g. v.x() = 4.5f;.
Definition: Vec3.h:97
bool isAffine(const Mat4< T > &m)
Definition: Mat4.h:1366
Definition: Exceptions.h:39
Vec3< T0 > transform(const Vec3< T0 > &v) const
Transform a Vec3 by post-multiplication, without homogenous division.
Definition: Mat4.h:1047
Vec4< T > row(int i) const
Get ith row, e.g. Vec4f v = m.row(1);.
Definition: Mat4.h:173
3x3 matrix class.
Definition: Mat3.h:54
void setToShear(Axis axis0, Axis axis1, T shearby)
Sets the matrix to a shear along axis0 by a fraction of axis1.
Definition: Mat4.h:1002
T value_type
Data type held by the matrix.
Definition: Mat4.h:61
const Mat4< T > & operator-=(const Mat4< S > &m1)
Returns m0, where for .
Definition: Mat4.h:452
T angle(const Vec2< T > &v1, const Vec2< T > &v2)
Definition: Vec2.h:446
Mat4< float > Mat4s
Definition: Mat4.h:1376
Mat3< T > getMat3() const
Definition: Mat4.h:339
T & operator()(int i, int j)
Definition: Mat4.h:210
void setToRotation(const Vec3< T > &v1, const Vec3< T > &v2)
Sets the matrix to a rotation that maps v1 onto v2 about the cross product of v1 and v2...
Definition: Mat4.h:834
Mat< 4, T > MyBase
Definition: Mat4.h:63
const Mat4< T > & operator*=(S scalar)
Return m, where for .
Definition: Mat4.h:397
Definition: Math.h:859
void setTranslation(const Vec3< T > &t)
Definition: Mat4.h:356
Mat4< double > Mat4d
Definition: Mat4.h:1377
Vec4< T0 > pretransform(const Vec4< T0 > &v) const
Transform a Vec4 by pre-multiplication.
Definition: Mat4.h:1054
T & y()
Definition: Vec3.h:98
void setToRotation(const Vec3< T > &axis, T angle)
Sets the matrix to a rotation about an arbitrary axis.
Definition: Mat4.h:830
static Mat4 translation(const Vec3d &v)
Sets the matrix to a matrix that translates by v.
Definition: Mat4.h:709
static const Mat4< T > & zero()
Predefined constant for zero matrix.
Definition: Mat4.h:157
bool hasTranslation(const Mat4< T > &m)
Definition: Mat4.h:1371
const T * asPointer() const
Definition: Mat4.h:205
Mat4()
Trivial constructor, the matrix is NOT initialized.
Definition: Mat4.h:66
Vec4< typename promote< T, MT >::type > operator*(const Mat4< MT > &_m, const Vec4< T > &_v)
Returns v, where for .
Definition: Mat4.h:1193
Vec3< T0 > transform3x3(const Vec3< T0 > &v) const
Transform a Vec3 by post-multiplication, without translation.
Definition: Mat4.h:1111
Mat4< typename promote< T0, T1 >::type > operator*(const Mat4< T0 > &m0, const Mat4< T1 > &m1)
Returns M, where for .
Definition: Mat4.h:1276
Vec3< typename promote< T, MT >::type > operator*(const Mat4< MT > &_m, const Vec3< T > &_v)
Returns v, where for .
Definition: Mat4.h:1224
Mat4< typename promote< T0, T1 >::type > operator+(const Mat4< T0 > &m0, const Mat4< T1 > &m1)
Returns M, where for .
Definition: Mat4.h:1253
void setToTranslation(const Vec3< T0 > &v)
Sets the matrix to a matrix that translates by v.
Definition: Mat4.h:720
Axis
Definition: Math.h:856
void setCol(int j, const Vec4< T > &v)
Set jth column to vector v.
Definition: Mat4.h:180
Mat4< typename promote< S, T >::type > operator*(const Mat4< T > &m, S scalar)
Returns M, where for .
Definition: Mat4.h:1182
T det() const
Determinant of matrix.
Definition: Mat3.h:523
Vec4< T > col(int j) const
Get jth column, e.g. Vec4f v = m.col(0);.
Definition: Mat4.h:190
Mat4< T > operator-() const
Negation operator, for e.g. m1 = -m2;.
Definition: Mat4.h:385
Vec4< T0 > transform(const Vec4< T0 > &v) const
Transform a Vec4 by post-multiplication.
Definition: Mat4.h:1040
#define OPENVDB_USE_VERSION_NAMESPACE
Definition: version.h:71
Definition: Mat.h:52
void setRows(const Vec4< T > &v1, const Vec4< T > &v2, const Vec4< T > &v3, const Vec4< T > &v4)
Set the rows of "this" matrix to the vectors v1, v2, v3, v4.
Definition: Mat4.h:228
Vec3< T > getTranslation() const
Return the translation component.
Definition: Mat4.h:351
Mat4(const Mat4< Source > &m)
Conversion constructor.
Definition: Mat4.h:142
void preScale(const Vec3< T0 > &v)
Definition: Mat4.h:778
void postTranslate(const Vec3< T0 > &tr)
Right multiplies by the specified translation matrix, i.e. (*this) * Trans.
Definition: Mat4.h:756
Definition: Mat4.h:51
bool operator!=(const Mat4< T0 > &m0, const Mat4< T1 > &m1)
Inequality operator, does exact floating point comparisons.
Definition: Mat4.h:1169
T * asPointer()
Definition: Mat3.h:208
Vec3< T1 > transformNormal(const Mat4< T0 > &m, const Vec3< T1 > &n)
Definition: Mat4.h:1288