MDL SDK API nvidia_logo_transpbg.gif Up
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
matrix.h
Go to the documentation of this file.
1 /***************************************************************************************************
2  * Copyright 2022 NVIDIA Corporation. All rights reserved.
3  **************************************************************************************************/
9 
10 #ifndef MI_MATH_MATRIX_H
11 #define MI_MATH_MATRIX_H
12 
13 #include <mi/base/types.h>
14 #include <mi/math/assert.h>
15 #include <mi/math/function.h>
16 #include <mi/math/vector.h>
17 
18 namespace mi {
19 
20 namespace math {
21 
48 //------- POD struct that provides storage for matrix elements --------
49 
88 template <typename T, Size ROW, Size COL>
90 {
91  T elements[ROW*COL];
92 };
93 
95 template <typename T> struct Matrix_struct<T,1,1>
96 {
97  T xx;
98 };
99 
101 template <typename T> struct Matrix_struct<T,2,1>
102 {
103  T xx;
104  T yx;
105 };
106 
108 template <typename T> struct Matrix_struct<T,3,1>
109 {
110  T xx;
111  T yx;
112  T zx;
113 };
114 
116 template <typename T> struct Matrix_struct<T,4,1>
117 {
118  T xx;
119  T yx;
120  T zx;
121  T wx;
122 };
123 
125 template <typename T> struct Matrix_struct<T,1,2>
126 {
127  T xx;
128  T xy;
129 };
130 
132 template <typename T> struct Matrix_struct<T,2,2>
133 {
134  T xx;
135  T xy;
136  T yx;
137  T yy;
138 };
139 
141 template <typename T> struct Matrix_struct<T,3,2>
142 {
143  T xx;
144  T xy;
145  T yx;
146  T yy;
147  T zx;
148  T zy;
149 };
150 
152 template <typename T> struct Matrix_struct<T,4,2>
153 {
154  T xx;
155  T xy;
156  T yx;
157  T yy;
158  T zx;
159  T zy;
160  T wx;
161  T wy;
162 };
163 
165 template <typename T> struct Matrix_struct<T,1,3>
166 {
167  T xx;
168  T xy;
169  T xz;
170 };
171 
173 template <typename T> struct Matrix_struct<T,2,3>
174 {
175  T xx;
176  T xy;
177  T xz;
178  T yx;
179  T yy;
180  T yz;
181 };
182 
184 template <typename T> struct Matrix_struct<T,3,3>
185 {
186  T xx;
187  T xy;
188  T xz;
189  T yx;
190  T yy;
191  T yz;
192  T zx;
193  T zy;
194  T zz;
195 };
196 
198 template <typename T> struct Matrix_struct<T,4,3>
199 {
200  T xx;
201  T xy;
202  T xz;
203  T yx;
204  T yy;
205  T yz;
206  T zx;
207  T zy;
208  T zz;
209  T wx;
210  T wy;
211  T wz;
212 };
213 
215 template <typename T> struct Matrix_struct<T,1,4>
216 {
217  T xx;
218  T xy;
219  T xz;
220  T xw;
221 };
222 
224 template <typename T> struct Matrix_struct<T,2,4>
225 {
226  T xx;
227  T xy;
228  T xz;
229  T xw;
230  T yx;
231  T yy;
232  T yz;
233  T yw;
234 };
235 
237 template <typename T> struct Matrix_struct<T,3,4>
238 {
239  T xx;
240  T xy;
241  T xz;
242  T xw;
243  T yx;
244  T yy;
245  T yz;
246  T yw;
247  T zx;
248  T zy;
249  T zz;
250  T zw;
251 };
252 
254 template <typename T> struct Matrix_struct<T,4,4>
255 {
256  T xx;
257  T xy;
258  T xz;
259  T xw;
260  T yx;
261  T yy;
262  T yz;
263  T yw;
264  T zx;
265  T zy;
266  T zz;
267  T zw;
268  T wx;
269  T wy;
270  T wz;
271  T ww;
272 };
273 
274 //------ Indirect access to matrix storage base ptr to keep Matrix_struct a POD --
275 
276 // Helper class to determine the matrix struct base pointer of the generic
277 // #Matrix_struct (\c specialized==\c false).
278 template <typename T, class Matrix, bool specialized>
279 struct Matrix_struct_get_base_pointer
280 {
281  static inline T* get_base_ptr( Matrix& m) { return m.elements; }
282  static inline const T* get_base_ptr( const Matrix& m) { return m.elements; }
283 };
284 
285 // Specialization of helper class to determine the matrix struct base pointer
286 // of a specialized #Matrix_struct (\c specialized==\c true).
287 template <typename T, class Matrix>
288 struct Matrix_struct_get_base_pointer<T,Matrix,true>
289 {
290  static inline T* get_base_ptr( Matrix& m) { return &m.xx; }
291  static inline const T* get_base_ptr( const Matrix& m) { return &m.xx; }
292 };
293 
294 
296 template <typename T, Size ROW, Size COL>
298 {
299  return Matrix_struct_get_base_pointer<T,Matrix_struct<T,ROW,COL>,
300  (ROW<=4 && COL<=4)>::get_base_ptr( mat);
301 }
302 
304 template <typename T, Size ROW, Size COL>
305 inline const T* matrix_base_ptr( const Matrix_struct<T,ROW,COL>& mat)
306 {
307  return Matrix_struct_get_base_pointer<T,Matrix_struct<T,ROW,COL>,
308  (ROW<=4 && COL<=4)>::get_base_ptr( mat);
309 }
310  // end group mi_math_matrix_struct
312 
313 
365 template <typename T, Size ROW, Size COL>
366 class Matrix : public Matrix_struct<T,ROW,COL> //-V690 PVS
367 {
368 public:
371 
372  typedef T value_type;
373  typedef Size size_type;
375  typedef T * pointer;
376  typedef const T * const_pointer;
377  typedef T & reference;
378  typedef const T & const_reference;
379 
382 
385 
386  static const Size ROWS = ROW;
387  static const Size COLUMNS = COL;
388  static const Size SIZE = ROW*COL;
389 
391  static inline Size size() { return SIZE; }
392 
394  static inline Size max_size() { return SIZE; }
395 
402  };
403 
405  inline T * begin() { return mi::math::matrix_base_ptr( *this); }
406 
408  inline T const * begin() const { return mi::math::matrix_base_ptr( *this); }
409 
413  inline T * end() { return begin() + SIZE; }
414 
418  inline T const * end() const { return begin() + SIZE; }
419 
422  {
423  mi_math_assert_msg( row < ROW, "precondition");
424  return reinterpret_cast<Row_vector&>(begin()[row * COL]);
425  }
426 
428  inline const Row_vector & operator[] (Size row) const
429  {
430  mi_math_assert_msg( row < ROW, "precondition");
431  return reinterpret_cast<const Row_vector&>(begin()[row * COL]);
432  }
433 
435  inline Matrix() { }
436 
437 #if (__cplusplus >= 201103L)
438  Matrix( const Matrix<T,ROW,COL>& other ) = default;
440 #endif
441 
443  inline Matrix( const Matrix_struct<T,ROW,COL>& other)
444  {
445  for( Size i(0u); i < SIZE; ++i)
446  begin()[i] = mi::math::matrix_base_ptr( other)[i];
447  }
448 
452  explicit inline Matrix( T diag)
453  {
454  for( Size i(0u); i < SIZE; ++i)
455  begin()[i] = T(0);
456  const Size MIN_DIM = (ROW < COL) ? ROW : COL;
457  for( Size k(0u); k < MIN_DIM; ++k)
458  begin()[k * COL + k] = diag;
459  }
460 
474  template <typename Iterator>
475  inline Matrix( From_iterator_tag, Iterator p)
476  {
477  for( Size i(0u); i < SIZE; ++i, ++p)
478  begin()[i] = *p;
479  }
480 
493  template <typename T2>
494  inline explicit Matrix( T2 const (& array)[SIZE])
495  {
496  for( Size i(0u); i < SIZE; ++i)
497  begin()[i] = array[i];
498  }
499 
502  template <typename T2>
503  inline explicit Matrix( const Matrix<T2,ROW,COL>& other)
504  {
505  for( Size i(0u); i < SIZE; ++i)
506  begin()[i] = T(other.begin()[i]);
507  }
508 
511  template <typename T2>
512  inline explicit Matrix( const Matrix_struct<T2,ROW,COL>& other)
513  {
514  for( Size i(0u); i < SIZE; ++i)
515  begin()[i] = T(mi::math::matrix_base_ptr( other)[i]);
516  }
517 
519  inline Matrix(
521  const Matrix<T,COL,ROW>& other)
522  {
523  for( Size i(0u); i < ROW; ++i)
524  for( Size j(0u); j < COL; ++j)
525  begin()[i * COL + j] = other.begin()[j * ROW + i];
526  }
527 
531  template <typename T2>
532  inline Matrix(
534  const Matrix<T2,COL,ROW>& other)
535  {
536  for( Size i(0u); i < ROW; ++i)
537  for( Size j(0u); j < COL; ++j)
538  begin()[i * COL + j] = T(other.begin()[j * ROW + i]);
539  }
540 
545  inline explicit Matrix(
546  const Row_vector& v0)
547  {
548  mi_static_assert( ROW == 1);
549  (*this)[0] = v0;
550  }
551 
556  inline Matrix(
557  const Row_vector& v0,
558  const Row_vector& v1)
559  {
560  mi_static_assert( ROW == 2);
561  (*this)[0] = v0;
562  (*this)[1] = v1;
563  }
564 
569  inline Matrix(
570  const Row_vector& v0,
571  const Row_vector& v1,
572  const Row_vector& v2)
573  {
574  mi_static_assert( ROW == 3);
575  (*this)[0] = v0;
576  (*this)[1] = v1;
577  (*this)[2] = v2;
578  }
579 
584  inline Matrix(
585  const Row_vector& v0,
586  const Row_vector& v1,
587  const Row_vector& v2,
588  const Row_vector& v3)
589  {
590  mi_static_assert( ROW == 4);
591  (*this)[0] = v0;
592  (*this)[1] = v1;
593  (*this)[2] = v2;
594  (*this)[3] = v3;
595  }
596 
600  inline Matrix( T m0, T m1)
601  {
602  mi_static_assert( SIZE == 2);
603  begin()[0] = m0; begin()[1] = m1;
604  }
605 
609  inline Matrix( T m0, T m1, T m2)
610  {
611  mi_static_assert( SIZE == 3);
612  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2;
613  }
614 
618  inline Matrix( T m0, T m1, T m2, T m3)
619  {
620  mi_static_assert( SIZE == 4);
621  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
622  }
623 
627  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5)
628  {
629  mi_static_assert( SIZE == 6);
630  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
631  begin()[4] = m4; begin()[5] = m5;
632  }
633 
637  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5, T m6, T m7)
638  {
639  mi_static_assert( SIZE == 8);
640  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
641  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
642  }
643 
647  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5, T m6, T m7, T m8)
648  {
649  mi_static_assert( SIZE == 9);
650  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
651  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
652  begin()[8] = m8;
653  }
654 
658  inline Matrix(
659  T m0, T m1, T m2, T m3,
660  T m4, T m5, T m6, T m7,
661  T m8, T m9, T m10, T m11)
662  {
663  mi_static_assert( SIZE == 12);
664  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
665  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
666  begin()[8] = m8; begin()[9] = m9; begin()[10] = m10; begin()[11] = m11;
667  }
668 
672  inline Matrix(
673  T m0, T m1, T m2, T m3,
674  T m4, T m5, T m6, T m7,
675  T m8, T m9, T m10, T m11,
676  T m12, T m13, T m14, T m15)
677  {
678  mi_static_assert( SIZE == 16);
679  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
680  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
681  begin()[8] = m8; begin()[9] = m9; begin()[10] = m10; begin()[11] = m11;
682  begin()[12] = m12; begin()[13] = m13; begin()[14] = m14; begin()[15] = m15;
683  }
684 
686  inline Matrix& operator=( const Matrix& other)
687  {
689  return *this;
690  }
691 
695  inline T& operator()( Size row, Size col)
696  {
697  mi_math_assert_msg( row < ROW, "precondition");
698  mi_math_assert_msg( col < COL, "precondition");
699  return begin()[row * COL + col];
700  }
701 
705  inline const T& operator()( Size row, Size col) const
706  {
707  mi_math_assert_msg( row < ROW, "precondition");
708  mi_math_assert_msg( col < COL, "precondition");
709  return begin()[row * COL + col];
710  }
711 
715  inline T get( Size i) const
716  {
717  mi_math_assert_msg( i < ROW*COL, "precondition");
718  return begin()[i];
719  }
720 
724  inline T get( Size row, Size col) const
725  {
726  mi_math_assert_msg( row < ROW, "precondition");
727  mi_math_assert_msg( col < COL, "precondition");
728  return begin()[row * COL + col];
729  }
730 
735  inline void set( Size i, T value)
736  {
737  mi_math_assert_msg( i < ROW*COL, "precondition");
738  begin()[i] = value;
739  }
740 
745  inline void set( Size row, Size col, T value)
746  {
747  mi_math_assert_msg( row < ROW, "precondition");
748  mi_math_assert_msg( col < COL, "precondition");
749  begin()[row * COL + col] = value;
750  }
751 
755  inline T det33() const
756  {
757  mi_static_assert( (ROW==3 || ROW==4) && (COL==3 || COL==4) );
758  return this->xx * this->yy * this->zz
759  + this->xy * this->yz * this->zx
760  + this->xz * this->yx * this->zy
761  - this->xx * this->zy * this->yz
762  - this->xy * this->zz * this->yx
763  - this->xz * this->zx * this->yy;
764  }
765 
769  inline bool invert();
770 
776  inline void transpose()
777  {
778  mi_static_assert( ROW==COL);
779  for( Size i=0; i < ROW-1; ++i) {
780  for( Size j=i+1; j < COL; ++j) {
781  T tmp = get(i,j);
782  set(i,j, get(j,i));
783  set(j,i, tmp);
784  }
785  }
786  }
787 
791  inline void translate( T x, T y, T z)
792  {
793  this->wx += x;
794  this->wy += y;
795  this->wz += z;
796  }
797 
801  inline void translate( const Vector<Float32,3>& vector)
802  {
803  this->wx += T( vector.x);
804  this->wy += T( vector.y);
805  this->wz += T( vector.z);
806  }
807 
811  inline void translate( const Vector<Float64,3>& vector)
812  {
813  this->wx += T( vector.x);
814  this->wy += T( vector.y);
815  this->wz += T( vector.z);
816  }
817 
821  inline void set_translation( T dx, T dy, T dz)
822  {
823  this->wx = dx;
824  this->wy = dy;
825  this->wz = dz;
826  }
827 
831  inline void set_translation( const Vector<Float32,3>& vector)
832  {
833  this->wx = T( vector.x);
834  this->wy = T( vector.y);
835  this->wz = T( vector.z);
836  }
837 
841  inline void set_translation( const Vector<Float64,3>& vector)
842  {
843  this->wx = T( vector.x);
844  this->wy = T( vector.y);
845  this->wz = T( vector.z);
846  }
847 
852  inline void rotate( T xangle, T yangle, T zangle)
853  {
854  Matrix<T,4,4> tmp( T( 1));
855  tmp.set_rotation( xangle, yangle, zangle);
856  (*this) *= tmp;
857  }
858 
863  inline void rotate( const Vector<Float32,3>& angles)
864  {
865  Matrix<T,4,4> tmp( T( 1));
866  tmp.set_rotation( T( angles.x), T( angles.y), T( angles.z));
867  (*this) *= tmp;
868  }
869 
874  inline void rotate( const Vector<Float64,3>& angles)
875  {
876  Matrix<T,4,4> tmp( T( 1));
877  tmp.set_rotation( T( angles.x), T( angles.y), T( angles.z));
878  (*this) *= tmp;
879  }
880 
886  inline void set_rotation( T x_angle, T y_angle, T z_angle);
887 
893  inline void set_rotation( const Vector<Float32,3>& angles)
894  {
895  set_rotation( T( angles.x), T( angles.y), T( angles.z));
896  }
897 
903  inline void set_rotation( const Vector<Float64,3>& angles)
904  {
905  set_rotation( T( angles.x), T( angles.y), T( angles.z));
906  }
907 
914  inline void set_rotation( const Vector<Float32,3>& axis, Float64 angle);
915 
922  inline void set_rotation( const Vector<Float64,3>& axis, Float64 angle);
923 
928  inline void lookat(
929  const Vector<Float32,3>& position,
930  const Vector<Float32,3>& target,
931  const Vector<Float32,3>& up);
932 
938  inline void lookat(
939  const Vector<Float64,3>& position,
940  const Vector<Float64,3>& target,
941  const Vector<Float64,3>& up);
942 };
943 
944 
945 //------ Free comparison operators ==, !=, <, <=, >, >= for matrices -----------
946 
948 template <typename T, Size ROW, Size COL>
949 inline bool operator==(
950  const Matrix<T,ROW,COL>& lhs,
951  const Matrix<T,ROW,COL>& rhs)
952 {
953  return is_equal( lhs, rhs);
954 }
955 
957 template <typename T, Size ROW, Size COL>
958 inline bool operator!=(
959  const Matrix<T,ROW,COL>& lhs,
960  const Matrix<T,ROW,COL>& rhs)
961 {
962  return is_not_equal( lhs, rhs);
963 }
964 
968 template <typename T, Size ROW, Size COL>
969 inline bool operator<(
970  const Matrix<T,ROW,COL>& lhs,
971  const Matrix<T,ROW,COL>& rhs)
972 {
973  return lexicographically_less( lhs, rhs);
974 }
975 
979 template <typename T, Size ROW, Size COL>
980 inline bool operator<=(
981  const Matrix<T,ROW,COL>& lhs,
982  const Matrix<T,ROW,COL>& rhs)
983 {
984  return lexicographically_less_or_equal( lhs, rhs);
985 }
986 
990 template <typename T, Size ROW, Size COL>
991 inline bool operator>(
992  const Matrix<T,ROW,COL>& lhs,
993  const Matrix<T,ROW,COL>& rhs)
994 {
995  return lexicographically_greater( lhs, rhs);
996 }
997 
1001 template <typename T, Size ROW, Size COL>
1002 inline bool operator>=(
1003  const Matrix<T,ROW,COL>& lhs,
1004  const Matrix<T,ROW,COL>& rhs)
1005 {
1006  return lexicographically_greater_or_equal( lhs, rhs);
1007 }
1008 
1009 //------ Operator declarations for Matrix -------------------------------------
1010 
1012 template <typename T, Size ROW, Size COL>
1013 Matrix<T,ROW,COL>& operator+=( Matrix<T,ROW,COL>& lhs, const Matrix<T,ROW,COL>& rhs);
1014 
1016 template <typename T, Size ROW, Size COL>
1017 Matrix<T,ROW,COL>& operator-=( Matrix<T,ROW,COL>& lhs, const Matrix<T,ROW,COL>& rhs);
1018 
1020 template <typename T, Size ROW, Size COL>
1022  const Matrix<T,ROW,COL>& lhs,
1023  const Matrix<T,ROW,COL>& rhs)
1024 {
1025  Matrix<T,ROW,COL> temp( lhs);
1026  temp += rhs;
1027  return temp;
1028 }
1029 
1031 template <typename T, Size ROW, Size COL>
1033  const Matrix<T,ROW,COL>& lhs,
1034  const Matrix<T,ROW,COL>& rhs)
1035 {
1036  Matrix<T,ROW,COL> temp( lhs);
1037  temp -= rhs;
1038  return temp;
1039 }
1040 
1042 template <typename T, Size ROW, Size COL>
1044  const Matrix<T,ROW,COL>& mat)
1045 {
1046  Matrix<T,ROW,COL> temp;
1047  for( Size i(0u); i < ROW*COL; ++i)
1048  temp.begin()[i] = -mat.get(i);
1049  return temp;
1050 }
1051 
1058 template <typename T, Size ROW, Size COL>
1060  Matrix<T,ROW,COL>& lhs,
1061  const Matrix<T,COL,COL>& rhs)
1062 {
1063  // There are more efficient ways of implementing this. Its a default solution. Specializations
1064  // exist below.
1065  Matrix<T,ROW,COL> old( lhs);
1066  for( Size rrow = 0; rrow < ROW; ++rrow) {
1067  for( Size rcol = 0; rcol < COL; ++rcol) {
1068  lhs( rrow, rcol) = T(0);
1069  for( Size k = 0; k < COL; ++k)
1070  lhs( rrow, rcol) += old( rrow, k) * rhs( k, rcol);
1071  }
1072  }
1073  return lhs;
1074 }
1075 
1081 template <typename T, Size ROW1, Size COL1, Size ROW2, Size COL2>
1083  const Matrix<T,ROW1,COL1>& lhs,
1084  const Matrix<T,ROW2,COL2>& rhs)
1085 {
1086  // There are more efficient ways of implementing this. Its a default solution. Specializations
1087  // exist below.
1088  mi_static_assert( COL1 == ROW2);
1089  Matrix<T,ROW1,COL2> result;
1090  for( Size rrow = 0; rrow < ROW1; ++rrow) {
1091  for( Size rcol = 0; rcol < COL2; ++rcol) {
1092  result( rrow, rcol) = T(0);
1093  for( Size k = 0; k < COL1; ++k)
1094  result( rrow, rcol) += lhs( rrow, k) * rhs( k, rcol);
1095  }
1096  }
1097  return result;
1098 }
1099 
1107 template <typename T, Size ROW, Size COL, Size DIM>
1109  const Matrix<T,ROW,COL>& mat,
1110  const Vector<T,DIM>& vec)
1111 {
1112  mi_static_assert( COL == DIM);
1113  Vector<T,ROW> result;
1114  for( Size row = 0; row < ROW; ++row) {
1115  result[row] = T(0);
1116  for( Size col = 0; col < COL; ++col)
1117  result[row] += mat( row, col) * vec[col];
1118  }
1119  return result;
1120 }
1121 
1128 template <Size DIM, typename T, Size ROW, Size COL>
1130  const Vector<T,DIM>& vec,
1131  const Matrix<T,ROW,COL>& mat)
1132 {
1133  mi_static_assert( DIM == ROW);
1134  Vector<T,COL> result;
1135  for( Size col = 0; col < COL; ++col) {
1136  result[col] = T(0);
1137  for( Size row = 0; row < ROW; ++row)
1138  result[col] += mat( row, col) * vec[row];
1139  }
1140  return result;
1141 }
1142 
1145 template <typename T, Size ROW, Size COL>
1147 {
1148  for( Size i=0; i < ROW*COL; ++i)
1149  mat.begin()[i] *= factor;
1150  return mat;
1151 }
1152 
1154 template <typename T, Size ROW, Size COL>
1155 inline Matrix<T,ROW,COL> operator*( const Matrix<T,ROW,COL>& mat, T factor)
1156 {
1157  Matrix<T,ROW,COL> temp( mat);
1158  temp *= factor;
1159  return temp;
1160 }
1161 
1163 template <typename T, Size ROW, Size COL>
1164 inline Matrix<T,ROW,COL> operator*( T factor, const Matrix<T,ROW,COL>& mat)
1165 {
1166  Matrix<T,ROW,COL> temp( mat);
1167  temp *= factor; // * on T should be commutative (IEEE-754)
1168  return temp;
1169 }
1170 
1171 
1172 //------ Free Functions for Matrix --------------------------------------------
1173 
1177 template <Size NEW_ROW, Size NEW_COL, typename T, Size ROW, Size COL>
1179  const Matrix<T,ROW,COL>& mat)
1180 {
1181  mi_static_assert( NEW_ROW <= ROW);
1182  mi_static_assert( NEW_COL <= COL);
1184  for( Size i=0; i < NEW_ROW; ++i)
1185  for( Size j=0; j < NEW_COL; ++j)
1186  result( i, j) = mat( i, j);
1187  return result;
1188 }
1189 
1190 
1191 
1193 template <typename T, Size ROW, Size COL>
1195  const Matrix<T,ROW,COL>& mat)
1196 {
1197  Matrix<T,COL,ROW> result;
1198  for( Size i=0; i < ROW; ++i)
1199  for( Size j=0; j < COL; ++j)
1200  result( j, i) = mat( i, j);
1201  return result;
1202 }
1203 
1209 template <typename T, typename U>
1211  const Matrix<T,4,4>& mat,
1212  const U& point)
1213 {
1214  const T w = T(mat.xw * point + mat.ww);
1215  if( w == T(0) || w == T(1))
1216  return U(mat.xx * point + mat.wx);
1217  else
1218  return U((mat.xx * point + mat.wx) / w);
1219 }
1220 
1226 template <typename T, typename U>
1228  const Matrix<T,4,4>& mat,
1229  const Vector<U,2>& point)
1230 {
1231  T w = T(mat.xw * point.x + mat.yw * point.y + mat.ww);
1232  if( w == T(0) || w == T(1))
1233  return Vector<U, 2>(
1234  U(mat.xx * point.x + mat.yx * point.y + mat.wx),
1235  U(mat.xy * point.x + mat.yy * point.y + mat.wy));
1236  else {
1237  w = T(1)/w; // optimization
1238  return Vector<U, 2>(
1239  U((mat.xx * point.x + mat.yx * point.y + mat.wx) * w),
1240  U((mat.xy * point.x + mat.yy * point.y + mat.wy) * w));
1241  }
1242 }
1243 
1249 template <typename T, typename U>
1251  const Matrix<T,4,3>& mat,
1252  const Vector<U,3>& point)
1253 {
1254  return Vector<U,3>(
1255  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx),
1256  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy),
1257  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz));
1258 }
1259 
1265 template <typename T, typename U>
1267  const Matrix<T,4,4>& mat,
1268  const Vector<U,3>& point)
1269 {
1270  T w = T(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww);
1271  if( w == T(0) || w == T(1)) // avoids temporary and division
1272  return Vector<U,3>(
1273  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx),
1274  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy),
1275  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz));
1276  else {
1277  w = T(1)/w; // optimization
1278  return Vector<U,3>(
1279  U((mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx) * w),
1280  U((mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy) * w),
1281  U((mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz) * w));
1282  }
1283 }
1284 
1290 template <typename T, typename U>
1292  const Matrix<T,4,4>& mat,
1293  const Vector<U,4>& point)
1294 {
1295  return Vector<U, 4>(
1296  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx * point.w),
1297  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy * point.w),
1298  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz * point.w),
1299  U(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww * point.w));
1300 }
1301 
1307 template <typename T, typename U>
1309  const Matrix<T,4,4>& mat,
1310  const U& vector)
1311 {
1312  return U(mat.xx * vector);
1313 }
1314 
1320 template <typename T, typename U>
1322  const Matrix<T,4,4>& mat,
1323  const Vector<U,2>& vector)
1324 {
1325  return Vector<U,2>(
1326  U(mat.xx * vector.x + mat.yx * vector.y),
1327  U(mat.xy * vector.x + mat.yy * vector.y));
1328 }
1329 
1335 template <typename T, typename U>
1337  const Matrix<T,3,3>& mat,
1338  const Vector<U,3>& vector)
1339 {
1340  return Vector<U,3>(
1341  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1342  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1343  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1344 }
1345 
1351 template <typename T, typename U>
1353  const Matrix<T,4,3>& mat,
1354  const Vector<U,3>& vector)
1355 {
1356  return Vector<U,3>(
1357  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1358  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1359  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1360 }
1361 
1367 template <typename T, typename U>
1369  const Matrix<T,4,4>& mat,
1370  const Vector<U,3>& vector)
1371 {
1372  return Vector<U,3>(
1373  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1374  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1375  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1376 }
1377 
1389 template <typename T, typename U>
1391  const Matrix<T,3,3>& inv_mat,
1392  const Vector<U,3>& normal)
1393 {
1394  return Vector<U,3>(
1395  U(inv_mat.xx * normal.x + inv_mat.xy * normal.y + inv_mat.xz * normal.z),
1396  U(inv_mat.yx * normal.x + inv_mat.yy * normal.y + inv_mat.yz * normal.z),
1397  U(inv_mat.zx * normal.x + inv_mat.zy * normal.y + inv_mat.zz * normal.z));
1398 }
1399 
1411 template <typename T, typename U>
1413  const Matrix<T,4,4>& inv_mat,
1414  const Vector<U,3>& normal)
1415 {
1416  return Vector<U,3>(
1417  U(inv_mat.xx * normal.x + inv_mat.xy * normal.y + inv_mat.xz * normal.z),
1418  U(inv_mat.yx * normal.x + inv_mat.yy * normal.y + inv_mat.yz * normal.z),
1419  U(inv_mat.zx * normal.x + inv_mat.zy * normal.y + inv_mat.zz * normal.z));
1420 }
1421 
1435 template <typename T, typename U>
1437  const Matrix<T,4,4>& mat,
1438  const Vector<U,3>& normal)
1439 {
1440  Matrix<T,3,3> sub_mat( mat.xx, mat.xy, mat.xz,
1441  mat.yx, mat.yy, mat.yz,
1442  mat.zx, mat.zy, mat.zz);
1443  bool inverted = sub_mat.invert();
1444  if( !inverted)
1445  return normal;
1446  return Vector<U,3>(
1447  U(sub_mat.xx * normal.x + sub_mat.xy * normal.y + sub_mat.xz * normal.z),
1448  U(sub_mat.yx * normal.x + sub_mat.yy * normal.y + sub_mat.yz * normal.z),
1449  U(sub_mat.zx * normal.x + sub_mat.zy * normal.y + sub_mat.zz * normal.z));
1450 }
1451 
1452 
1453 //------ Definitions of non-inline function -----------------------------------
1454 
1455 template <typename T, Size ROW, Size COL>
1457  Matrix<T,ROW,COL>& lhs,
1458  const Matrix<T,ROW,COL>& rhs)
1459 {
1460  for( Size i=0; i < ROW*COL; ++i)
1461  lhs.begin()[i] += rhs.get(i);
1462  return lhs;
1463 }
1464 
1465 template <typename T, Size ROW, Size COL>
1467  Matrix<T,ROW,COL>& lhs,
1468  const Matrix<T,ROW,COL>& rhs)
1469 {
1470  for( Size i=0; i < ROW*COL; ++i)
1471  lhs.begin()[i] -= rhs.get(i);
1472  return lhs;
1473 }
1474 
1475 #ifndef MI_FOR_DOXYGEN_ONLY
1476 
1477 template <typename T, Size ROW, Size COL>
1478 inline void Matrix<T,ROW,COL>::set_rotation( T xangle, T yangle, T zangle)
1479 {
1480  mi_static_assert( COL == 4 && ROW == 4);
1481  T tsx, tsy, tsz; // sine of [xyz]_angle
1482  T tcx, tcy, tcz; // cosine of [xyz]_angle
1483  T tmp;
1484  const T min_angle = T(0.00024f);
1485 
1486  if( abs( xangle) > min_angle) {
1487  tsx = sin( xangle);
1488  tcx = cos( xangle);
1489  } else {
1490  tsx = xangle;
1491  tcx = T(1);
1492  }
1493  if( abs( yangle) > min_angle) {
1494  tsy = sin( yangle);
1495  tcy = cos( yangle);
1496  } else {
1497  tsy = yangle;
1498  tcy = T(1);
1499  }
1500  if( abs(zangle) > min_angle) {
1501  tsz = sin( zangle);
1502  tcz = cos( zangle);
1503  } else {
1504  tsz = T(zangle);
1505  tcz = T(1);
1506  }
1507  this->xx = tcy * tcz;
1508  this->xy = tcy * tsz;
1509  this->xz = -tsy;
1510 
1511  tmp = tsx * tsy;
1512  this->yx = tmp * tcz - tcx * tsz;
1513  this->yy = tmp * tsz + tcx * tcz;
1514  this->yz = tsx * tcy;
1515 
1516  tmp = tcx * tsy;
1517  this->zx = tmp * tcz + tsx * tsz;
1518  this->zy = tmp * tsz - tsx * tcz;
1519  this->zz = tcx * tcy;
1520 }
1521 
1522 template <typename T, Size ROW, Size COL>
1523 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float32,3>& axis_v, Float64 angle)
1524 {
1525  mi_static_assert( COL == 4 && ROW == 4);
1526  Vector<T,3> axis( axis_v);
1527  const T min_angle = T(0.00024f);
1528 
1529  if( abs( T(angle)) < min_angle) {
1530  T xa = axis.x * T(angle);
1531  T ya = axis.y * T(angle);
1532  T za = axis.z * T(angle);
1533 
1534  this->xx = T(1);
1535  this->xy = za;
1536  this->xz = -ya;
1537  this->xw = T(0);
1538 
1539  this->yx = za;
1540  this->yy = T(1);
1541  this->yz = xa;
1542  this->yw = T(0);
1543 
1544  this->zx = -ya;
1545  this->zy = -xa;
1546  this->zz = T(1);
1547  this->zw = T(0);
1548  } else {
1549  T s = sin( T(angle));
1550  T c = cos( T(angle));
1551  T t = T(1) - c;
1552  T tmp;
1553 
1554  tmp = t * T(axis.x);
1555  this->xx = tmp * T(axis.x) + c;
1556  this->xy = tmp * T(axis.y) + s * T(axis.z);
1557  this->xz = tmp * T(axis.z) - s * T(axis.y);
1558  this->xw = T(0);
1559 
1560  tmp = t * T(axis.y);
1561  this->yx = tmp * T(axis.x) - s * T(axis.z);
1562  this->yy = tmp * T(axis.y) + c;
1563  this->yz = tmp * T(axis.z) + s * T(axis.x);
1564  this->yw = T(0);
1565 
1566  tmp = t * T(axis.z);
1567  this->zx = tmp * T(axis.x) + s * T(axis.y);
1568  this->zy = tmp * T(axis.y) - s * T(axis.x);
1569  this->zz = tmp * T(axis.z) + c;
1570  this->zw = T(0);
1571  }
1572  this->wx = this->wy = this->wz = T(0);
1573  this->ww = T(1);
1574 }
1575 
1576 template <typename T, Size ROW, Size COL>
1577 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float64,3>& axis_v, Float64 angle)
1578 {
1579  mi_static_assert( COL == 4 && ROW == 4);
1580  Vector<T,3> axis( axis_v);
1581  const T min_angle = T(0.00024f);
1582 
1583  if( abs(T(angle)) < min_angle) {
1584  T xa = axis.x * T(angle);
1585  T ya = axis.y * T(angle);
1586  T za = axis.z * T(angle);
1587 
1588  this->xx = T(1);
1589  this->xy = za;
1590  this->xz = -ya;
1591  this->xw = T(0);
1592 
1593  this->yx = za;
1594  this->yy = T(1);
1595  this->yz = xa;
1596  this->yw = T(0);
1597 
1598  this->zx = -ya;
1599  this->zy = -xa;
1600  this->zz = T(1);
1601  this->zw = T(0);
1602  } else {
1603  T s = sin( T(angle));
1604  T c = cos( T(angle));
1605  T t = T(1) - c;
1606  T tmp;
1607 
1608  tmp = t * T(axis.x);
1609  this->xx = tmp * T(axis.x) + c;
1610  this->xy = tmp * T(axis.y) + s * T(axis.z);
1611  this->xz = tmp * T(axis.z) - s * T(axis.y);
1612  this->xw = T(0);
1613 
1614  tmp = t * T(axis.y);
1615  this->yx = tmp * T(axis.x) - s * T(axis.z);
1616  this->yy = tmp * T(axis.y) + c;
1617  this->yz = tmp * T(axis.z) + s * T(axis.x);
1618  this->yw = T(0);
1619 
1620  tmp = t * T(axis.z);
1621  this->zx = tmp * T(axis.x) + s * T(axis.y);
1622  this->zy = tmp * T(axis.y) - s * T(axis.x);
1623  this->zz = tmp * T(axis.z) + c;
1624  this->zw = T(0);
1625  }
1626  this->wx = this->wy = this->wz = T(0);
1627  this->ww = T(1);
1628 }
1629 
1630 template <typename T, Size ROW, Size COL>
1631 inline void Matrix<T,ROW,COL>::lookat(
1632  const Vector<Float32,3>& position,
1633  const Vector<Float32,3>& target,
1634  const Vector<Float32,3>& up)
1635 {
1636  mi_static_assert( COL == 4 && ROW == 4);
1637  Vector<Float32,3> xaxis, yaxis, zaxis;
1638 
1639  // Z vector
1640  zaxis = position - target;
1641  zaxis.normalize();
1642 
1643  // X vector = up cross Z
1644  xaxis = cross( up, zaxis);
1645  xaxis.normalize();
1646 
1647  // Recompute Y = Z cross X
1648  yaxis = cross( zaxis, xaxis);
1649  yaxis.normalize();
1650 
1651  // Build rotation matrix
1652  Matrix<T,4,4> rot(
1653  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1654  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1655  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1656  T(0), T(0), T(0), T(1));
1657 
1658  // Compute the new position by multiplying the inverse position with the rotation matrix
1659  Matrix<T,4,4> trans(
1660  T(1), T(0), T(0), T(0),
1661  T(0), T(1), T(0), T(0),
1662  T(0), T(0), T(1), T(0),
1663  T(-position.x), T(-position.y), T(-position.z), T(1));
1664 
1665  *this = trans * rot;
1666 }
1667 
1668 template <typename T, Size ROW, Size COL>
1669 inline void Matrix<T,ROW,COL>::lookat(
1670  const Vector<Float64,3>& position,
1671  const Vector<Float64,3>& target,
1672  const Vector<Float64,3>& up)
1673 {
1674  mi_static_assert( COL == 4 && ROW == 4);
1675  Vector<Float64,3> xaxis, yaxis, zaxis;
1676 
1677  // Z vector
1678  zaxis = position - target;
1679  zaxis.normalize();
1680 
1681  // X vector = up cross Z
1682  xaxis = cross( up, zaxis);
1683  xaxis.normalize();
1684 
1685  // Recompute Y = Z cross X
1686  yaxis = cross( zaxis, xaxis);
1687  yaxis.normalize();
1688 
1689  // Build rotation matrix
1690  Matrix<T,4,4> rot(
1691  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1692  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1693  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1694  T(0), T(0), T(0), T(1));
1695 
1696  // Compute the new position by multiplying the inverse position with the rotation matrix
1697  Matrix<T,4,4> trans(
1698  T(1), T(0), T(0), T(0),
1699  T(0), T(1), T(0), T(0),
1700  T(0), T(0), T(1), T(0),
1701  T(-position.x), T(-position.y), T(-position.z), T(1));
1702 
1703  *this = trans * rot;
1704 }
1705 
1706 
1707 //------ Definition and helper for matrix inversion ---------------------------
1708 
1709 // Matrix inversion class, specialized for symmetric matrices and low dimensions
1710 template <class T, Size ROW, Size COL>
1711 class Matrix_inverter
1712 {
1713 public:
1714  typedef math::Matrix<T,ROW,COL> Matrix;
1715 
1716  // Inverts the matrix \c M if possible and returns \c true.
1717  //
1718  // If the matrix cannot be inverted or if \c ROW != \c COL, this function returns \c false.
1719  static inline bool invert( Matrix& /* M */) { return false; }
1720 };
1721 
1722 // Matrix inversion class, specialization for symmetric matrices
1723 template <class T, Size DIM>
1724 class Matrix_inverter<T,DIM,DIM>
1725 {
1726 public:
1727  typedef math::Matrix<T,DIM,DIM> Matrix;
1728  typedef math::Vector<T,DIM> Value_vector;
1729  typedef math::Vector<Size,DIM> Index_vector;
1730 
1731  // LU decomposition of matrix lu in place.
1732  //
1733  // Returns \c false if matrix cannot be decomposed, for example, if it is not invertible, and \c
1734  // true otherwise. Returns also a row permutation in indx.
1735  static bool lu_decomposition(
1736  Matrix& lu, // matrix to decompose, modified in place
1737  Index_vector& indx); // row permutation info indx[4]
1738 
1739  // Solves the equation lu * x = b for x. lu and indx are the results of lu_decomposition. The
1740  // solution is returned in b.
1741  static void lu_backsubstitution(
1742  const Matrix& lu, // LU decomposed matrix
1743  const Index_vector& indx, // permutation vector
1744  Value_vector& b); // right hand side vector b, modified in place
1745 
1746  static bool invert( Matrix& mat);
1747 };
1748 
1749 template <class T, Size DIM>
1750 bool Matrix_inverter<T,DIM,DIM>::lu_decomposition(
1751  Matrix& lu,
1752  Index_vector& indx)
1753 {
1754  Value_vector vv; // implicit scaling of each row
1755 
1756  for( Size i = 0; i < DIM; i++) { // get implicit scaling
1757  T big = T(0);
1758  for( Size j = 0; j < DIM; j++) {
1759  T temp = abs(lu.get(i,j));
1760  if( temp > big)
1761  big = temp;
1762  }
1763  if( big == T(0))
1764  return false;
1765  vv[i] = T(1) / big; // save scaling
1766  }
1767  //
1768  // loop over columns of Crout's method
1769  //
1770  Size imax = 0;
1771  for( Size j = 0; j < DIM; j++) {
1772  for( Size i = 0; i < j; i++) {
1773  T sum = lu.get(i,j);
1774  for( Size k = 0; k < i; k++)
1775  sum -= lu.get(i,k) * lu.get(k,j);
1776  lu.set(i, j, sum);
1777  }
1778  T big = 0; // init for pivot search
1779  for( Size i = j; i < DIM; i++) {
1780  T sum = lu.get(i,j);
1781  for( Size k = 0; k < j; k++)
1782  sum -= lu.get(i,k) * lu.get(k,j);
1783  lu.set(i, j, sum);
1784  T dum = vv[i] * abs(sum);
1785  if( dum >= big) { // pivot good?
1786  big = dum;
1787  imax = i;
1788  }
1789  }
1790  if( j != imax) { // interchange rows
1791  for( Size k = 0; k < DIM; k++) {
1792  T dum = lu.get(imax,k);
1793  lu.set(imax, k, lu.get(j,k));
1794  lu.set(j, k, dum);
1795  }
1796  vv[imax] = vv[j]; // interchange scale factor
1797  }
1798  indx[j] = imax;
1799  if( lu.get(j,j) == 0)
1800  return false;
1801  if( j != DIM-1) { // divide by pivot element
1802  T dum = T(1) / lu.get(j,j);
1803  for( Size i = j + 1; i < DIM; i++)
1804  lu.set(i, j, lu.get(i,j) * dum);
1805  }
1806  }
1807  return true;
1808 }
1809 
1810 template <class T, Size DIM>
1811 void Matrix_inverter<T,DIM,DIM>::lu_backsubstitution(
1812  const Matrix& lu,
1813  const Index_vector& indx,
1814  Value_vector& b)
1815 {
1816  // when ii != DIM+1, ii is index of first non-vanishing element of b
1817  Size ii = DIM+1;
1818  for( Size i = 0; i < DIM; i++) {
1819  Size ip = indx[i];
1820  T sum = b[ip];
1821  b[ip] = b[i];
1822  if( ii != DIM+1) {
1823  for( Size j = ii; j < i; j++) {
1824  sum -= lu.get(i,j) * b[j];
1825  }
1826  } else {
1827  if( sum != T(0)) { // non-zero element
1828  ii = i;
1829  }
1830  }
1831  b[i] = sum;
1832  }
1833  for( Size i2 = DIM; i2 > 0;) { // backsubstitution
1834  --i2;
1835  T sum = b[i2];
1836  for( Size j = i2+1; j < DIM; j++)
1837  sum -= lu.get(i2,j) * b[j];
1838  b[i2] = sum / lu.get(i2,i2); // store comp. of sol. vector
1839  }
1840 }
1841 
1842 template <class T, Size DIM>
1843 bool Matrix_inverter<T,DIM,DIM>::invert( Matrix& mat)
1844 {
1845  Matrix lu(mat); // working copy of matrix to invert
1846  Index_vector indx; // row permutation info
1847 
1848  // LU decompose, return if it fails
1849  if( !lu_decomposition(lu, indx))
1850  return false;
1851 
1852  // solve for each column vector of the I matrix
1853  for( Size j = 0; j < DIM; ++j) {
1854  Value_vector col(T(0)); // TODO: do that directly in the mat matrix
1855  col[j] = T(1);
1856  lu_backsubstitution( lu, indx, col);
1857  for( Size i = 0; i < DIM; ++i) {
1858  mat.set( i, j, col[i]);
1859  }
1860  }
1861  return true;
1862 }
1863 
1864 // Matrix inversion class, specialization for 1x1 matrix
1865 template <class T>
1866 class Matrix_inverter<T,1,1>
1867 {
1868 public:
1869  typedef math::Matrix<T,1,1> Matrix;
1870 
1871  static inline bool invert( Matrix& mat)
1872  {
1873  T s = mat.get( 0, 0);
1874  if( s == T(0))
1875  return false;
1876  mat.set( 0, 0, T(1) / s);
1877  return true;
1878  }
1879 };
1880 
1881 // Matrix inversion class, specialization for 2x2 matrix
1882 template <class T>
1883 class Matrix_inverter<T,2,2>
1884 {
1885 public:
1886  typedef math::Matrix<T,2,2> Matrix;
1887 
1888  static inline bool invert( Matrix& mat)
1889  {
1890  T a = mat.get( 0, 0);
1891  T b = mat.get( 0, 1);
1892  T c = mat.get( 1, 0);
1893  T d = mat.get( 1, 1);
1894  T det = a*d - b*c;
1895  if( det == T( 0))
1896  return false;
1897  T rdet = T(1) / det;
1898  mat.set( 0, 0, d * rdet);
1899  mat.set( 0, 1,-b * rdet);
1900  mat.set( 1, 0,-c * rdet);
1901  mat.set( 1, 1, a * rdet);
1902  return true;
1903  }
1904 };
1905 
1906 template <typename T, Size ROW, Size COL>
1907 inline bool Matrix<T,ROW,COL>::invert()
1908 {
1909  return Matrix_inverter<T,ROW,COL>::invert( *this);
1910 }
1911 
1912 
1913 
1914 //------ Specializations and (specialized) overloads --------------------------
1915 
1916 // Specialization of common matrix multiplication for 4x4 matrices.
1917 template <typename T>
1918 inline Matrix<T,4,4>& operator*=(
1919  Matrix<T,4,4>& lhs,
1920  const Matrix<T,4,4>& rhs)
1921 {
1922  Matrix<T,4,4> old( lhs);
1923 
1924  lhs.xx = old.xx * rhs.xx + old.xy * rhs.yx + old.xz * rhs.zx + old.xw * rhs.wx;
1925  lhs.xy = old.xx * rhs.xy + old.xy * rhs.yy + old.xz * rhs.zy + old.xw * rhs.wy;
1926  lhs.xz = old.xx * rhs.xz + old.xy * rhs.yz + old.xz * rhs.zz + old.xw * rhs.wz;
1927  lhs.xw = old.xx * rhs.xw + old.xy * rhs.yw + old.xz * rhs.zw + old.xw * rhs.ww;
1928 
1929  lhs.yx = old.yx * rhs.xx + old.yy * rhs.yx + old.yz * rhs.zx + old.yw * rhs.wx;
1930  lhs.yy = old.yx * rhs.xy + old.yy * rhs.yy + old.yz * rhs.zy + old.yw * rhs.wy;
1931  lhs.yz = old.yx * rhs.xz + old.yy * rhs.yz + old.yz * rhs.zz + old.yw * rhs.wz;
1932  lhs.yw = old.yx * rhs.xw + old.yy * rhs.yw + old.yz * rhs.zw + old.yw * rhs.ww;
1933 
1934  lhs.zx = old.zx * rhs.xx + old.zy * rhs.yx + old.zz * rhs.zx + old.zw * rhs.wx;
1935  lhs.zy = old.zx * rhs.xy + old.zy * rhs.yy + old.zz * rhs.zy + old.zw * rhs.wy;
1936  lhs.zz = old.zx * rhs.xz + old.zy * rhs.yz + old.zz * rhs.zz + old.zw * rhs.wz;
1937  lhs.zw = old.zx * rhs.xw + old.zy * rhs.yw + old.zz * rhs.zw + old.zw * rhs.ww;
1938 
1939  lhs.wx = old.wx * rhs.xx + old.wy * rhs.yx + old.wz * rhs.zx + old.ww * rhs.wx;
1940  lhs.wy = old.wx * rhs.xy + old.wy * rhs.yy + old.wz * rhs.zy + old.ww * rhs.wy;
1941  lhs.wz = old.wx * rhs.xz + old.wy * rhs.yz + old.wz * rhs.zz + old.ww * rhs.wz;
1942  lhs.ww = old.wx * rhs.xw + old.wy * rhs.yw + old.wz * rhs.zw + old.ww * rhs.ww;
1943 
1944  return lhs;
1945 }
1946 
1947 // Specialization of common matrix multiplication for 4x4 matrices.
1948 template <typename T>
1949 inline Matrix<T,4,4> operator*(
1950  const Matrix<T,4,4>& lhs,
1951  const Matrix<T,4,4>& rhs)
1952 {
1953  Matrix<T,4,4> temp( lhs);
1954  temp *= rhs;
1955  return temp;
1956 }
1957 
1958 #endif // MI_FOR_DOXYGEN_ONLY
1959  // end group mi_math_matrix
1961 
1962 } // namespace math
1963 
1964 } // namespace mi
1965 
1966 #endif // MI_MATH_MATRIX_H