Material Definition Language 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 2020 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>
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 
438  inline Matrix( const Matrix_struct<T,ROW,COL>& other)
439  {
440  for( Size i(0u); i < SIZE; ++i)
441  begin()[i] = mi::math::matrix_base_ptr( other)[i];
442  }
443 
447  explicit inline Matrix( T diag)
448  {
449  for( Size i(0u); i < SIZE; ++i)
450  begin()[i] = T(0);
451  const Size MIN_DIM = (ROW < COL) ? ROW : COL;
452  for( Size k(0u); k < MIN_DIM; ++k)
453  begin()[k * COL + k] = diag;
454  }
455 
469  template <typename Iterator>
470  inline Matrix( From_iterator_tag, Iterator p)
471  {
472  for( Size i(0u); i < SIZE; ++i, ++p)
473  begin()[i] = *p;
474  }
475 
488  template <typename T2>
489  inline explicit Matrix( T2 const (& array)[SIZE])
490  {
491  for( Size i(0u); i < SIZE; ++i)
492  begin()[i] = array[i];
493  }
494 
497  template <typename T2>
498  inline explicit Matrix( const Matrix<T2,ROW,COL>& other)
499  {
500  for( Size i(0u); i < SIZE; ++i)
501  begin()[i] = T(other.begin()[i]);
502  }
503 
506  template <typename T2>
507  inline explicit Matrix( const Matrix_struct<T2,ROW,COL>& other)
508  {
509  for( Size i(0u); i < SIZE; ++i)
510  begin()[i] = T(mi::math::matrix_base_ptr( other)[i]);
511  }
512 
514  inline Matrix(
516  const Matrix<T,COL,ROW>& other)
517  {
518  for( Size i(0u); i < ROW; ++i)
519  for( Size j(0u); j < COL; ++j)
520  begin()[i * COL + j] = other.begin()[j * ROW + i];
521  }
522 
526  template <typename T2>
527  inline Matrix(
529  const Matrix<T2,COL,ROW>& other)
530  {
531  for( Size i(0u); i < ROW; ++i)
532  for( Size j(0u); j < COL; ++j)
533  begin()[i * COL + j] = T(other.begin()[j * ROW + i]);
534  }
535 
540  inline explicit Matrix(
541  const Row_vector& v0)
542  {
543  mi_static_assert( ROW == 1);
544  (*this)[0] = v0;
545  }
546 
551  inline Matrix(
552  const Row_vector& v0,
553  const Row_vector& v1)
554  {
555  mi_static_assert( ROW == 2);
556  (*this)[0] = v0;
557  (*this)[1] = v1;
558  }
559 
564  inline Matrix(
565  const Row_vector& v0,
566  const Row_vector& v1,
567  const Row_vector& v2)
568  {
569  mi_static_assert( ROW == 3);
570  (*this)[0] = v0;
571  (*this)[1] = v1;
572  (*this)[2] = v2;
573  }
574 
579  inline Matrix(
580  const Row_vector& v0,
581  const Row_vector& v1,
582  const Row_vector& v2,
583  const Row_vector& v3)
584  {
585  mi_static_assert( ROW == 4);
586  (*this)[0] = v0;
587  (*this)[1] = v1;
588  (*this)[2] = v2;
589  (*this)[3] = v3;
590  }
591 
595  inline Matrix( T m0, T m1)
596  {
597  mi_static_assert( SIZE == 2);
598  begin()[0] = m0; begin()[1] = m1;
599  }
600 
604  inline Matrix( T m0, T m1, T m2)
605  {
606  mi_static_assert( SIZE == 3);
607  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2;
608  }
609 
613  inline Matrix( T m0, T m1, T m2, T m3)
614  {
615  mi_static_assert( SIZE == 4);
616  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
617  }
618 
622  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5)
623  {
624  mi_static_assert( SIZE == 6);
625  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
626  begin()[4] = m4; begin()[5] = m5;
627  }
628 
632  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5, T m6, T m7)
633  {
634  mi_static_assert( SIZE == 8);
635  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
636  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
637  }
638 
642  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5, T m6, T m7, T m8)
643  {
644  mi_static_assert( SIZE == 9);
645  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
646  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
647  begin()[8] = m8;
648  }
649 
653  inline Matrix(
654  T m0, T m1, T m2, T m3,
655  T m4, T m5, T m6, T m7,
656  T m8, T m9, T m10, T m11)
657  {
658  mi_static_assert( SIZE == 12);
659  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
660  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
661  begin()[8] = m8; begin()[9] = m9; begin()[10] = m10; begin()[11] = m11;
662  }
663 
667  inline Matrix(
668  T m0, T m1, T m2, T m3,
669  T m4, T m5, T m6, T m7,
670  T m8, T m9, T m10, T m11,
671  T m12, T m13, T m14, T m15)
672  {
673  mi_static_assert( SIZE == 16);
674  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
675  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
676  begin()[8] = m8; begin()[9] = m9; begin()[10] = m10; begin()[11] = m11;
677  begin()[12] = m12; begin()[13] = m13; begin()[14] = m14; begin()[15] = m15;
678  }
679 
681  inline Matrix& operator=( const Matrix& other)
682  {
684  return *this;
685  }
686 
690  inline T& operator()( Size row, Size col)
691  {
692  mi_math_assert_msg( row < ROW, "precondition");
693  mi_math_assert_msg( col < COL, "precondition");
694  return begin()[row * COL + col];
695  }
696 
700  inline const T& operator()( Size row, Size col) const
701  {
702  mi_math_assert_msg( row < ROW, "precondition");
703  mi_math_assert_msg( col < COL, "precondition");
704  return begin()[row * COL + col];
705  }
706 
710  inline T get( Size i) const
711  {
712  mi_math_assert_msg( i < ROW*COL, "precondition");
713  return begin()[i];
714  }
715 
719  inline T get( Size row, Size col) const
720  {
721  mi_math_assert_msg( row < ROW, "precondition");
722  mi_math_assert_msg( col < COL, "precondition");
723  return begin()[row * COL + col];
724  }
725 
730  inline void set( Size i, T value)
731  {
732  mi_math_assert_msg( i < ROW*COL, "precondition");
733  begin()[i] = value;
734  }
735 
740  inline void set( Size row, Size col, T value)
741  {
742  mi_math_assert_msg( row < ROW, "precondition");
743  mi_math_assert_msg( col < COL, "precondition");
744  begin()[row * COL + col] = value;
745  }
746 
750  inline T det33() const
751  {
752  mi_static_assert( (ROW==3 || ROW==4) && (COL==3 || COL==4) );
753  return this->xx * this->yy * this->zz
754  + this->xy * this->yz * this->zx
755  + this->xz * this->yx * this->zy
756  - this->xx * this->zy * this->yz
757  - this->xy * this->zz * this->yx
758  - this->xz * this->zx * this->yy;
759  }
760 
764  inline bool invert();
765 
771  inline void transpose()
772  {
773  mi_static_assert( ROW==COL);
774  for( Size i=0; i < ROW-1; ++i) {
775  for( Size j=i+1; j < COL; ++j) {
776  T tmp = get(i,j);
777  set(i,j, get(j,i));
778  set(j,i, tmp);
779  }
780  }
781  }
782 
786  inline void translate( T x, T y, T z)
787  {
788  this->wx += x;
789  this->wy += y;
790  this->wz += z;
791  }
792 
796  inline void translate( const Vector<Float32,3>& vector)
797  {
798  this->wx += T( vector.x);
799  this->wy += T( vector.y);
800  this->wz += T( vector.z);
801  }
802 
806  inline void translate( const Vector<Float64,3>& vector)
807  {
808  this->wx += T( vector.x);
809  this->wy += T( vector.y);
810  this->wz += T( vector.z);
811  }
812 
816  inline void set_translation( T dx, T dy, T dz)
817  {
818  this->wx = dx;
819  this->wy = dy;
820  this->wz = dz;
821  }
822 
826  inline void set_translation( const Vector<Float32,3>& vector)
827  {
828  this->wx = T( vector.x);
829  this->wy = T( vector.y);
830  this->wz = T( vector.z);
831  }
832 
836  inline void set_translation( const Vector<Float64,3>& vector)
837  {
838  this->wx = T( vector.x);
839  this->wy = T( vector.y);
840  this->wz = T( vector.z);
841  }
842 
847  inline void rotate( T xangle, T yangle, T zangle)
848  {
849  Matrix<T,4,4> tmp( T( 1));
850  tmp.set_rotation( xangle, yangle, zangle);
851  (*this) *= tmp;
852  }
853 
858  inline void rotate( const Vector<Float32,3>& angles)
859  {
860  Matrix<T,4,4> tmp( T( 1));
861  tmp.set_rotation( T( angles.x), T( angles.y), T( angles.z));
862  (*this) *= tmp;
863  }
864 
869  inline void rotate( const Vector<Float64,3>& angles)
870  {
871  Matrix<T,4,4> tmp( T( 1));
872  tmp.set_rotation( T( angles.x), T( angles.y), T( angles.z));
873  (*this) *= tmp;
874  }
875 
881  inline void set_rotation( T x_angle, T y_angle, T z_angle);
882 
888  inline void set_rotation( const Vector<Float32,3>& angles)
889  {
890  set_rotation( T( angles.x), T( angles.y), T( angles.z));
891  }
892 
898  inline void set_rotation( const Vector<Float64,3>& angles)
899  {
900  set_rotation( T( angles.x), T( angles.y), T( angles.z));
901  }
902 
909  inline void set_rotation( const Vector<Float32,3>& axis, Float64 angle);
910 
917  inline void set_rotation( const Vector<Float64,3>& axis, Float64 angle);
918 
923  inline void lookat(
924  const Vector<Float32,3>& position,
925  const Vector<Float32,3>& target,
926  const Vector<Float32,3>& up);
927 
933  inline void lookat(
934  const Vector<Float64,3>& position,
935  const Vector<Float64,3>& target,
936  const Vector<Float64,3>& up);
937 };
938 
939 
940 //------ Free comparison operators ==, !=, <, <=, >, >= for matrices -----------
941 
943 template <typename T, Size ROW, Size COL>
944 inline bool operator==(
945  const Matrix<T,ROW,COL>& lhs,
946  const Matrix<T,ROW,COL>& rhs)
947 {
948  return is_equal( lhs, rhs);
949 }
950 
952 template <typename T, Size ROW, Size COL>
953 inline bool operator!=(
954  const Matrix<T,ROW,COL>& lhs,
955  const Matrix<T,ROW,COL>& rhs)
956 {
957  return is_not_equal( lhs, rhs);
958 }
959 
963 template <typename T, Size ROW, Size COL>
964 inline bool operator<(
965  const Matrix<T,ROW,COL>& lhs,
966  const Matrix<T,ROW,COL>& rhs)
967 {
968  return lexicographically_less( lhs, rhs);
969 }
970 
974 template <typename T, Size ROW, Size COL>
975 inline bool operator<=(
976  const Matrix<T,ROW,COL>& lhs,
977  const Matrix<T,ROW,COL>& rhs)
978 {
979  return lexicographically_less_or_equal( lhs, rhs);
980 }
981 
985 template <typename T, Size ROW, Size COL>
986 inline bool operator>(
987  const Matrix<T,ROW,COL>& lhs,
988  const Matrix<T,ROW,COL>& rhs)
989 {
990  return lexicographically_greater( lhs, rhs);
991 }
992 
996 template <typename T, Size ROW, Size COL>
997 inline bool operator>=(
998  const Matrix<T,ROW,COL>& lhs,
999  const Matrix<T,ROW,COL>& rhs)
1000 {
1001  return lexicographically_greater_or_equal( lhs, rhs);
1002 }
1003 
1004 //------ Operator declarations for Matrix -------------------------------------
1005 
1007 template <typename T, Size ROW, Size COL>
1008 Matrix<T,ROW,COL>& operator+=( Matrix<T,ROW,COL>& lhs, const Matrix<T,ROW,COL>& rhs);
1009 
1011 template <typename T, Size ROW, Size COL>
1012 Matrix<T,ROW,COL>& operator-=( Matrix<T,ROW,COL>& lhs, const Matrix<T,ROW,COL>& rhs);
1013 
1015 template <typename T, Size ROW, Size COL>
1017  const Matrix<T,ROW,COL>& lhs,
1018  const Matrix<T,ROW,COL>& rhs)
1019 {
1020  Matrix<T,ROW,COL> temp( lhs);
1021  temp += rhs;
1022  return temp;
1023 }
1024 
1026 template <typename T, Size ROW, Size COL>
1028  const Matrix<T,ROW,COL>& lhs,
1029  const Matrix<T,ROW,COL>& rhs)
1030 {
1031  Matrix<T,ROW,COL> temp( lhs);
1032  temp -= rhs;
1033  return temp;
1034 }
1035 
1037 template <typename T, Size ROW, Size COL>
1039  const Matrix<T,ROW,COL>& mat)
1040 {
1041  Matrix<T,ROW,COL> temp;
1042  for( Size i(0u); i < ROW*COL; ++i)
1043  temp.begin()[i] = -mat.get(i);
1044  return temp;
1045 }
1046 
1053 template <typename T, Size ROW, Size COL>
1055  Matrix<T,ROW,COL>& lhs,
1056  const Matrix<T,COL,COL>& rhs)
1057 {
1058  // There are more efficient ways of implementing this. Its a default solution. Specializations
1059  // exist below.
1060  Matrix<T,ROW,COL> old( lhs);
1061  for( Size rrow = 0; rrow < ROW; ++rrow) {
1062  for( Size rcol = 0; rcol < COL; ++rcol) {
1063  lhs( rrow, rcol) = T(0);
1064  for( Size k = 0; k < COL; ++k)
1065  lhs( rrow, rcol) += old( rrow, k) * rhs( k, rcol);
1066  }
1067  }
1068  return lhs;
1069 }
1070 
1076 template <typename T, Size ROW1, Size COL1, Size ROW2, Size COL2>
1078  const Matrix<T,ROW1,COL1>& lhs,
1079  const Matrix<T,ROW2,COL2>& rhs)
1080 {
1081  // There are more efficient ways of implementing this. Its a default solution. Specializations
1082  // exist below.
1083  mi_static_assert( COL1 == ROW2);
1084  Matrix<T,ROW1,COL2> result;
1085  for( Size rrow = 0; rrow < ROW1; ++rrow) {
1086  for( Size rcol = 0; rcol < COL2; ++rcol) {
1087  result( rrow, rcol) = T(0);
1088  for( Size k = 0; k < COL1; ++k)
1089  result( rrow, rcol) += lhs( rrow, k) * rhs( k, rcol);
1090  }
1091  }
1092  return result;
1093 }
1094 
1102 template <typename T, Size ROW, Size COL, Size DIM>
1104  const Matrix<T,ROW,COL>& mat,
1105  const Vector<T,DIM>& vec)
1106 {
1107  mi_static_assert( COL == DIM);
1108  Vector<T,ROW> result;
1109  for( Size row = 0; row < ROW; ++row) {
1110  result[row] = T(0);
1111  for( Size col = 0; col < COL; ++col)
1112  result[row] += mat( row, col) * vec[col];
1113  }
1114  return result;
1115 }
1116 
1123 template <Size DIM, typename T, Size ROW, Size COL>
1125  const Vector<T,DIM>& vec,
1126  const Matrix<T,ROW,COL>& mat)
1127 {
1128  mi_static_assert( DIM == ROW);
1129  Vector<T,COL> result;
1130  for( Size col = 0; col < COL; ++col) {
1131  result[col] = T(0);
1132  for( Size row = 0; row < ROW; ++row)
1133  result[col] += mat( row, col) * vec[row];
1134  }
1135  return result;
1136 }
1137 
1140 template <typename T, Size ROW, Size COL>
1142 {
1143  for( Size i=0; i < ROW*COL; ++i)
1144  mat.begin()[i] *= factor;
1145  return mat;
1146 }
1147 
1149 template <typename T, Size ROW, Size COL>
1150 inline Matrix<T,ROW,COL> operator*( const Matrix<T,ROW,COL>& mat, T factor)
1151 {
1152  Matrix<T,ROW,COL> temp( mat);
1153  temp *= factor;
1154  return temp;
1155 }
1156 
1158 template <typename T, Size ROW, Size COL>
1159 inline Matrix<T,ROW,COL> operator*( T factor, const Matrix<T,ROW,COL>& mat)
1160 {
1161  Matrix<T,ROW,COL> temp( mat);
1162  temp *= factor; // * on T should be commutative (IEEE-754)
1163  return temp;
1164 }
1165 
1166 
1167 //------ Free Functions for Matrix --------------------------------------------
1168 
1172 template <Size NEW_ROW, Size NEW_COL, typename T, Size ROW, Size COL>
1174  const Matrix<T,ROW,COL>& mat)
1175 {
1176  mi_static_assert( NEW_ROW <= ROW);
1177  mi_static_assert( NEW_COL <= COL);
1179  for( Size i=0; i < NEW_ROW; ++i)
1180  for( Size j=0; j < NEW_COL; ++j)
1181  result( i, j) = mat( i, j);
1182  return result;
1183 }
1184 
1185 
1186 
1188 template <typename T, Size ROW, Size COL>
1190  const Matrix<T,ROW,COL>& mat)
1191 {
1192  Matrix<T,COL,ROW> result;
1193  for( Size i=0; i < ROW; ++i)
1194  for( Size j=0; j < COL; ++j)
1195  result( j, i) = mat( i, j);
1196  return result;
1197 }
1198 
1204 template <typename T, typename U>
1206  const Matrix<T,4,4>& mat,
1207  const U& point)
1208 {
1209  const T w = T(mat.xw * point + mat.ww);
1210  if( w == T(0) || w == T(1))
1211  return U(mat.xx * point + mat.wx);
1212  else
1213  return U((mat.xx * point + mat.wx) / w);
1214 }
1215 
1221 template <typename T, typename U>
1223  const Matrix<T,4,4>& mat,
1224  const Vector<U,2>& point)
1225 {
1226  T w = T(mat.xw * point.x + mat.yw * point.y + mat.ww);
1227  if( w == T(0) || w == T(1))
1228  return Vector<U, 2>(
1229  U(mat.xx * point.x + mat.yx * point.y + mat.wx),
1230  U(mat.xy * point.x + mat.yy * point.y + mat.wy));
1231  else {
1232  w = T(1)/w; // optimization
1233  return Vector<U, 2>(
1234  U((mat.xx * point.x + mat.yx * point.y + mat.wx) * w),
1235  U((mat.xy * point.x + mat.yy * point.y + mat.wy) * w));
1236  }
1237 }
1238 
1244 template <typename T, typename U>
1246  const Matrix<T,4,3>& mat,
1247  const Vector<U,3>& point)
1248 {
1249  return Vector<U,3>(
1250  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx),
1251  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy),
1252  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz));
1253 }
1254 
1260 template <typename T, typename U>
1262  const Matrix<T,4,4>& mat,
1263  const Vector<U,3>& point)
1264 {
1265  T w = T(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww);
1266  if( w == T(0) || w == T(1)) // avoids temporary and division
1267  return Vector<U,3>(
1268  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx),
1269  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy),
1270  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz));
1271  else {
1272  w = T(1)/w; // optimization
1273  return Vector<U,3>(
1274  U((mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx) * w),
1275  U((mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy) * w),
1276  U((mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz) * w));
1277  }
1278 }
1279 
1285 template <typename T, typename U>
1287  const Matrix<T,4,4>& mat,
1288  const Vector<U,4>& point)
1289 {
1290  return Vector<U, 4>(
1291  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx * point.w),
1292  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy * point.w),
1293  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz * point.w),
1294  U(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww * point.w));
1295 }
1296 
1302 template <typename T, typename U>
1304  const Matrix<T,4,4>& mat,
1305  const U& vector)
1306 {
1307  return U(mat.xx * vector);
1308 }
1309 
1315 template <typename T, typename U>
1317  const Matrix<T,4,4>& mat,
1318  const Vector<U,2>& vector)
1319 {
1320  return Vector<U,2>(
1321  U(mat.xx * vector.x + mat.yx * vector.y),
1322  U(mat.xy * vector.x + mat.yy * vector.y));
1323 }
1324 
1330 template <typename T, typename U>
1332  const Matrix<T,3,3>& mat,
1333  const Vector<U,3>& vector)
1334 {
1335  return Vector<U,3>(
1336  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1337  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1338  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1339 }
1340 
1346 template <typename T, typename U>
1348  const Matrix<T,4,3>& mat,
1349  const Vector<U,3>& vector)
1350 {
1351  return Vector<U,3>(
1352  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1353  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1354  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1355 }
1356 
1362 template <typename T, typename U>
1364  const Matrix<T,4,4>& mat,
1365  const Vector<U,3>& vector)
1366 {
1367  return Vector<U,3>(
1368  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1369  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1370  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1371 }
1372 
1384 template <typename T, typename U>
1386  const Matrix<T,3,3>& inv_mat,
1387  const Vector<U,3>& normal)
1388 {
1389  return Vector<U,3>(
1390  U(inv_mat.xx * normal.x + inv_mat.xy * normal.y + inv_mat.xz * normal.z),
1391  U(inv_mat.yx * normal.x + inv_mat.yy * normal.y + inv_mat.yz * normal.z),
1392  U(inv_mat.zx * normal.x + inv_mat.zy * normal.y + inv_mat.zz * normal.z));
1393 }
1394 
1406 template <typename T, typename U>
1408  const Matrix<T,4,4>& inv_mat,
1409  const Vector<U,3>& normal)
1410 {
1411  return Vector<U,3>(
1412  U(inv_mat.xx * normal.x + inv_mat.xy * normal.y + inv_mat.xz * normal.z),
1413  U(inv_mat.yx * normal.x + inv_mat.yy * normal.y + inv_mat.yz * normal.z),
1414  U(inv_mat.zx * normal.x + inv_mat.zy * normal.y + inv_mat.zz * normal.z));
1415 }
1416 
1430 template <typename T, typename U>
1432  const Matrix<T,4,4>& mat,
1433  const Vector<U,3>& normal)
1434 {
1435  Matrix<T,3,3> sub_mat( mat.xx, mat.xy, mat.xz,
1436  mat.yx, mat.yy, mat.yz,
1437  mat.zx, mat.zy, mat.zz);
1438  bool inverted = sub_mat.invert();
1439  if( !inverted)
1440  return normal;
1441  return Vector<U,3>(
1442  U(sub_mat.xx * normal.x + sub_mat.xy * normal.y + sub_mat.xz * normal.z),
1443  U(sub_mat.yx * normal.x + sub_mat.yy * normal.y + sub_mat.yz * normal.z),
1444  U(sub_mat.zx * normal.x + sub_mat.zy * normal.y + sub_mat.zz * normal.z));
1445 }
1446 
1447 
1448 //------ Definitions of non-inline function -----------------------------------
1449 
1450 template <typename T, Size ROW, Size COL>
1452  Matrix<T,ROW,COL>& lhs,
1453  const Matrix<T,ROW,COL>& rhs)
1454 {
1455  for( Size i=0; i < ROW*COL; ++i)
1456  lhs.begin()[i] += rhs.get(i);
1457  return lhs;
1458 }
1459 
1460 template <typename T, Size ROW, Size COL>
1462  Matrix<T,ROW,COL>& lhs,
1463  const Matrix<T,ROW,COL>& rhs)
1464 {
1465  for( Size i=0; i < ROW*COL; ++i)
1466  lhs.begin()[i] -= rhs.get(i);
1467  return lhs;
1468 }
1469 
1470 #ifndef MI_FOR_DOXYGEN_ONLY
1471 
1472 template <typename T, Size ROW, Size COL>
1473 inline void Matrix<T,ROW,COL>::set_rotation( T xangle, T yangle, T zangle)
1474 {
1475  mi_static_assert( COL == 4 && ROW == 4);
1476  T tsx, tsy, tsz; // sine of [xyz]_angle
1477  T tcx, tcy, tcz; // cosine of [xyz]_angle
1478  T tmp;
1479  const T min_angle = T(0.00024f);
1480 
1481  if( abs( xangle) > min_angle) {
1482  tsx = sin( xangle);
1483  tcx = cos( xangle);
1484  } else {
1485  tsx = xangle;
1486  tcx = T(1);
1487  }
1488  if( abs( yangle) > min_angle) {
1489  tsy = sin( yangle);
1490  tcy = cos( yangle);
1491  } else {
1492  tsy = yangle;
1493  tcy = T(1);
1494  }
1495  if( abs(zangle) > min_angle) {
1496  tsz = sin( zangle);
1497  tcz = cos( zangle);
1498  } else {
1499  tsz = T(zangle);
1500  tcz = T(1);
1501  }
1502  this->xx = tcy * tcz;
1503  this->xy = tcy * tsz;
1504  this->xz = -tsy;
1505 
1506  tmp = tsx * tsy;
1507  this->yx = tmp * tcz - tcx * tsz;
1508  this->yy = tmp * tsz + tcx * tcz;
1509  this->yz = tsx * tcy;
1510 
1511  tmp = tcx * tsy;
1512  this->zx = tmp * tcz + tsx * tsz;
1513  this->zy = tmp * tsz - tsx * tcz;
1514  this->zz = tcx * tcy;
1515 }
1516 
1517 template <typename T, Size ROW, Size COL>
1518 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float32,3>& axis_v, Float64 angle)
1519 {
1520  mi_static_assert( COL == 4 && ROW == 4);
1521  Vector<T,3> axis( axis_v);
1522  const T min_angle = T(0.00024f);
1523 
1524  if( abs( T(angle)) < min_angle) {
1525  T xa = axis.x * T(angle);
1526  T ya = axis.y * T(angle);
1527  T za = axis.z * T(angle);
1528 
1529  this->xx = T(1);
1530  this->xy = za;
1531  this->xz = -ya;
1532  this->xw = T(0);
1533 
1534  this->yx = za;
1535  this->yy = T(1);
1536  this->yz = xa;
1537  this->yw = T(0);
1538 
1539  this->zx = -ya;
1540  this->zy = -xa;
1541  this->zz = T(1);
1542  this->zw = T(0);
1543  } else {
1544  T s = sin( T(angle));
1545  T c = cos( T(angle));
1546  T t = T(1) - c;
1547  T tmp;
1548 
1549  tmp = t * T(axis.x);
1550  this->xx = tmp * T(axis.x) + c;
1551  this->xy = tmp * T(axis.y) + s * T(axis.z);
1552  this->xz = tmp * T(axis.z) - s * T(axis.y);
1553  this->xw = T(0);
1554 
1555  tmp = t * T(axis.y);
1556  this->yx = tmp * T(axis.x) - s * T(axis.z);
1557  this->yy = tmp * T(axis.y) + c;
1558  this->yz = tmp * T(axis.z) + s * T(axis.x);
1559  this->yw = T(0);
1560 
1561  tmp = t * T(axis.z);
1562  this->zx = tmp * T(axis.x) + s * T(axis.y);
1563  this->zy = tmp * T(axis.y) - s * T(axis.x);
1564  this->zz = tmp * T(axis.z) + c;
1565  this->zw = T(0);
1566  }
1567  this->wx = this->wy = this->wz = T(0);
1568  this->ww = T(1);
1569 }
1570 
1571 template <typename T, Size ROW, Size COL>
1572 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float64,3>& axis_v, Float64 angle)
1573 {
1574  mi_static_assert( COL == 4 && ROW == 4);
1575  Vector<T,3> axis( axis_v);
1576  const T min_angle = T(0.00024f);
1577 
1578  if( abs(T(angle)) < min_angle) {
1579  T xa = axis.x * T(angle);
1580  T ya = axis.y * T(angle);
1581  T za = axis.z * T(angle);
1582 
1583  this->xx = T(1);
1584  this->xy = za;
1585  this->xz = -ya;
1586  this->xw = T(0);
1587 
1588  this->yx = za;
1589  this->yy = T(1);
1590  this->yz = xa;
1591  this->yw = T(0);
1592 
1593  this->zx = -ya;
1594  this->zy = -xa;
1595  this->zz = T(1);
1596  this->zw = T(0);
1597  } else {
1598  T s = sin( T(angle));
1599  T c = cos( T(angle));
1600  T t = T(1) - c;
1601  T tmp;
1602 
1603  tmp = t * T(axis.x);
1604  this->xx = tmp * T(axis.x) + c;
1605  this->xy = tmp * T(axis.y) + s * T(axis.z);
1606  this->xz = tmp * T(axis.z) - s * T(axis.y);
1607  this->xw = T(0);
1608 
1609  tmp = t * T(axis.y);
1610  this->yx = tmp * T(axis.x) - s * T(axis.z);
1611  this->yy = tmp * T(axis.y) + c;
1612  this->yz = tmp * T(axis.z) + s * T(axis.x);
1613  this->yw = T(0);
1614 
1615  tmp = t * T(axis.z);
1616  this->zx = tmp * T(axis.x) + s * T(axis.y);
1617  this->zy = tmp * T(axis.y) - s * T(axis.x);
1618  this->zz = tmp * T(axis.z) + c;
1619  this->zw = T(0);
1620  }
1621  this->wx = this->wy = this->wz = T(0);
1622  this->ww = T(1);
1623 }
1624 
1625 template <typename T, Size ROW, Size COL>
1626 inline void Matrix<T,ROW,COL>::lookat(
1627  const Vector<Float32,3>& position,
1628  const Vector<Float32,3>& target,
1629  const Vector<Float32,3>& up)
1630 {
1631  mi_static_assert( COL == 4 && ROW == 4);
1632  Vector<Float32,3> xaxis, yaxis, zaxis;
1633 
1634  // Z vector
1635  zaxis = position - target;
1636  zaxis.normalize();
1637 
1638  // X vector = up cross Z
1639  xaxis = cross( up, zaxis);
1640  xaxis.normalize();
1641 
1642  // Recompute Y = Z cross X
1643  yaxis = cross( zaxis, xaxis);
1644  yaxis.normalize();
1645 
1646  // Build rotation matrix
1647  Matrix<T,4,4> rot(
1648  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1649  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1650  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1651  T(0), T(0), T(0), T(1));
1652 
1653  // Compute the new position by multiplying the inverse position with the rotation matrix
1654  Matrix<T,4,4> trans(
1655  T(1), T(0), T(0), T(0),
1656  T(0), T(1), T(0), T(0),
1657  T(0), T(0), T(1), T(0),
1658  T(-position.x), T(-position.y), T(-position.z), T(1));
1659 
1660  *this = trans * rot;
1661 }
1662 
1663 template <typename T, Size ROW, Size COL>
1664 inline void Matrix<T,ROW,COL>::lookat(
1665  const Vector<Float64,3>& position,
1666  const Vector<Float64,3>& target,
1667  const Vector<Float64,3>& up)
1668 {
1669  mi_static_assert( COL == 4 && ROW == 4);
1670  Vector<Float64,3> xaxis, yaxis, zaxis;
1671 
1672  // Z vector
1673  zaxis = position - target;
1674  zaxis.normalize();
1675 
1676  // X vector = up cross Z
1677  xaxis = cross( up, zaxis);
1678  xaxis.normalize();
1679 
1680  // Recompute Y = Z cross X
1681  yaxis = cross( zaxis, xaxis);
1682  yaxis.normalize();
1683 
1684  // Build rotation matrix
1685  Matrix<T,4,4> rot(
1686  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1687  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1688  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1689  T(0), T(0), T(0), T(1));
1690 
1691  // Compute the new position by multiplying the inverse position with the rotation matrix
1692  Matrix<T,4,4> trans(
1693  T(1), T(0), T(0), T(0),
1694  T(0), T(1), T(0), T(0),
1695  T(0), T(0), T(1), T(0),
1696  T(-position.x), T(-position.y), T(-position.z), T(1));
1697 
1698  *this = trans * rot;
1699 }
1700 
1701 
1702 //------ Definition and helper for matrix inversion ---------------------------
1703 
1704 // Matrix inversion class, specialized for symmetric matrices and low dimensions
1705 template <class T, Size ROW, Size COL>
1706 class Matrix_inverter
1707 {
1708 public:
1709  typedef math::Matrix<T,ROW,COL> Matrix;
1710 
1711  // Inverts the matrix \c M if possible and returns \c true.
1712  //
1713  // If the matrix cannot be inverted or if \c ROW != \c COL, this function returns \c false.
1714  static inline bool invert( Matrix& /* M */) { return false; }
1715 };
1716 
1717 // Matrix inversion class, specialization for symmetric matrices
1718 template <class T, Size DIM>
1719 class Matrix_inverter<T,DIM,DIM>
1720 {
1721 public:
1722  typedef math::Matrix<T,DIM,DIM> Matrix;
1723  typedef math::Vector<T,DIM> Value_vector;
1724  typedef math::Vector<Size,DIM> Index_vector;
1725 
1726  // LU decomposition of matrix lu in place.
1727  //
1728  // Returns \c false if matrix cannot be decomposed, for example, if it is not invertible, and \c
1729  // true otherwise. Returns also a row permutation in indx.
1730  static bool lu_decomposition(
1731  Matrix& lu, // matrix to decompose, modified in place
1732  Index_vector& indx); // row permutation info indx[4]
1733 
1734  // Solves the equation lu * x = b for x. lu and indx are the results of lu_decomposition. The
1735  // solution is returned in b.
1736  static void lu_backsubstitution(
1737  const Matrix& lu, // LU decomposed matrix
1738  const Index_vector& indx, // permutation vector
1739  Value_vector& b); // right hand side vector b, modified in place
1740 
1741  static bool invert( Matrix& mat);
1742 };
1743 
1744 template <class T, Size DIM>
1745 bool Matrix_inverter<T,DIM,DIM>::lu_decomposition(
1746  Matrix& lu,
1747  Index_vector& indx)
1748 {
1749  Value_vector vv; // implicit scaling of each row
1750 
1751  for( Size i = 0; i < DIM; i++) { // get implicit scaling
1752  T big = T(0);
1753  for( Size j = 0; j < DIM; j++) {
1754  T temp = abs(lu.get(i,j));
1755  if( temp > big)
1756  big = temp;
1757  }
1758  if( big == T(0))
1759  return false;
1760  vv[i] = T(1) / big; // save scaling
1761  }
1762  //
1763  // loop over columns of Crout's method
1764  //
1765  Size imax = 0;
1766  for( Size j = 0; j < DIM; j++) {
1767  for( Size i = 0; i < j; i++) {
1768  T sum = lu.get(i,j);
1769  for( Size k = 0; k < i; k++)
1770  sum -= lu.get(i,k) * lu.get(k,j);
1771  lu.set(i, j, sum);
1772  }
1773  T big = 0; // init for pivot search
1774  for( Size i = j; i < DIM; i++) {
1775  T sum = lu.get(i,j);
1776  for( Size k = 0; k < j; k++)
1777  sum -= lu.get(i,k) * lu.get(k,j);
1778  lu.set(i, j, sum);
1779  T dum = vv[i] * abs(sum);
1780  if( dum >= big) { // pivot good?
1781  big = dum;
1782  imax = i;
1783  }
1784  }
1785  if( j != imax) { // interchange rows
1786  for( Size k = 0; k < DIM; k++) {
1787  T dum = lu.get(imax,k);
1788  lu.set(imax, k, lu.get(j,k));
1789  lu.set(j, k, dum);
1790  }
1791  vv[imax] = vv[j]; // interchange scale factor
1792  }
1793  indx[j] = imax;
1794  if( lu.get(j,j) == 0)
1795  return false;
1796  if( j != DIM-1) { // divide by pivot element
1797  T dum = T(1) / lu.get(j,j);
1798  for( Size i = j + 1; i < DIM; i++)
1799  lu.set(i, j, lu.get(i,j) * dum);
1800  }
1801  }
1802  return true;
1803 }
1804 
1805 template <class T, Size DIM>
1806 void Matrix_inverter<T,DIM,DIM>::lu_backsubstitution(
1807  const Matrix& lu,
1808  const Index_vector& indx,
1809  Value_vector& b)
1810 {
1811  // when ii != DIM+1, ii is index of first non-vanishing element of b
1812  Size ii = DIM+1;
1813  for( Size i = 0; i < DIM; i++) {
1814  Size ip = indx[i];
1815  T sum = b[ip];
1816  b[ip] = b[i];
1817  if( ii != DIM+1) {
1818  for( Size j = ii; j < i; j++) {
1819  sum -= lu.get(i,j) * b[j];
1820  }
1821  } else {
1822  if( sum != T(0)) { // non-zero element
1823  ii = i;
1824  }
1825  }
1826  b[i] = sum;
1827  }
1828  for( Size i2 = DIM; i2 > 0;) { // backsubstitution
1829  --i2;
1830  T sum = b[i2];
1831  for( Size j = i2+1; j < DIM; j++)
1832  sum -= lu.get(i2,j) * b[j];
1833  b[i2] = sum / lu.get(i2,i2); // store comp. of sol. vector
1834  }
1835 }
1836 
1837 template <class T, Size DIM>
1838 bool Matrix_inverter<T,DIM,DIM>::invert( Matrix& mat)
1839 {
1840  Matrix lu(mat); // working copy of matrix to invert
1841  Index_vector indx; // row permutation info
1842 
1843  // LU decompose, return if it fails
1844  if( !lu_decomposition(lu, indx))
1845  return false;
1846 
1847  // solve for each column vector of the I matrix
1848  for( Size j = 0; j < DIM; ++j) {
1849  Value_vector col(T(0)); // TODO: do that directly in the mat matrix
1850  col[j] = T(1);
1851  lu_backsubstitution( lu, indx, col);
1852  for( Size i = 0; i < DIM; ++i) {
1853  mat.set( i, j, col[i]);
1854  }
1855  }
1856  return true;
1857 }
1858 
1859 // Matrix inversion class, specialization for 1x1 matrix
1860 template <class T>
1861 class Matrix_inverter<T,1,1>
1862 {
1863 public:
1864  typedef math::Matrix<T,1,1> Matrix;
1865 
1866  static inline bool invert( Matrix& mat)
1867  {
1868  T s = mat.get( 0, 0);
1869  if( s == T(0))
1870  return false;
1871  mat.set( 0, 0, T(1) / s);
1872  return true;
1873  }
1874 };
1875 
1876 // Matrix inversion class, specialization for 2x2 matrix
1877 template <class T>
1878 class Matrix_inverter<T,2,2>
1879 {
1880 public:
1881  typedef math::Matrix<T,2,2> Matrix;
1882 
1883  static inline bool invert( Matrix& mat)
1884  {
1885  T a = mat.get( 0, 0);
1886  T b = mat.get( 0, 1);
1887  T c = mat.get( 1, 0);
1888  T d = mat.get( 1, 1);
1889  T det = a*d - b*c;
1890  if( det == T( 0))
1891  return false;
1892  T rdet = T(1) / det;
1893  mat.set( 0, 0, d * rdet);
1894  mat.set( 0, 1,-b * rdet);
1895  mat.set( 1, 0,-c * rdet);
1896  mat.set( 1, 1, a * rdet);
1897  return true;
1898  }
1899 };
1900 
1901 template <typename T, Size ROW, Size COL>
1902 inline bool Matrix<T,ROW,COL>::invert()
1903 {
1904  return Matrix_inverter<T,ROW,COL>::invert( *this);
1905 }
1906 
1907 
1908 
1909 //------ Specializations and (specialized) overloads --------------------------
1910 
1911 // Specialization of common matrix multiplication for 4x4 matrices.
1912 template <typename T>
1913 inline Matrix<T,4,4>& operator*=(
1914  Matrix<T,4,4>& lhs,
1915  const Matrix<T,4,4>& rhs)
1916 {
1917  Matrix<T,4,4> old( lhs);
1918 
1919  lhs.xx = old.xx * rhs.xx + old.xy * rhs.yx + old.xz * rhs.zx + old.xw * rhs.wx;
1920  lhs.xy = old.xx * rhs.xy + old.xy * rhs.yy + old.xz * rhs.zy + old.xw * rhs.wy;
1921  lhs.xz = old.xx * rhs.xz + old.xy * rhs.yz + old.xz * rhs.zz + old.xw * rhs.wz;
1922  lhs.xw = old.xx * rhs.xw + old.xy * rhs.yw + old.xz * rhs.zw + old.xw * rhs.ww;
1923 
1924  lhs.yx = old.yx * rhs.xx + old.yy * rhs.yx + old.yz * rhs.zx + old.yw * rhs.wx;
1925  lhs.yy = old.yx * rhs.xy + old.yy * rhs.yy + old.yz * rhs.zy + old.yw * rhs.wy;
1926  lhs.yz = old.yx * rhs.xz + old.yy * rhs.yz + old.yz * rhs.zz + old.yw * rhs.wz;
1927  lhs.yw = old.yx * rhs.xw + old.yy * rhs.yw + old.yz * rhs.zw + old.yw * rhs.ww;
1928 
1929  lhs.zx = old.zx * rhs.xx + old.zy * rhs.yx + old.zz * rhs.zx + old.zw * rhs.wx;
1930  lhs.zy = old.zx * rhs.xy + old.zy * rhs.yy + old.zz * rhs.zy + old.zw * rhs.wy;
1931  lhs.zz = old.zx * rhs.xz + old.zy * rhs.yz + old.zz * rhs.zz + old.zw * rhs.wz;
1932  lhs.zw = old.zx * rhs.xw + old.zy * rhs.yw + old.zz * rhs.zw + old.zw * rhs.ww;
1933 
1934  lhs.wx = old.wx * rhs.xx + old.wy * rhs.yx + old.wz * rhs.zx + old.ww * rhs.wx;
1935  lhs.wy = old.wx * rhs.xy + old.wy * rhs.yy + old.wz * rhs.zy + old.ww * rhs.wy;
1936  lhs.wz = old.wx * rhs.xz + old.wy * rhs.yz + old.wz * rhs.zz + old.ww * rhs.wz;
1937  lhs.ww = old.wx * rhs.xw + old.wy * rhs.yw + old.wz * rhs.zw + old.ww * rhs.ww;
1938 
1939  return lhs;
1940 }
1941 
1942 // Specialization of common matrix multiplication for 4x4 matrices.
1943 template <typename T>
1944 inline Matrix<T,4,4> operator*(
1945  const Matrix<T,4,4>& lhs,
1946  const Matrix<T,4,4>& rhs)
1947 {
1948  Matrix<T,4,4> temp( lhs);
1949  temp *= rhs;
1950  return temp;
1951 }
1952 
1953 #endif // MI_FOR_DOXYGEN_ONLY
1954  // end group mi_math_matrix
1956 
1957 } // namespace math
1958 
1959 } // namespace mi
1960 
1961 #endif // MI_MATH_MATRIX_H