NVIDIA Iray: Math API Home  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 1986, 2016 NVIDIA Corporation. All rights reserved.
3 //*****************************************************************************
9 //*****************************************************************************
10 
11 #ifndef MI_MATH_MATRIX_H
12 #define MI_MATH_MATRIX_H
13 
14 #include <mi/base/types.h>
15 #include <mi/math/assert.h>
16 #include <mi/math/function.h>
17 #include <mi/math/vector.h>
18 
19 namespace mi {
20 
21 namespace math {
22 
49 //------- POD struct that provides storage for matrix elements --------
50 
88 template <typename T, Size ROW, Size COL>
89 struct Matrix_struct
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>
297 inline T* matrix_base_ptr( Matrix_struct<T,ROW,COL>& mat)
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 
364 template <typename T, Size ROW, Size COL>
365 class Matrix : public Matrix_struct<T,ROW,COL>
366 {
367 public:
368  typedef Matrix_struct<T,ROW,COL> Pod_type;
369  typedef Matrix_struct<T,ROW,COL> storage_type;
370 
371  typedef T value_type;
372  typedef Size size_type;
373  typedef Difference difference_type;
374  typedef T * pointer;
375  typedef const T * const_pointer;
376  typedef T & reference;
377  typedef const T & const_reference;
378 
380  typedef Vector<T,COL> Row_vector;
381 
383  typedef Vector<T,ROW> Column_vector;
384 
385  static const Size ROWS = ROW;
386  static const Size COLUMNS = COL;
387  static const Size SIZE = ROW*COL;
388 
390  static inline Size size() { return SIZE; }
391 
393  static inline Size max_size() { return SIZE; }
394 
397  enum Transposed_copy_tag {
400  TRANSPOSED_COPY_TAG
401  };
402 
404  inline T * begin() { return mi::math::matrix_base_ptr( *this); }
405 
407  inline T const * begin() const { return mi::math::matrix_base_ptr( *this); }
408 
412  inline T * end() { return begin() + SIZE; }
413 
417  inline T const * end() const { return begin() + SIZE; }
418 
420  inline Row_vector & operator[] (Size row)
421  {
422  mi_math_assert_msg( row < ROW, "precondition");
423  return reinterpret_cast<Row_vector&>(begin()[row * COL]);
424  }
425 
427  inline const Row_vector & operator[] (Size row) const
428  {
429  mi_math_assert_msg( row < ROW, "precondition");
430  return reinterpret_cast<const Row_vector&>(begin()[row * COL]);
431  }
432 
434  inline Matrix() { }
435 
437  inline Matrix( const Matrix_struct<T,ROW,COL>& other)
438  {
439  for( Size i(0u); i < ROW * COL; ++i)
440  begin()[i] = mi::math::matrix_base_ptr( other)[i];
441  }
442 
446  explicit inline Matrix( T diag)
447  {
448  for( Size i(0u); i < SIZE; ++i)
449  begin()[i] = T(0);
450  const Size MIN_DIM = (ROW < COL) ? ROW : COL;
451  for( Size k(0u); k < MIN_DIM; ++k)
452  begin()[k * COL + k] = diag;
453  }
454 
468  template <typename Iterator>
469  inline Matrix( From_iterator_tag, Iterator p)
470  {
471  for( Size i(0u); i < SIZE; ++i, ++p)
472  begin()[i] = *p;
473  }
474 
487  template <typename T2>
488  inline explicit Matrix( T2 const (& array)[SIZE])
489  {
490  for( Size i(0u); i < SIZE; ++i)
491  begin()[i] = array[i];
492  }
493 
496  template <typename T2>
497  inline explicit Matrix( const Matrix<T2,ROW,COL>& other)
498  {
499  for( Size i(0u); i < SIZE; ++i)
500  begin()[i] = T(other.begin()[i]);
501  }
502 
504  inline Matrix(
505  Transposed_copy_tag,
506  const Matrix<T,COL,ROW>& other)
507  {
508  for( Size i(0u); i < ROW; ++i)
509  for( Size j(0u); j < COL; ++j)
510  begin()[i * COL + j] = other.begin()[j * ROW + i];
511  }
512 
516  template <typename T2>
517  inline Matrix(
518  Transposed_copy_tag,
519  const Matrix<T2,COL,ROW>& other)
520  {
521  for( Size i(0u); i < ROW; ++i)
522  for( Size j(0u); j < COL; ++j)
523  begin()[i * COL + j] = T(other.begin()[j * ROW + i]);
524  }
525 
530  inline explicit Matrix(
531  const Row_vector& v0)
532  {
533  mi_static_assert( ROW == 1);
534  (*this)[0] = v0;
535  }
536 
541  inline Matrix(
542  const Row_vector& v0,
543  const Row_vector& v1)
544  {
545  mi_static_assert( ROW == 2);
546  (*this)[0] = v0;
547  (*this)[1] = v1;
548  }
549 
554  inline Matrix(
555  const Row_vector& v0,
556  const Row_vector& v1,
557  const Row_vector& v2)
558  {
559  mi_static_assert( ROW == 3);
560  (*this)[0] = v0;
561  (*this)[1] = v1;
562  (*this)[2] = v2;
563  }
564 
569  inline Matrix(
570  const Row_vector& v0,
571  const Row_vector& v1,
572  const Row_vector& v2,
573  const Row_vector& v3)
574  {
575  mi_static_assert( ROW == 4);
576  (*this)[0] = v0;
577  (*this)[1] = v1;
578  (*this)[2] = v2;
579  (*this)[3] = v3;
580  }
581 
585  inline Matrix( T m0, T m1)
586  {
587  mi_static_assert( SIZE == 2);
588  begin()[0] = m0; begin()[1] = m1;
589  }
590 
594  inline Matrix( T m0, T m1, T m2)
595  {
596  mi_static_assert( SIZE == 3);
597  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2;
598  }
599 
603  inline Matrix( T m0, T m1, T m2, T m3)
604  {
605  mi_static_assert( SIZE == 4);
606  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
607  }
608 
612  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5)
613  {
614  mi_static_assert( SIZE == 6);
615  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
616  begin()[4] = m4; begin()[5] = m5;
617  }
618 
622  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5, T m6, T m7)
623  {
624  mi_static_assert( SIZE == 8);
625  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
626  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
627  }
628 
632  inline Matrix( T m0, T m1, T m2, T m3, T m4, T m5, T m6, T m7, T m8)
633  {
634  mi_static_assert( SIZE == 9);
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  begin()[8] = m8;
638  }
639 
643  inline Matrix(
644  T m0, T m1, T m2, T m3,
645  T m4, T m5, T m6, T m7,
646  T m8, T m9, T m10, T m11)
647  {
648  mi_static_assert( SIZE == 12);
649  begin()[0] = m0; begin()[1] = m1; begin()[2] = m2; begin()[3] = m3;
650  begin()[4] = m4; begin()[5] = m5; begin()[6] = m6; begin()[7] = m7;
651  begin()[8] = m8; begin()[9] = m9; begin()[10] = m10; begin()[11] = m11;
652  }
653 
657  inline Matrix(
658  T m0, T m1, T m2, T m3,
659  T m4, T m5, T m6, T m7,
660  T m8, T m9, T m10, T m11,
661  T m12, T m13, T m14, T m15)
662  {
663  mi_static_assert( SIZE == 16);
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  begin()[12] = m12; begin()[13] = m13; begin()[14] = m14; begin()[15] = m15;
668  }
669 
671  inline Matrix& operator=( const Matrix& other)
672  {
673  Matrix_struct<T,ROW,COL>::operator=( other);
674  return *this;
675  }
676 
680  inline T& operator()( Size row, Size col)
681  {
682  mi_math_assert_msg( row < ROW, "precondition");
683  mi_math_assert_msg( col < COL, "precondition");
684  return begin()[row * COL + col];
685  }
686 
690  inline const T& operator()( Size row, Size col) const
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 T get( Size i) const
701  {
702  mi_math_assert_msg( i < (ROW*COL), "precondition");
703  return begin()[i];
704  }
705 
709  inline T get( Size row, Size col) const
710  {
711  mi_math_assert_msg( row < ROW, "precondition");
712  mi_math_assert_msg( col < COL, "precondition");
713  return begin()[row * COL + col];
714  }
715 
720  inline void set( Size i, T value)
721  {
722  mi_math_assert_msg( i < (ROW*COL), "precondition");
723  begin()[i] = value;
724  }
725 
730  inline void set( Size row, Size col, T value)
731  {
732  mi_math_assert_msg( row < ROW, "precondition");
733  mi_math_assert_msg( col < COL, "precondition");
734  begin()[row * COL + col] = value;
735  }
736 
740  inline T det33() const
741  {
742  mi_static_assert( (ROW==3 || ROW==4) && (COL==3 || COL==4) );
743  return this->xx * this->yy * this->zz
744  + this->xy * this->yz * this->zx
745  + this->xz * this->yx * this->zy
746  - this->xx * this->zy * this->yz
747  - this->xy * this->zz * this->yx
748  - this->xz * this->zx * this->yy;
749  }
750 
754  inline bool invert();
755 
761  inline void transpose()
762  {
763  mi_static_assert( ROW==COL);
764  for( Size i=0; i < ROW-1; ++i) {
765  for( Size j=i+1; j < COL; ++j) {
766  T tmp = get(i,j);
767  set(i,j, get(j,i));
768  set(j,i, tmp);
769  }
770  }
771  }
772 
776  inline void translate( T x, T y, T z)
777  {
778  this->wx += x;
779  this->wy += y;
780  this->wz += z;
781  }
782 
786  inline void translate( const Vector<Float32,3>& vector)
787  {
788  this->wx += T( vector.x);
789  this->wy += T( vector.y);
790  this->wz += T( vector.z);
791  }
792 
796  inline void translate( const Vector<Float64,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 set_translation( T dx, T dy, T dz)
807  {
808  this->wx = dx;
809  this->wy = dy;
810  this->wz = dz;
811  }
812 
816  inline void set_translation( const Vector<Float32,3>& vector)
817  {
818  this->wx = T( vector.x);
819  this->wy = T( vector.y);
820  this->wz = T( vector.z);
821  }
822 
826  inline void set_translation( const Vector<Float64,3>& vector)
827  {
828  this->wx = T( vector.x);
829  this->wy = T( vector.y);
830  this->wz = T( vector.z);
831  }
832 
837  inline void rotate( T xangle, T yangle, T zangle)
838  {
839  Matrix<T,4,4> tmp( T( 1));
840  tmp.set_rotation( xangle, yangle, zangle);
841  (*this) *= tmp;
842  }
843 
848  inline void rotate( const Vector<Float32,3>& angles)
849  {
850  Matrix<T,4,4> tmp( T( 1));
851  tmp.set_rotation( T( angles.x), T( angles.y), T( angles.z));
852  (*this) *= tmp;
853  }
854 
859  inline void rotate( const Vector<Float64,3>& angles)
860  {
861  Matrix<T,4,4> tmp( T( 1));
862  tmp.set_rotation( T( angles.x), T( angles.y), T( angles.z));
863  (*this) *= tmp;
864  }
865 
871  inline void set_rotation( T x_angle, T y_angle, T z_angle);
872 
878  inline void set_rotation( const Vector<Float32,3>& angles)
879  {
880  set_rotation( T( angles.x), T( angles.y), T( angles.z));
881  }
882 
888  inline void set_rotation( const Vector<Float64,3>& angles)
889  {
890  set_rotation( T( angles.x), T( angles.y), T( angles.z));
891  }
892 
899  inline void set_rotation( const Vector<Float32,3>& axis, Float64 angle);
900 
907  inline void set_rotation( const Vector<Float64,3>& axis, Float64 angle);
908 
913  inline void lookat(
914  const Vector<Float32,3>& position,
915  const Vector<Float32,3>& target,
916  const Vector<Float32,3>& up);
917 
923  inline void lookat(
924  const Vector<Float64,3>& position,
925  const Vector<Float64,3>& target,
926  const Vector<Float64,3>& up);
927 };
928 
929 
930 //------ Free comparison operators ==, !=, <, <=, >, >= for matrices -----------
931 
933 template <typename T, Size ROW, Size COL>
934 inline bool operator==(
935  const Matrix<T,ROW,COL>& lhs,
936  const Matrix<T,ROW,COL>& rhs)
937 {
938  return is_equal( lhs, rhs);
939 }
940 
942 template <typename T, Size ROW, Size COL>
943 inline bool operator!=(
944  const Matrix<T,ROW,COL>& lhs,
945  const Matrix<T,ROW,COL>& rhs)
946 {
947  return is_not_equal( lhs, rhs);
948 }
949 
953 template <typename T, Size ROW, Size COL>
954 inline bool operator<(
955  const Matrix<T,ROW,COL>& lhs,
956  const Matrix<T,ROW,COL>& rhs)
957 {
958  return lexicographically_less( lhs, rhs);
959 }
960 
964 template <typename T, Size ROW, Size COL>
965 inline bool operator<=(
966  const Matrix<T,ROW,COL>& lhs,
967  const Matrix<T,ROW,COL>& rhs)
968 {
969  return lexicographically_less_or_equal( lhs, rhs);
970 }
971 
975 template <typename T, Size ROW, Size COL>
976 inline bool operator>(
977  const Matrix<T,ROW,COL>& lhs,
978  const Matrix<T,ROW,COL>& rhs)
979 {
980  return lexicographically_greater( lhs, rhs);
981 }
982 
986 template <typename T, Size ROW, Size COL>
987 inline bool operator>=(
988  const Matrix<T,ROW,COL>& lhs,
989  const Matrix<T,ROW,COL>& rhs)
990 {
991  return lexicographically_greater_or_equal( lhs, rhs);
992 }
993 
994 //------ Operator declarations for Matrix -------------------------------------
995 
997 template <typename T, Size ROW, Size COL>
998 Matrix<T,ROW,COL>& operator+=( Matrix<T,ROW,COL>& lhs, const Matrix<T,ROW,COL>& rhs);
999 
1001 template <typename T, Size ROW, Size COL>
1002 Matrix<T,ROW,COL>& operator-=( Matrix<T,ROW,COL>& lhs, const Matrix<T,ROW,COL>& rhs);
1003 
1005 template <typename T, Size ROW, Size COL>
1006 inline Matrix<T,ROW,COL> operator+(
1007  const Matrix<T,ROW,COL>& lhs,
1008  const Matrix<T,ROW,COL>& rhs)
1009 {
1010  Matrix<T,ROW,COL> temp( lhs);
1011  temp += rhs;
1012  return temp;
1013 }
1014 
1016 template <typename T, Size ROW, Size COL>
1017 inline Matrix<T,ROW,COL> operator-(
1018  const Matrix<T,ROW,COL>& lhs,
1019  const Matrix<T,ROW,COL>& rhs)
1020 {
1021  Matrix<T,ROW,COL> temp( lhs);
1022  temp -= rhs;
1023  return temp;
1024 }
1025 
1027 template <typename T, Size ROW, Size COL>
1028 inline Matrix<T,ROW,COL> operator-(
1029  const Matrix<T,ROW,COL>& mat)
1030 {
1031  Matrix<T,ROW,COL> temp;
1032  for( Size i(0u); i < (ROW*COL); ++i)
1033  temp.begin()[i] = -mat.get(i);
1034  return temp;
1035 }
1036 
1043 template <typename T, Size ROW, Size COL>
1044 inline Matrix<T,ROW,COL>& operator*=(
1045  Matrix<T,ROW,COL>& lhs,
1046  const Matrix<T,COL,COL>& rhs)
1047 {
1048  // There are more efficient ways of implementing this. Its a default solution. Specializations
1049  // exist below.
1050  Matrix<T,ROW,COL> old( lhs);
1051  for( Size rrow = 0; rrow < ROW; ++rrow) {
1052  for( Size rcol = 0; rcol < COL; ++rcol) {
1053  lhs( rrow, rcol) = T(0);
1054  for( Size k = 0; k < COL; ++k)
1055  lhs( rrow, rcol) += old( rrow, k) * rhs( k, rcol);
1056  }
1057  }
1058  return lhs;
1059 }
1060 
1066 template <typename T, Size ROW1, Size COL1, Size ROW2, Size COL2>
1067 inline Matrix<T,ROW1,COL2> operator*(
1068  const Matrix<T,ROW1,COL1>& lhs,
1069  const Matrix<T,ROW2,COL2>& rhs)
1070 {
1071  // There are more efficient ways of implementing this. Its a default solution. Specializations
1072  // exist below.
1073  mi_static_assert( COL1 == ROW2);
1074  Matrix<T,ROW1,COL2> result;
1075  for( Size rrow = 0; rrow < ROW1; ++rrow) {
1076  for( Size rcol = 0; rcol < COL2; ++rcol) {
1077  result( rrow, rcol) = T(0);
1078  for( Size k = 0; k < COL1; ++k)
1079  result( rrow, rcol) += lhs( rrow, k) * rhs( k, rcol);
1080  }
1081  }
1082  return result;
1083 }
1084 
1092 template <typename T, Size ROW, Size COL, Size DIM>
1093 inline Vector<T,ROW> operator*(
1094  const Matrix<T,ROW,COL>& mat,
1095  const Vector<T,DIM>& vec)
1096 {
1097  mi_static_assert( COL == DIM);
1098  Vector<T,ROW> result;
1099  for( Size row = 0; row < ROW; ++row) {
1100  result[row] = T(0);
1101  for( Size col = 0; col < COL; ++col)
1102  result[row] += mat( row, col) * vec[col];
1103  }
1104  return result;
1105 }
1106 
1113 template <Size DIM, typename T, Size ROW, Size COL>
1114 inline Vector<T,COL> operator*(
1115  const Vector<T,DIM>& vec,
1116  const Matrix<T,ROW,COL>& mat)
1117 {
1118  mi_static_assert( DIM == ROW);
1119  Vector<T,COL> result;
1120  for( Size col = 0; col < COL; ++col) {
1121  result[col] = T(0);
1122  for( Size row = 0; row < ROW; ++row)
1123  result[col] += mat( row, col) * vec[row];
1124  }
1125  return result;
1126 }
1127 
1130 template <typename T, Size ROW, Size COL>
1131 inline Matrix<T,ROW,COL>& operator*=( Matrix<T,ROW,COL>& mat, T factor)
1132 {
1133  for( Size i=0; i < (ROW*COL); ++i)
1134  mat.begin()[i] *= factor;
1135  return mat;
1136 }
1137 
1139 template <typename T, Size ROW, Size COL>
1140 inline Matrix<T,ROW,COL> operator*( const Matrix<T,ROW,COL>& mat, T factor)
1141 {
1142  Matrix<T,ROW,COL> temp( mat);
1143  temp *= factor;
1144  return temp;
1145 }
1146 
1148 template <typename T, Size ROW, Size COL>
1149 inline Matrix<T,ROW,COL> operator*( T factor, const Matrix<T,ROW,COL>& mat)
1150 {
1151  Matrix<T,ROW,COL> temp( mat);
1152  temp *= factor; // * on T should be commutative (IEEE-754)
1153  return temp;
1154 }
1155 
1156 
1157 //------ Free Functions for Matrix --------------------------------------------
1158 
1162 template <Size NEW_ROW, Size NEW_COL, typename T, Size ROW, Size COL>
1163 inline Matrix<T,NEW_ROW,NEW_COL> sub_matrix(
1164  const Matrix<T,ROW,COL>& mat)
1165 {
1166  mi_static_assert( NEW_ROW <= ROW);
1167  mi_static_assert( NEW_COL <= COL);
1168  Matrix<T,NEW_ROW,NEW_COL> result;
1169  for( Size i=0; i < NEW_ROW; ++i)
1170  for( Size j=0; j < NEW_COL; ++j)
1171  result( i, j) = mat( i, j);
1172  return result;
1173 }
1174 
1175 
1176 
1178 template <typename T, Size ROW, Size COL>
1179 inline Matrix<T,COL,ROW> transpose(
1180  const Matrix<T,ROW,COL>& mat)
1181 {
1182  Matrix<T,COL,ROW> result;
1183  for( Size i=0; i < ROW; ++i)
1184  for( Size j=0; j < COL; ++j)
1185  result( j, i) = mat( i, j);
1186  return result;
1187 }
1188 
1194 template <typename T, typename U>
1195 inline U transform_point(
1196  const Matrix<T,4,4>& mat,
1197  const U& point)
1198 {
1199  const T w = T(mat.xw * point + mat.ww);
1200  if( w == T(0) || w == T(1))
1201  return U(mat.xx * point + mat.wx);
1202  else
1203  return U((mat.xx * point + mat.wx) / w);
1204 }
1205 
1211 template <typename T, typename U>
1212 inline Vector<U,2> transform_point(
1213  const Matrix<T,4,4>& mat,
1214  const Vector<U,2>& point)
1215 {
1216  T w = T(mat.xw * point.x + mat.yw * point.y + mat.ww);
1217  if( w == T(0) || w == T(1))
1218  return Vector<U, 2>(
1219  U(mat.xx * point.x + mat.yx * point.y + mat.wx),
1220  U(mat.xy * point.x + mat.yy * point.y + mat.wy));
1221  else {
1222  w = T(1)/w; // optimization
1223  return Vector<U, 2>(
1224  U((mat.xx * point.x + mat.yx * point.y + mat.wx) * w),
1225  U((mat.xy * point.x + mat.yy * point.y + mat.wy) * w));
1226  }
1227 }
1228 
1234 template <typename T, typename U>
1235 inline Vector<U,3> transform_point(
1236  const Matrix<T,4,4>& mat,
1237  const Vector<U,3>& point)
1238 {
1239  T w = T(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww);
1240  if( w == T(0) || w == T(1)) // avoids temporary and division
1241  return Vector<U,3>(
1242  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx),
1243  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy),
1244  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz));
1245  else {
1246  w = T(1)/w; // optimization
1247  return Vector<U,3>(
1248  U((mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx) * w),
1249  U((mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy) * w),
1250  U((mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz) * w));
1251  }
1252 }
1253 
1259 template <typename T, typename U>
1260 inline Vector<U,4> transform_point(
1261  const Matrix<T,4,4>& mat,
1262  const Vector<U,4>& point)
1263 {
1264  return Vector<U, 4>(
1265  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx * point.w),
1266  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy * point.w),
1267  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz * point.w),
1268  U(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww * point.w));
1269 }
1270 
1276 template <typename T, typename U>
1277 inline U transform_vector(
1278  const Matrix<T,4,4>& mat,
1279  const U& vector)
1280 {
1281  return U(mat.xx * vector);
1282 }
1283 
1289 template <typename T, typename U>
1290 inline Vector<U,2> transform_vector(
1291  const Matrix<T,4,4>& mat,
1292  const Vector<U,2>& vector)
1293 {
1294  return Vector<U,2>(
1295  U(mat.xx * vector.x + mat.yx * vector.y),
1296  U(mat.xy * vector.x + mat.yy * vector.y));
1297 }
1298 
1304 template <typename T, typename U>
1305 inline Vector<U,3> transform_vector(
1306  const Matrix<T,4,4>& mat,
1307  const Vector<U,3>& vector)
1308 {
1309  return Vector<U,3>(
1310  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1311  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1312  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1313 }
1314 
1326 template <typename T, typename U>
1327 inline Vector<U,3> transform_normal_inv(
1328  const Matrix<T,4,4>& inv_mat,
1329  const Vector<U,3>& normal)
1330 {
1331  return Vector<U,3>(
1332  U(inv_mat.xx * normal.x + inv_mat.xy * normal.y + inv_mat.xz * normal.z),
1333  U(inv_mat.yx * normal.x + inv_mat.yy * normal.y + inv_mat.yz * normal.z),
1334  U(inv_mat.zx * normal.x + inv_mat.zy * normal.y + inv_mat.zz * normal.z));
1335 }
1336 
1350 template <typename T, typename U>
1351 inline Vector<U,3> transform_normal(
1352  const Matrix<T,4,4>& mat,
1353  const Vector<U,3>& normal)
1354 {
1355  Matrix<T,3,3> sub_mat( mat.xx, mat.xy, mat.xz,
1356  mat.yx, mat.yy, mat.yz,
1357  mat.zx, mat.zy, mat.zz);
1358  bool inverted = sub_mat.invert();
1359  if( !inverted)
1360  return normal;
1361  return Vector<U,3>(
1362  U(sub_mat.xx * normal.x + sub_mat.xy * normal.y + sub_mat.xz * normal.z),
1363  U(sub_mat.yx * normal.x + sub_mat.yy * normal.y + sub_mat.yz * normal.z),
1364  U(sub_mat.zx * normal.x + sub_mat.zy * normal.y + sub_mat.zz * normal.z));
1365 }
1366 
1367 
1368 //------ Definitions of non-inline function -----------------------------------
1369 
1370 template <typename T, Size ROW, Size COL>
1371 Matrix<T,ROW,COL>& operator+=(
1372  Matrix<T,ROW,COL>& lhs,
1373  const Matrix<T,ROW,COL>& rhs)
1374 {
1375  for( Size i=0; i < ROW*COL; ++i)
1376  lhs.begin()[i] += rhs.get(i);
1377  return lhs;
1378 }
1379 
1380 template <typename T, Size ROW, Size COL>
1381 Matrix<T,ROW,COL>& operator-=(
1382  Matrix<T,ROW,COL>& lhs,
1383  const Matrix<T,ROW,COL>& rhs)
1384 {
1385  for( Size i=0; i < ROW*COL; ++i)
1386  lhs.begin()[i] -= rhs.get(i);
1387  return lhs;
1388 }
1389 
1390 #ifndef MI_FOR_DOXYGEN_ONLY
1391 
1392 template <typename T, Size ROW, Size COL>
1393 inline void Matrix<T,ROW,COL>::set_rotation( T xangle, T yangle, T zangle)
1394 {
1395  mi_static_assert( COL == 4 && ROW == 4);
1396  T tsx, tsy, tsz; // sine of [xyz]_angle
1397  T tcx, tcy, tcz; // cosine of [xyz]_angle
1398  T tmp;
1399  const T min_angle = T(0.00024f);
1400 
1401  if( abs( xangle) > min_angle) {
1402  tsx = sin( xangle);
1403  tcx = cos( xangle);
1404  } else {
1405  tsx = xangle;
1406  tcx = T(1);
1407  }
1408  if( abs( yangle) > min_angle) {
1409  tsy = sin( yangle);
1410  tcy = cos( yangle);
1411  } else {
1412  tsy = yangle;
1413  tcy = T(1);
1414  }
1415  if( abs(zangle) > min_angle) {
1416  tsz = sin( zangle);
1417  tcz = cos( zangle);
1418  } else {
1419  tsz = T(zangle);
1420  tcz = T(1);
1421  }
1422  this->xx = tcy * tcz;
1423  this->xy = tcy * tsz;
1424  this->xz = -tsy;
1425 
1426  tmp = tsx * tsy;
1427  this->yx = tmp * tcz - tcx * tsz;
1428  this->yy = tmp * tsz + tcx * tcz;
1429  this->yz = tsx * tcy;
1430 
1431  tmp = tcx * tsy;
1432  this->zx = tmp * tcz + tsx * tsz;
1433  this->zy = tmp * tsz - tsx * tcz;
1434  this->zz = tcx * tcy;
1435 }
1436 
1437 template <typename T, Size ROW, Size COL>
1438 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float32,3>& axis_v, Float64 angle)
1439 {
1440  mi_static_assert( COL == 4 && ROW == 4);
1441  Vector<T,3> axis( axis_v);
1442  const T min_angle = T(0.00024f);
1443 
1444  if( abs( T(angle)) < min_angle) {
1445  T xa = axis.x * T(angle);
1446  T ya = axis.y * T(angle);
1447  T za = axis.z * T(angle);
1448 
1449  this->xx = T(1);
1450  this->xy = za;
1451  this->xz = -ya;
1452  this->xw = T(0);
1453 
1454  this->yx = za;
1455  this->yy = T(1);
1456  this->yz = xa;
1457  this->yw = T(0);
1458 
1459  this->zx = -ya;
1460  this->zy = -xa;
1461  this->zz = T(1);
1462  this->zw = T(0);
1463  } else {
1464  T s = sin( T(angle));
1465  T c = cos( T(angle));
1466  T t = T(1) - c;
1467  T tmp;
1468 
1469  tmp = t * T(axis.x);
1470  this->xx = tmp * T(axis.x) + c;
1471  this->xy = tmp * T(axis.y) + s * T(axis.z);
1472  this->xz = tmp * T(axis.z) - s * T(axis.y);
1473  this->xw = T(0);
1474 
1475  tmp = t * T(axis.y);
1476  this->yx = tmp * T(axis.x) - s * T(axis.z);
1477  this->yy = tmp * T(axis.y) + c;
1478  this->yz = tmp * T(axis.z) + s * T(axis.x);
1479  this->yw = T(0);
1480 
1481  tmp = t * T(axis.z);
1482  this->zx = tmp * T(axis.x) + s * T(axis.y);
1483  this->zy = tmp * T(axis.y) - s * T(axis.x);
1484  this->zz = tmp * T(axis.z) + c;
1485  this->zw = T(0);
1486  }
1487  this->wx = this->wy = this->wz = T(0);
1488  this->ww = T(1);
1489 }
1490 
1491 template <typename T, Size ROW, Size COL>
1492 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float64,3>& axis_v, Float64 angle)
1493 {
1494  mi_static_assert( COL == 4 && ROW == 4);
1495  Vector<T,3> axis( axis_v);
1496  const T min_angle = T(0.00024f);
1497 
1498  if( abs(T(angle)) < min_angle) {
1499  T xa = axis.x * T(angle);
1500  T ya = axis.y * T(angle);
1501  T za = axis.z * T(angle);
1502 
1503  this->xx = T(1);
1504  this->xy = za;
1505  this->xz = -ya;
1506  this->xw = T(0);
1507 
1508  this->yx = za;
1509  this->yy = T(1);
1510  this->yz = xa;
1511  this->yw = T(0);
1512 
1513  this->zx = -ya;
1514  this->zy = -xa;
1515  this->zz = T(1);
1516  this->zw = T(0);
1517  } else {
1518  T s = sin( T(angle));
1519  T c = cos( T(angle));
1520  T t = T(1) - c;
1521  T tmp;
1522 
1523  tmp = t * T(axis.x);
1524  this->xx = tmp * T(axis.x) + c;
1525  this->xy = tmp * T(axis.y) + s * T(axis.z);
1526  this->xz = tmp * T(axis.z) - s * T(axis.y);
1527  this->xw = T(0);
1528 
1529  tmp = t * T(axis.y);
1530  this->yx = tmp * T(axis.x) - s * T(axis.z);
1531  this->yy = tmp * T(axis.y) + c;
1532  this->yz = tmp * T(axis.z) + s * T(axis.x);
1533  this->yw = T(0);
1534 
1535  tmp = t * T(axis.z);
1536  this->zx = tmp * T(axis.x) + s * T(axis.y);
1537  this->zy = tmp * T(axis.y) - s * T(axis.x);
1538  this->zz = tmp * T(axis.z) + c;
1539  this->zw = T(0);
1540  }
1541  this->wx = this->wy = this->wz = T(0);
1542  this->ww = T(1);
1543 }
1544 
1545 template <typename T, Size ROW, Size COL>
1546 inline void Matrix<T,ROW,COL>::lookat(
1547  const Vector<Float32,3>& position,
1548  const Vector<Float32,3>& target,
1549  const Vector<Float32,3>& up)
1550 {
1551  mi_static_assert( COL == 4 && ROW == 4);
1552  Vector<Float32,3> xaxis, yaxis, zaxis;
1553 
1554  // Z vector
1555  zaxis = position - target;
1556  zaxis.normalize();
1557 
1558  // X vector = up cross Z
1559  xaxis = cross( up, zaxis);
1560  xaxis.normalize();
1561 
1562  // Recompute Y = Z cross X
1563  yaxis = cross( zaxis, xaxis);
1564  yaxis.normalize();
1565 
1566  // Build rotation matrix
1567  Matrix<T,4,4> rot(
1568  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1569  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1570  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1571  T(0), T(0), T(0), T(1));
1572 
1573  // Compute the new position by multiplying the inverse position with the rotation matrix
1574  Matrix<T,4,4> trans(
1575  T(1), T(0), T(0), T(0),
1576  T(0), T(1), T(0), T(0),
1577  T(0), T(0), T(1), T(0),
1578  T(-position.x), T(-position.y), T(-position.z), T(1));
1579 
1580  *this = trans * rot;
1581 }
1582 
1583 template <typename T, Size ROW, Size COL>
1584 inline void Matrix<T,ROW,COL>::lookat(
1585  const Vector<Float64,3>& position,
1586  const Vector<Float64,3>& target,
1587  const Vector<Float64,3>& up)
1588 {
1589  mi_static_assert( COL == 4 && ROW == 4);
1590  Vector<Float64,3> xaxis, yaxis, zaxis;
1591 
1592  // Z vector
1593  zaxis = position - target;
1594  zaxis.normalize();
1595 
1596  // X vector = up cross Z
1597  xaxis = cross( up, zaxis);
1598  xaxis.normalize();
1599 
1600  // Recompute Y = Z cross X
1601  yaxis = cross( zaxis, xaxis);
1602  yaxis.normalize();
1603 
1604  // Build rotation matrix
1605  Matrix<T,4,4> rot(
1606  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1607  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1608  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1609  T(0), T(0), T(0), T(1));
1610 
1611  // Compute the new position by multiplying the inverse position with the rotation matrix
1612  Matrix<T,4,4> trans(
1613  T(1), T(0), T(0), T(0),
1614  T(0), T(1), T(0), T(0),
1615  T(0), T(0), T(1), T(0),
1616  T(-position.x), T(-position.y), T(-position.z), T(1));
1617 
1618  *this = trans * rot;
1619 }
1620 
1621 
1622 //------ Definition and helper for matrix inversion ---------------------------
1623 
1624 // Matrix inversion class, specialized for symmetric matrices and low dimensions
1625 template <class T, Size ROW, Size COL>
1626 class Matrix_inverter
1627 {
1628 public:
1629  typedef math::Matrix<T,ROW,COL> Matrix;
1630 
1631  // Inverts the matrix \c M if possible and returns \c true.
1632  //
1633  // If the matrix cannot be inverted or if \c ROW != \c COL, this function returns \c false.
1634  static inline bool invert( Matrix& /* M */) { return false; }
1635 };
1636 
1637 // Matrix inversion class, specialization for symmetric matrices
1638 template <class T, Size DIM>
1639 class Matrix_inverter<T,DIM,DIM>
1640 {
1641 public:
1642  typedef math::Matrix<T,DIM,DIM> Matrix;
1643  typedef math::Vector<T,DIM> Vector;
1644  typedef math::Vector<Size,DIM> Index_vector;
1645 
1646  // LU decomposition of matrix lu in place.
1647  //
1648  // Returns \c false if matrix cannot be decomposed, for example, if it is not invertible, and \c
1649  // true otherwise. Returns also a row permutation in indx.
1650  static bool lu_decomposition(
1651  Matrix& lu, // matrix to decompose, modified in place
1652  Index_vector& indx); // row permutation info indx[4]
1653 
1654  // Solves the equation lu * x = b for x. lu and indx are the results of lu_decomposition. The
1655  // solution is returned in b.
1656  static void lu_backsubstitution(
1657  const Matrix& lu, // LU decomposed matrix
1658  const Index_vector& indx, // permutation vector
1659  Vector& b); // right hand side vector b, modified in place
1660 
1661  static bool invert( Matrix& mat);
1662 };
1663 
1664 template <class T, Size DIM>
1665 bool Matrix_inverter<T,DIM,DIM>::lu_decomposition(
1666  Matrix& lu,
1667  Index_vector& indx)
1668 {
1669  Vector vv; // implicit scaling of each row
1670 
1671  for( Size i = 0; i < DIM; i++) { // get implicit scaling
1672  T big = T(0);
1673  for( Size j = 0; j < DIM; j++) {
1674  T temp = abs(lu.get(i,j));
1675  if( temp > big)
1676  big = temp;
1677  }
1678  if( big == T(0))
1679  return false;
1680  vv[i] = T(1) / big; // save scaling
1681  }
1682  //
1683  // loop over columns of Crout's method
1684  //
1685  Size imax = 0;
1686  for( Size j = 0; j < DIM; j++) {
1687  for( Size i = 0; i < j; i++) {
1688  T sum = lu.get(i,j);
1689  for( Size k = 0; k < i; k++)
1690  sum -= lu.get(i,k) * lu.get(k,j);
1691  lu.set(i, j, sum);
1692  }
1693  T big = 0; // init for pivot search
1694  for( Size i = j; i < DIM; i++) {
1695  T sum = lu.get(i,j);
1696  for( Size k = 0; k < j; k++)
1697  sum -= lu.get(i,k) * lu.get(k,j);
1698  lu.set(i, j, sum);
1699  T dum = vv[i] * abs(sum);
1700  if( dum >= big) { // pivot good?
1701  big = dum;
1702  imax = i;
1703  }
1704  }
1705  if( j != imax) { // interchange rows
1706  for( Size k = 0; k < DIM; k++) {
1707  T dum = lu.get(imax,k);
1708  lu.set(imax, k, lu.get(j,k));
1709  lu.set(j, k, dum);
1710  }
1711  vv[imax] = vv[j]; // interchange scale factor
1712  }
1713  indx[j] = imax;
1714  if( lu.get(j,j) == 0)
1715  return false;
1716  if( j != DIM-1) { // divide by pivot element
1717  T dum = T(1) / lu.get(j,j);
1718  for( Size i = j + 1; i < DIM; i++)
1719  lu.set(i, j, lu.get(i,j) * dum);
1720  }
1721  }
1722  return true;
1723 }
1724 
1725 template <class T, Size DIM>
1726 void Matrix_inverter<T,DIM,DIM>::lu_backsubstitution(
1727  const Matrix& lu,
1728  const Index_vector& indx,
1729  Vector& b)
1730 {
1731  // when ii != DIM+1, ii is index of first non-vanishing element of b
1732  Size ii = DIM+1;
1733  for( Size i = 0; i < DIM; i++) {
1734  Size ip = indx[i];
1735  T sum = b[ip];
1736  b[ip] = b[i];
1737  if( ii != DIM+1) {
1738  for( Size j = ii; j < i; j++) {
1739  sum -= lu.get(i,j) * b[j];
1740  }
1741  } else {
1742  if( sum != T(0)) { // non-zero element
1743  ii = i;
1744  }
1745  }
1746  b[i] = sum;
1747  }
1748  for( Size i2 = DIM; i2 > 0;) { // backsubstitution
1749  --i2;
1750  T sum = b[i2];
1751  for( Size j = i2+1; j < DIM; j++)
1752  sum -= lu.get(i2,j) * b[j];
1753  b[i2] = sum / lu.get(i2,i2); // store comp. of sol. vector
1754  }
1755 }
1756 
1757 template <class T, Size DIM>
1758 bool Matrix_inverter<T,DIM,DIM>::invert( Matrix& mat)
1759 {
1760  Matrix lu(mat); // working copy of matrix to invert
1761  Index_vector indx; // row permutation info
1762 
1763  // LU decompose, return if it fails
1764  if( !lu_decomposition(lu, indx))
1765  return false;
1766 
1767  // solve for each column vector of the I matrix
1768  for( Size j = 0; j < DIM; ++j) {
1769  Vector col(T(0)); // TODO: do that directly in the mat matrix
1770  col[j] = T(1);
1771  lu_backsubstitution( lu, indx, col);
1772  for( Size i = 0; i < DIM; ++i) {
1773  mat.set( i, j, col[i]);
1774  }
1775  }
1776  return true;
1777 }
1778 
1779 // Matrix inversion class, specialization for 1x1 matrix
1780 template <class T>
1781 class Matrix_inverter<T,1,1>
1782 {
1783 public:
1784  typedef math::Matrix<T,1,1> Matrix;
1785 
1786  static inline bool invert( Matrix& mat)
1787  {
1788  T s = mat.get( 0, 0);
1789  if( s == T(0))
1790  return false;
1791  mat.set( 0, 0, T(1) / s);
1792  return true;
1793  }
1794 };
1795 
1796 // Matrix inversion class, specialization for 2x2 matrix
1797 template <class T>
1798 class Matrix_inverter<T,2,2>
1799 {
1800 public:
1801  typedef math::Matrix<T,2,2> Matrix;
1802 
1803  static inline bool invert( Matrix& mat)
1804  {
1805  T a = mat.get( 0, 0);
1806  T b = mat.get( 0, 1);
1807  T c = mat.get( 1, 0);
1808  T d = mat.get( 1, 1);
1809  T det = a*d - b*c;
1810  if( det == T( 0))
1811  return false;
1812  T rdet = T(1) / det;
1813  mat.set( 0, 0, d * rdet);
1814  mat.set( 0, 1,-b * rdet);
1815  mat.set( 1, 0,-c * rdet);
1816  mat.set( 1, 1, a * rdet);
1817  return true;
1818  }
1819 };
1820 
1821 template <typename T, Size ROW, Size COL>
1822 inline bool Matrix<T,ROW,COL>::invert()
1823 {
1824  return Matrix_inverter<T,ROW,COL>::invert( *this);
1825 }
1826 
1827 
1828 
1829 //------ Specializations and (specialized) overloads --------------------------
1830 
1831 // Specialization of common matrix multiplication for 4x4 matrices.
1832 template <typename T>
1833 inline Matrix<T,4,4>& operator*=(
1834  Matrix<T,4,4>& lhs,
1835  const Matrix<T,4,4>& rhs)
1836 {
1837  Matrix<T,4,4> old( lhs);
1838 
1839  lhs.xx = old.xx * rhs.xx + old.xy * rhs.yx + old.xz * rhs.zx + old.xw * rhs.wx;
1840  lhs.xy = old.xx * rhs.xy + old.xy * rhs.yy + old.xz * rhs.zy + old.xw * rhs.wy;
1841  lhs.xz = old.xx * rhs.xz + old.xy * rhs.yz + old.xz * rhs.zz + old.xw * rhs.wz;
1842  lhs.xw = old.xx * rhs.xw + old.xy * rhs.yw + old.xz * rhs.zw + old.xw * rhs.ww;
1843 
1844  lhs.yx = old.yx * rhs.xx + old.yy * rhs.yx + old.yz * rhs.zx + old.yw * rhs.wx;
1845  lhs.yy = old.yx * rhs.xy + old.yy * rhs.yy + old.yz * rhs.zy + old.yw * rhs.wy;
1846  lhs.yz = old.yx * rhs.xz + old.yy * rhs.yz + old.yz * rhs.zz + old.yw * rhs.wz;
1847  lhs.yw = old.yx * rhs.xw + old.yy * rhs.yw + old.yz * rhs.zw + old.yw * rhs.ww;
1848 
1849  lhs.zx = old.zx * rhs.xx + old.zy * rhs.yx + old.zz * rhs.zx + old.zw * rhs.wx;
1850  lhs.zy = old.zx * rhs.xy + old.zy * rhs.yy + old.zz * rhs.zy + old.zw * rhs.wy;
1851  lhs.zz = old.zx * rhs.xz + old.zy * rhs.yz + old.zz * rhs.zz + old.zw * rhs.wz;
1852  lhs.zw = old.zx * rhs.xw + old.zy * rhs.yw + old.zz * rhs.zw + old.zw * rhs.ww;
1853 
1854  lhs.wx = old.wx * rhs.xx + old.wy * rhs.yx + old.wz * rhs.zx + old.ww * rhs.wx;
1855  lhs.wy = old.wx * rhs.xy + old.wy * rhs.yy + old.wz * rhs.zy + old.ww * rhs.wy;
1856  lhs.wz = old.wx * rhs.xz + old.wy * rhs.yz + old.wz * rhs.zz + old.ww * rhs.wz;
1857  lhs.ww = old.wx * rhs.xw + old.wy * rhs.yw + old.wz * rhs.zw + old.ww * rhs.ww;
1858 
1859  return lhs;
1860 }
1861 
1862 // Specialization of common matrix multiplication for 4x4 matrices.
1863 template <typename T>
1864 inline Matrix<T,4,4> operator*(
1865  const Matrix<T,4,4>& lhs,
1866  const Matrix<T,4,4>& rhs)
1867 {
1868  Matrix<T,4,4> temp( lhs);
1869  temp *= rhs;
1870  return temp;
1871 }
1872 
1873 #endif // MI_FOR_DOXYGEN_ONLY
1874  // end group mi_math_matrix
1876 
1877 } // namespace math
1878 
1879 } // namespace mi
1880 
1881 #endif // MI_MATH_MATRIX_H