31 #ifndef OPENVDB_MATH_MAT4_H_HAS_BEEN_INCLUDED
32 #define OPENVDB_MATH_MAT4_H_HAS_BEEN_INCLUDED
34 #include <openvdb/Exceptions.h>
35 #include <openvdb/Platform.h>
51 template<
typename T>
class Vec4;
75 template<
typename Source>
80 for (i = 0; i < 16; i++) {
92 template<
typename Source>
93 Mat4(Source a, Source b, Source c, Source d,
94 Source e, Source f, Source g, Source h,
95 Source i, Source j, Source k, Source l,
96 Source m, Source n, Source o, Source p)
120 template<
typename Source>
124 setBasis(v1, v2, v3, v4);
130 for (
int i = 0; i < 4; ++i) {
131 for (
int j = 0; j < 4; ++j) {
132 MyBase::mm[i*4 + j] = m[i][j];
138 template<
typename Source>
143 for (
int i=0; i<16; ++i) {
144 MyBase::mm[i] =
static_cast<T
>(src[i]);
163 MyBase::mm[i4+0] = v[0];
164 MyBase::mm[i4+1] = v[1];
165 MyBase::mm[i4+2] = v[2];
166 MyBase::mm[i4+3] = v[3];
173 return Vec4<T>((*this)(i,0), (*
this)(i,1), (*
this)(i,2), (*
this)(i,3));
180 MyBase::mm[ 0+j] = v[0];
181 MyBase::mm[ 4+j] = v[1];
182 MyBase::mm[ 8+j] = v[2];
183 MyBase::mm[12+j] = v[3];
190 return Vec4<T>((*this)(0,j), (*
this)(1,j), (*
this)(2,j), (*
this)(3,j));
197 const T*
operator[](
int i)
const {
return &(MyBase::mm[i<<2]); }
207 T& operator()(
int i,
int j)
211 return MyBase::mm[4*i+j];
217 T operator()(
int i,
int j)
const
221 return MyBase::mm[4*i+j];
228 MyBase::mm[ 0] = v1[0];
229 MyBase::mm[ 1] = v1[1];
230 MyBase::mm[ 2] = v1[2];
231 MyBase::mm[ 3] = v1[3];
233 MyBase::mm[ 4] = v2[0];
234 MyBase::mm[ 5] = v2[1];
235 MyBase::mm[ 6] = v2[2];
236 MyBase::mm[ 7] = v2[3];
238 MyBase::mm[ 8] = v3[0];
239 MyBase::mm[ 9] = v3[1];
240 MyBase::mm[10] = v3[2];
241 MyBase::mm[11] = v3[3];
243 MyBase::mm[12] = v4[0];
244 MyBase::mm[13] = v4[1];
245 MyBase::mm[14] = v4[2];
246 MyBase::mm[15] = v4[3];
299 for (
int i = 0; i < 3; i++)
300 for (
int j=0; j < 3; j++)
301 MyBase::mm[i*4+j] = m[i][j];
308 for (
int i = 0; i < 3; i++)
309 for (
int j = 0; j < 3; j++)
310 m[i][j] = MyBase::mm[i*4+j];
318 return Vec3<T>(MyBase::mm[12], MyBase::mm[13], MyBase::mm[14]);
323 MyBase::mm[12] = t[0];
324 MyBase::mm[13] = t[1];
325 MyBase::mm[14] = t[2];
329 template<
typename Source>
335 std::copy(src, (src + this->numElements()), MyBase::mm);
340 bool eq(
const Mat4 &m, T eps=1.0e-8)
const
342 for (
int i = 0; i < 16; i++) {
353 -MyBase::mm[ 0], -MyBase::mm[ 1], -MyBase::mm[ 2], -MyBase::mm[ 3],
354 -MyBase::mm[ 4], -MyBase::mm[ 5], -MyBase::mm[ 6], -MyBase::mm[ 7],
355 -MyBase::mm[ 8], -MyBase::mm[ 9], -MyBase::mm[10], -MyBase::mm[11],
356 -MyBase::mm[12], -MyBase::mm[13], -MyBase::mm[14], -MyBase::mm[15]
361 template <
typename S>
364 MyBase::mm[ 0] *= scalar;
365 MyBase::mm[ 1] *= scalar;
366 MyBase::mm[ 2] *= scalar;
367 MyBase::mm[ 3] *= scalar;
369 MyBase::mm[ 4] *= scalar;
370 MyBase::mm[ 5] *= scalar;
371 MyBase::mm[ 6] *= scalar;
372 MyBase::mm[ 7] *= scalar;
374 MyBase::mm[ 8] *= scalar;
375 MyBase::mm[ 9] *= scalar;
376 MyBase::mm[10] *= scalar;
377 MyBase::mm[11] *= scalar;
379 MyBase::mm[12] *= scalar;
380 MyBase::mm[13] *= scalar;
381 MyBase::mm[14] *= scalar;
382 MyBase::mm[15] *= scalar;
387 template <
typename S>
392 MyBase::mm[ 0] += s[ 0];
393 MyBase::mm[ 1] += s[ 1];
394 MyBase::mm[ 2] += s[ 2];
395 MyBase::mm[ 3] += s[ 3];
397 MyBase::mm[ 4] += s[ 4];
398 MyBase::mm[ 5] += s[ 5];
399 MyBase::mm[ 6] += s[ 6];
400 MyBase::mm[ 7] += s[ 7];
402 MyBase::mm[ 8] += s[ 8];
403 MyBase::mm[ 9] += s[ 9];
404 MyBase::mm[10] += s[10];
405 MyBase::mm[11] += s[11];
407 MyBase::mm[12] += s[12];
408 MyBase::mm[13] += s[13];
409 MyBase::mm[14] += s[14];
410 MyBase::mm[15] += s[15];
416 template <
typename S>
421 MyBase::mm[ 0] -= s[ 0];
422 MyBase::mm[ 1] -= s[ 1];
423 MyBase::mm[ 2] -= s[ 2];
424 MyBase::mm[ 3] -= s[ 3];
426 MyBase::mm[ 4] -= s[ 4];
427 MyBase::mm[ 5] -= s[ 5];
428 MyBase::mm[ 6] -= s[ 6];
429 MyBase::mm[ 7] -= s[ 7];
431 MyBase::mm[ 8] -= s[ 8];
432 MyBase::mm[ 9] -= s[ 9];
433 MyBase::mm[10] -= s[10];
434 MyBase::mm[11] -= s[11];
436 MyBase::mm[12] -= s[12];
437 MyBase::mm[13] -= s[13];
438 MyBase::mm[14] -= s[14];
439 MyBase::mm[15] -= s[15];
445 template <
typename S>
453 for (
int i = 0; i < 4; i++) {
455 MyBase::mm[i4+0] =
static_cast<T
>(s0[i4+0] * s1[ 0] +
460 MyBase::mm[i4+1] =
static_cast<T
>(s0[i4+0] * s1[ 1] +
465 MyBase::mm[i4+2] =
static_cast<T
>(s0[i4+0] * s1[ 2] +
470 MyBase::mm[i4+3] =
static_cast<T
>(s0[i4+0] * s1[ 3] +
482 MyBase::mm[ 0], MyBase::mm[ 4], MyBase::mm[ 8], MyBase::mm[12],
483 MyBase::mm[ 1], MyBase::mm[ 5], MyBase::mm[ 9], MyBase::mm[13],
484 MyBase::mm[ 2], MyBase::mm[ 6], MyBase::mm[10], MyBase::mm[14],
485 MyBase::mm[ 3], MyBase::mm[ 7], MyBase::mm[11], MyBase::mm[15]
515 T m0011 = m[0][0] * m[1][1];
516 T m0012 = m[0][0] * m[1][2];
517 T m0110 = m[0][1] * m[1][0];
518 T m0210 = m[0][2] * m[1][0];
519 T m0120 = m[0][1] * m[2][0];
520 T m0220 = m[0][2] * m[2][0];
522 T detA = m0011 * m[2][2] - m0012 * m[2][1] - m0110 * m[2][2]
523 + m0210 * m[2][1] + m0120 * m[1][2] - m0220 * m[1][1];
525 bool hasPerspective =
532 if (hasPerspective) {
533 det = m[0][3] * det3(m, 1,2,3, 0,2,1)
534 + m[1][3] * det3(m, 2,0,3, 0,2,1)
535 + m[2][3] * det3(m, 3,0,1, 0,2,1)
538 det = detA * m[3][3];
558 inv[0][0] = detA * ( m[1][1] * m[2][2] - m[1][2] * m[2][1]);
559 inv[0][1] = detA * (-m[0][1] * m[2][2] + m[0][2] * m[2][1]);
560 inv[0][2] = detA * ( m[0][1] * m[1][2] - m[0][2] * m[1][1]);
562 inv[1][0] = detA * (-m[1][0] * m[2][2] + m[1][2] * m[2][0]);
563 inv[1][1] = detA * ( m[0][0] * m[2][2] - m0220);
564 inv[1][2] = detA * ( m0210 - m0012);
566 inv[2][0] = detA * ( m[1][0] * m[2][1] - m[1][1] * m[2][0]);
567 inv[2][1] = detA * ( m0120 - m[0][0] * m[2][1]);
568 inv[2][2] = detA * ( m0011 - m0110);
570 if (hasPerspective) {
575 r[0] = m[3][0] * inv[0][0] + m[3][1] * inv[1][0]
576 + m[3][2] * inv[2][0];
577 r[1] = m[3][0] * inv[0][1] + m[3][1] * inv[1][1]
578 + m[3][2] * inv[2][1];
579 r[2] = m[3][0] * inv[0][2] + m[3][1] * inv[1][2]
580 + m[3][2] * inv[2][2];
583 p[0] = inv[0][0] * m[0][3] + inv[0][1] * m[1][3]
584 + inv[0][2] * m[2][3];
585 p[1] = inv[1][0] * m[0][3] + inv[1][1] * m[1][3]
586 + inv[1][2] * m[2][3];
587 p[2] = inv[2][0] * m[0][3] + inv[2][1] * m[1][3]
588 + inv[2][2] * m[2][3];
590 T h = m[3][3] - p.
dot(
Vec3<T>(m[3][0],m[3][1],m[3][2]));
601 inv[3][0] = -h * r[0];
602 inv[3][1] = -h * r[1];
603 inv[3][2] = -h * r[2];
605 inv[0][3] = -h * p[0];
606 inv[1][3] = -h * p[1];
607 inv[2][3] = -h * p[2];
613 inv[0][0] += p[0] * r[0];
614 inv[0][1] += p[0] * r[1];
615 inv[0][2] += p[0] * r[2];
616 inv[1][0] += p[1] * r[0];
617 inv[1][1] += p[1] * r[1];
618 inv[1][2] += p[1] * r[2];
619 inv[2][0] += p[2] * r[0];
620 inv[2][1] += p[2] * r[1];
621 inv[2][2] += p[2] * r[2];
625 inv[3][0] = - (m[3][0] * inv[0][0] + m[3][1] * inv[1][0]
626 + m[3][2] * inv[2][0]);
627 inv[3][1] = - (m[3][0] * inv[0][1] + m[3][1] * inv[1][1]
628 + m[3][2] * inv[2][1]);
629 inv[3][2] = - (m[3][0] * inv[0][2] + m[3][1] * inv[1][2]
630 + m[3][2] * inv[2][2]);
654 for (i = 0; i < 4; i++) {
655 ap = &MyBase::mm[ 0];
657 for (j = 0; j < 4; j++) {
658 for (k = 0; k < 4; k++) {
659 if ((k != i) && (j != 0)) {
666 det += sign * MyBase::mm[i] * submat.
det();
678 {
return snapBasis(*
this, axis, direction);}
684 T(1), T(0), T(0), T(0),
685 T(0), T(1), T(0), T(0),
686 T(0), T(0), T(1), T(0),
687 T(v.
x()), T(v.
y()),T(v.
z()), T(1));
691 template <
typename T0>
709 MyBase::mm[12] = v.
x();
710 MyBase::mm[13] = v.
y();
711 MyBase::mm[14] = v.
z();
716 template <
typename T0>
722 *
this = Tr * (*this);
727 template <
typename T0>
733 *
this = (*this) * Tr;
739 template <
typename T0>
743 MyBase::mm[ 0] = v.
x();
744 MyBase::mm[ 5] = v.
y();
745 MyBase::mm[10] = v.
z();
749 template <
typename T0>
752 MyBase::mm[ 0] *= v.
x();
753 MyBase::mm[ 1] *= v.
x();
754 MyBase::mm[ 2] *= v.
x();
755 MyBase::mm[ 3] *= v.
x();
757 MyBase::mm[ 4] *= v.
y();
758 MyBase::mm[ 5] *= v.
y();
759 MyBase::mm[ 6] *= v.
y();
760 MyBase::mm[ 7] *= v.
y();
762 MyBase::mm[ 8] *= v.
z();
763 MyBase::mm[ 9] *= v.
z();
764 MyBase::mm[10] *= v.
z();
765 MyBase::mm[11] *= v.
z();
771 template <
typename T0>
775 MyBase::mm[ 0] *= v.
x();
776 MyBase::mm[ 1] *= v.
y();
777 MyBase::mm[ 2] *= v.
z();
779 MyBase::mm[ 4] *= v.
x();
780 MyBase::mm[ 5] *= v.
y();
781 MyBase::mm[ 6] *= v.
z();
783 MyBase::mm[ 8] *= v.
x();
784 MyBase::mm[ 9] *= v.
y();
785 MyBase::mm[10] *= v.
z();
787 MyBase::mm[12] *= v.
x();
788 MyBase::mm[13] *= v.
y();
789 MyBase::mm[14] *= v.
z();
814 T c =
static_cast<T
>(cos(angle));
815 T s = -
static_cast<T
>(sin(angle));
822 a4 = c * MyBase::mm[ 4] - s * MyBase::mm[ 8];
823 a5 = c * MyBase::mm[ 5] - s * MyBase::mm[ 9];
824 a6 = c * MyBase::mm[ 6] - s * MyBase::mm[10];
825 a7 = c * MyBase::mm[ 7] - s * MyBase::mm[11];
828 MyBase::mm[ 8] = s * MyBase::mm[ 4] + c * MyBase::mm[ 8];
829 MyBase::mm[ 9] = s * MyBase::mm[ 5] + c * MyBase::mm[ 9];
830 MyBase::mm[10] = s * MyBase::mm[ 6] + c * MyBase::mm[10];
831 MyBase::mm[11] = s * MyBase::mm[ 7] + c * MyBase::mm[11];
844 a0 = c * MyBase::mm[ 0] + s * MyBase::mm[ 8];
845 a1 = c * MyBase::mm[ 1] + s * MyBase::mm[ 9];
846 a2 = c * MyBase::mm[ 2] + s * MyBase::mm[10];
847 a3 = c * MyBase::mm[ 3] + s * MyBase::mm[11];
849 MyBase::mm[ 8] = -s * MyBase::mm[ 0] + c * MyBase::mm[ 8];
850 MyBase::mm[ 9] = -s * MyBase::mm[ 1] + c * MyBase::mm[ 9];
851 MyBase::mm[10] = -s * MyBase::mm[ 2] + c * MyBase::mm[10];
852 MyBase::mm[11] = -s * MyBase::mm[ 3] + c * MyBase::mm[11];
866 a0 = c * MyBase::mm[ 0] - s * MyBase::mm[ 4];
867 a1 = c * MyBase::mm[ 1] - s * MyBase::mm[ 5];
868 a2 = c * MyBase::mm[ 2] - s * MyBase::mm[ 6];
869 a3 = c * MyBase::mm[ 3] - s * MyBase::mm[ 7];
871 MyBase::mm[ 4] = s * MyBase::mm[ 0] + c * MyBase::mm[ 4];
872 MyBase::mm[ 5] = s * MyBase::mm[ 1] + c * MyBase::mm[ 5];
873 MyBase::mm[ 6] = s * MyBase::mm[ 2] + c * MyBase::mm[ 6];
874 MyBase::mm[ 7] = s * MyBase::mm[ 3] + c * MyBase::mm[ 7];
894 T c =
static_cast<T
>(cos(angle));
895 T s = -
static_cast<T
>(sin(angle));
904 a2 = c * MyBase::mm[ 2] - s * MyBase::mm[ 1];
905 a6 = c * MyBase::mm[ 6] - s * MyBase::mm[ 5];
906 a10 = c * MyBase::mm[10] - s * MyBase::mm[ 9];
907 a14 = c * MyBase::mm[14] - s * MyBase::mm[13];
910 MyBase::mm[ 1] = c * MyBase::mm[ 1] + s * MyBase::mm[ 2];
911 MyBase::mm[ 5] = c * MyBase::mm[ 5] + s * MyBase::mm[ 6];
912 MyBase::mm[ 9] = c * MyBase::mm[ 9] + s * MyBase::mm[10];
913 MyBase::mm[13] = c * MyBase::mm[13] + s * MyBase::mm[14];
917 MyBase::mm[10] = a10;
918 MyBase::mm[14] = a14;
926 a2 = c * MyBase::mm[ 2] + s * MyBase::mm[ 0];
927 a6 = c * MyBase::mm[ 6] + s * MyBase::mm[ 4];
928 a10 = c * MyBase::mm[10] + s * MyBase::mm[ 8];
929 a14 = c * MyBase::mm[14] + s * MyBase::mm[12];
931 MyBase::mm[ 0] = c * MyBase::mm[ 0] - s * MyBase::mm[ 2];
932 MyBase::mm[ 4] = c * MyBase::mm[ 4] - s * MyBase::mm[ 6];
933 MyBase::mm[ 8] = c * MyBase::mm[ 8] - s * MyBase::mm[10];
934 MyBase::mm[12] = c * MyBase::mm[12] - s * MyBase::mm[14];
938 MyBase::mm[10] = a10;
939 MyBase::mm[14] = a14;
947 a1 = c * MyBase::mm[ 1] - s * MyBase::mm[ 0];
948 a5 = c * MyBase::mm[ 5] - s * MyBase::mm[ 4];
949 a9 = c * MyBase::mm[ 9] - s * MyBase::mm[ 8];
950 a13 = c * MyBase::mm[13] - s * MyBase::mm[12];
952 MyBase::mm[ 0] = c * MyBase::mm[ 0] + s * MyBase::mm[ 1];
953 MyBase::mm[ 4] = c * MyBase::mm[ 4] + s * MyBase::mm[ 5];
954 MyBase::mm[ 8] = c * MyBase::mm[ 8] + s * MyBase::mm[ 9];
955 MyBase::mm[12] = c * MyBase::mm[12] + s * MyBase::mm[13];
960 MyBase::mm[13] = a13;
976 *
this = shear<Mat4<T> >(axis0, axis1, shearby);
984 int index0 =
static_cast<int>(axis0);
985 int index1 =
static_cast<int>(axis1);
988 MyBase::mm[index1 * 4 + 0] += shear * MyBase::mm[index0 * 4 + 0];
989 MyBase::mm[index1 * 4 + 1] += shear * MyBase::mm[index0 * 4 + 1];
990 MyBase::mm[index1 * 4 + 2] += shear * MyBase::mm[index0 * 4 + 2];
991 MyBase::mm[index1 * 4 + 3] += shear * MyBase::mm[index0 * 4 + 3];
999 int index0 =
static_cast<int>(axis0);
1000 int index1 =
static_cast<int>(axis1);
1003 MyBase::mm[index0 + 0] += shear * MyBase::mm[index1 + 0];
1004 MyBase::mm[index0 + 4] += shear * MyBase::mm[index1 + 4];
1005 MyBase::mm[index0 + 8] += shear * MyBase::mm[index1 + 8];
1006 MyBase::mm[index0 + 12] += shear * MyBase::mm[index1 + 12];
1011 template<
typename T0>
1014 return static_cast< Vec4<T0> >(v * *
this);
1018 template<
typename T0>
1021 return static_cast< Vec3<T0> >(v * *
this);
1025 template<
typename T0>
1028 return static_cast< Vec4<T0> >(*
this * v);
1032 template<
typename T0>
1035 return static_cast< Vec3<T0> >(*
this * v);
1039 template<
typename T0>
1045 w = p[0] * MyBase::mm[ 3] + p[1] * MyBase::mm[ 7] + p[2] * MyBase::mm[11] + MyBase::mm[15];
1048 return Vec3<T0>(
static_cast<T0
>((p[0] * MyBase::mm[ 0] + p[1] * MyBase::mm[ 4] +
1049 p[2] * MyBase::mm[ 8] + MyBase::mm[12]) / w),
1050 static_cast<T0
>((p[0] * MyBase::mm[ 1] + p[1] * MyBase::mm[ 5] +
1051 p[2] * MyBase::mm[ 9] + MyBase::mm[13]) / w),
1052 static_cast<T0
>((p[0] * MyBase::mm[ 2] + p[1] * MyBase::mm[ 6] +
1053 p[2] * MyBase::mm[10] + MyBase::mm[14]) / w));
1060 template<
typename T0>
1066 w = p[0] * MyBase::mm[12] + p[1] * MyBase::mm[13] + p[2] * MyBase::mm[14] + MyBase::mm[15];
1069 return Vec3<T0>(
static_cast<T0
>((p[0] * MyBase::mm[ 0] + p[1] * MyBase::mm[ 1] +
1070 p[2] * MyBase::mm[ 2] + MyBase::mm[ 3]) / w),
1071 static_cast<T0
>((p[0] * MyBase::mm[ 4] + p[1] * MyBase::mm[ 5] +
1072 p[2] * MyBase::mm[ 6] + MyBase::mm[ 7]) / w),
1073 static_cast<T0
>((p[0] * MyBase::mm[ 8] + p[1] * MyBase::mm[ 9] +
1074 p[2] * MyBase::mm[10] + MyBase::mm[11]) / w));
1081 template<
typename T0>
1085 static_cast<T0
>(v[0] * MyBase::mm[ 0] + v[1] * MyBase::mm[ 4] + v[2] * MyBase::mm[ 8]),
1086 static_cast<T0>(v[0] * MyBase::mm[ 1] + v[1] * MyBase::mm[ 5] + v[2] * MyBase::mm[ 9]),
1087 static_cast<T0
>(v[0] * MyBase::mm[ 2] + v[1] * MyBase::mm[ 6] + v[2] * MyBase::mm[10]));
1094 T det2(
const Mat4<T> &a,
int i0,
int i1,
int j0,
int j1)
const {
1097 return a.
mm[i0row+j0]*a.
mm[i1row+j1] - a.
mm[i0row+j1]*a.
mm[i1row+j0];
1100 T det3(
const Mat4<T> &a,
int i0,
int i1,
int i2,
1101 int j0,
int j1,
int j2)
const {
1103 return a.mm[i0row+j0]*det2(a, i1,i2, j1,j2) +
1104 a.mm[i0row+j1]*det2(a, i1,i2, j2,j0) +
1105 a.mm[i0row+j2]*det2(a, i1,i2, j0,j1);
1108 static const Mat4<T> sIdentity;
1109 static const Mat4<T> sZero;
1113 template <
typename T>
1114 const Mat4<T> Mat4<T>::sIdentity = Mat4<T>(1, 0, 0, 0,
1119 template <
typename T>
1120 const Mat4<T> Mat4<T>::sZero = Mat4<T>(0, 0, 0, 0,
1127 template <
typename T0,
typename T1>
1133 for (
int i=0; i<16; ++i)
if (!
isExactlyEqual(t0[i], t1[i]))
return false;
1139 template <
typename T0,
typename T1>
1144 template <
typename S,
typename T>
1152 template <
typename S,
typename T>
1162 template<
typename T,
typename MT>
1169 _v[0]*m[0] + _v[1]*m[1] + _v[2]*m[2] + _v[3]*m[3],
1170 _v[0]*m[4] + _v[1]*m[5] + _v[2]*m[6] + _v[3]*m[7],
1171 _v[0]*m[8] + _v[1]*m[9] + _v[2]*m[10] + _v[3]*m[11],
1172 _v[0]*m[12] + _v[1]*m[13] + _v[2]*m[14] + _v[3]*m[15]);
1177 template<
typename T,
typename MT>
1184 _v[0]*m[0] + _v[1]*m[4] + _v[2]*m[8] + _v[3]*m[12],
1185 _v[0]*m[1] + _v[1]*m[5] + _v[2]*m[9] + _v[3]*m[13],
1186 _v[0]*m[2] + _v[1]*m[6] + _v[2]*m[10] + _v[3]*m[14],
1187 _v[0]*m[3] + _v[1]*m[7] + _v[2]*m[11] + _v[3]*m[15]);
1193 template<
typename T,
typename MT>
1200 _v[0]*m[0] + _v[1]*m[1] + _v[2]*m[2] + m[3],
1201 _v[0]*m[4] + _v[1]*m[5] + _v[2]*m[6] + m[7],
1202 _v[0]*m[8] + _v[1]*m[9] + _v[2]*m[10] + m[11]);
1208 template<
typename T,
typename MT>
1215 _v[0]*m[0] + _v[1]*m[4] + _v[2]*m[8] + m[12],
1216 _v[0]*m[1] + _v[1]*m[5] + _v[2]*m[9] + m[13],
1217 _v[0]*m[2] + _v[1]*m[6] + _v[2]*m[10] + m[14]);
1222 template <
typename T0,
typename T1>
1233 template <
typename T0,
typename T1>
1245 template <
typename T0,
typename T1>
1258 template<
typename T0,
typename T1>
1262 static_cast<T1
>(m[0][0]*n[0] + m[0][1]*n[1] + m[0][2]*n[2]),
1263 static_cast<T1>(m[1][0]*n[0] + m[1][1]*n[1] + m[1][2]*n[2]),
1264 static_cast<T1
>(m[2][0]*n[0] + m[2][1]*n[1] + m[2][2]*n[2]));
1269 template<
typename T>
1270 bool Mat4<T>::invert(Mat4<T> &inverse, T tolerance)
const
1272 Mat4<T> temp(*
this);
1273 inverse.setIdentity();
1277 for (
int i = 0; i < 4; ++i) {
1279 double max = fabs(temp[i][i]);
1281 for (
int k = i+1; k < 4; ++k) {
1282 if (fabs(temp[k][i]) > max) {
1284 max = fabs(temp[k][i]);
1293 for (
int k = 0; k < 4; ++k) {
1294 std::swap(temp[row][k], temp[i][k]);
1295 std::swap(inverse[row][k], inverse[i][k]);
1299 double pivot = temp[i][i];
1303 for (
int k = 0; k < 4; ++k) {
1304 temp[i][k] /= pivot;
1305 inverse[i][k] /= pivot;
1309 for (
int j = i+1; j < 4; ++j) {
1310 double t = temp[j][i];
1313 for (
int k = 0; k < 4; ++k) {
1314 temp[j][k] -= temp[i][k] * t;
1315 inverse[j][k] -= inverse[i][k] * t;
1322 for (
int i = 3; i > 0; --i) {
1323 for (
int j = 0; j < i; ++j) {
1324 double t = temp[j][i];
1327 for (
int k = 0; k < 4; ++k) {
1328 inverse[j][k] -= inverse[i][k]*t;
1333 return det*det >= tolerance*tolerance;
1336 template <
typename T>
1341 template <
typename T>
1350 #if DWREAL_IS_DOUBLE == 1
1354 #endif // DWREAL_IS_DOUBLE
1359 template<>
inline math::Mat4s zeroVal<math::Mat4s>() {
return math::Mat4s::identity(); }
1360 template<>
inline math::Mat4d zeroVal<math::Mat4d>() {
return math::Mat4d::identity(); }
1365 #endif // OPENVDB_UTIL_MAT4_H_HAS_BEEN_INCLUDED