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 2019 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,4>& mat,
1247  const Vector<U,3>& point)
1248 {
1249  T w = T(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww);
1250  if( w == T(0) || w == T(1)) // avoids temporary and division
1251  return Vector<U,3>(
1252  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx),
1253  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy),
1254  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz));
1255  else {
1256  w = T(1)/w; // optimization
1257  return Vector<U,3>(
1258  U((mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx) * w),
1259  U((mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy) * w),
1260  U((mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz) * w));
1261  }
1262 }
1263 
1269 template <typename T, typename U>
1271  const Matrix<T,4,4>& mat,
1272  const Vector<U,4>& point)
1273 {
1274  return Vector<U, 4>(
1275  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx * point.w),
1276  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy * point.w),
1277  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz * point.w),
1278  U(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww * point.w));
1279 }
1280 
1286 template <typename T, typename U>
1288  const Matrix<T,4,4>& mat,
1289  const U& vector)
1290 {
1291  return U(mat.xx * vector);
1292 }
1293 
1299 template <typename T, typename U>
1301  const Matrix<T,4,4>& mat,
1302  const Vector<U,2>& vector)
1303 {
1304  return Vector<U,2>(
1305  U(mat.xx * vector.x + mat.yx * vector.y),
1306  U(mat.xy * vector.x + mat.yy * vector.y));
1307 }
1308 
1314 template <typename T, typename U>
1316  const Matrix<T,4,4>& mat,
1317  const Vector<U,3>& vector)
1318 {
1319  return Vector<U,3>(
1320  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1321  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1322  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1323 }
1324 
1336 template <typename T, typename U>
1338  const Matrix<T,4,4>& inv_mat,
1339  const Vector<U,3>& normal)
1340 {
1341  return Vector<U,3>(
1342  U(inv_mat.xx * normal.x + inv_mat.xy * normal.y + inv_mat.xz * normal.z),
1343  U(inv_mat.yx * normal.x + inv_mat.yy * normal.y + inv_mat.yz * normal.z),
1344  U(inv_mat.zx * normal.x + inv_mat.zy * normal.y + inv_mat.zz * normal.z));
1345 }
1346 
1360 template <typename T, typename U>
1362  const Matrix<T,4,4>& mat,
1363  const Vector<U,3>& normal)
1364 {
1365  Matrix<T,3,3> sub_mat( mat.xx, mat.xy, mat.xz,
1366  mat.yx, mat.yy, mat.yz,
1367  mat.zx, mat.zy, mat.zz);
1368  bool inverted = sub_mat.invert();
1369  if( !inverted)
1370  return normal;
1371  return Vector<U,3>(
1372  U(sub_mat.xx * normal.x + sub_mat.xy * normal.y + sub_mat.xz * normal.z),
1373  U(sub_mat.yx * normal.x + sub_mat.yy * normal.y + sub_mat.yz * normal.z),
1374  U(sub_mat.zx * normal.x + sub_mat.zy * normal.y + sub_mat.zz * normal.z));
1375 }
1376 
1377 
1378 //------ Definitions of non-inline function -----------------------------------
1379 
1380 template <typename T, Size ROW, Size COL>
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 template <typename T, Size ROW, Size COL>
1392  Matrix<T,ROW,COL>& lhs,
1393  const Matrix<T,ROW,COL>& rhs)
1394 {
1395  for( Size i=0; i < ROW*COL; ++i)
1396  lhs.begin()[i] -= rhs.get(i);
1397  return lhs;
1398 }
1399 
1400 #ifndef MI_FOR_DOXYGEN_ONLY
1401 
1402 template <typename T, Size ROW, Size COL>
1403 inline void Matrix<T,ROW,COL>::set_rotation( T xangle, T yangle, T zangle)
1404 {
1405  mi_static_assert( COL == 4 && ROW == 4);
1406  T tsx, tsy, tsz; // sine of [xyz]_angle
1407  T tcx, tcy, tcz; // cosine of [xyz]_angle
1408  T tmp;
1409  const T min_angle = T(0.00024f);
1410 
1411  if( abs( xangle) > min_angle) {
1412  tsx = sin( xangle);
1413  tcx = cos( xangle);
1414  } else {
1415  tsx = xangle;
1416  tcx = T(1);
1417  }
1418  if( abs( yangle) > min_angle) {
1419  tsy = sin( yangle);
1420  tcy = cos( yangle);
1421  } else {
1422  tsy = yangle;
1423  tcy = T(1);
1424  }
1425  if( abs(zangle) > min_angle) {
1426  tsz = sin( zangle);
1427  tcz = cos( zangle);
1428  } else {
1429  tsz = T(zangle);
1430  tcz = T(1);
1431  }
1432  this->xx = tcy * tcz;
1433  this->xy = tcy * tsz;
1434  this->xz = -tsy;
1435 
1436  tmp = tsx * tsy;
1437  this->yx = tmp * tcz - tcx * tsz;
1438  this->yy = tmp * tsz + tcx * tcz;
1439  this->yz = tsx * tcy;
1440 
1441  tmp = tcx * tsy;
1442  this->zx = tmp * tcz + tsx * tsz;
1443  this->zy = tmp * tsz - tsx * tcz;
1444  this->zz = tcx * tcy;
1445 }
1446 
1447 template <typename T, Size ROW, Size COL>
1448 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float32,3>& axis_v, Float64 angle)
1449 {
1450  mi_static_assert( COL == 4 && ROW == 4);
1451  Vector<T,3> axis( axis_v);
1452  const T min_angle = T(0.00024f);
1453 
1454  if( abs( T(angle)) < min_angle) {
1455  T xa = axis.x * T(angle);
1456  T ya = axis.y * T(angle);
1457  T za = axis.z * T(angle);
1458 
1459  this->xx = T(1);
1460  this->xy = za;
1461  this->xz = -ya;
1462  this->xw = T(0);
1463 
1464  this->yx = za;
1465  this->yy = T(1);
1466  this->yz = xa;
1467  this->yw = T(0);
1468 
1469  this->zx = -ya;
1470  this->zy = -xa;
1471  this->zz = T(1);
1472  this->zw = T(0);
1473  } else {
1474  T s = sin( T(angle));
1475  T c = cos( T(angle));
1476  T t = T(1) - c;
1477  T tmp;
1478 
1479  tmp = t * T(axis.x);
1480  this->xx = tmp * T(axis.x) + c;
1481  this->xy = tmp * T(axis.y) + s * T(axis.z);
1482  this->xz = tmp * T(axis.z) - s * T(axis.y);
1483  this->xw = T(0);
1484 
1485  tmp = t * T(axis.y);
1486  this->yx = tmp * T(axis.x) - s * T(axis.z);
1487  this->yy = tmp * T(axis.y) + c;
1488  this->yz = tmp * T(axis.z) + s * T(axis.x);
1489  this->yw = T(0);
1490 
1491  tmp = t * T(axis.z);
1492  this->zx = tmp * T(axis.x) + s * T(axis.y);
1493  this->zy = tmp * T(axis.y) - s * T(axis.x);
1494  this->zz = tmp * T(axis.z) + c;
1495  this->zw = T(0);
1496  }
1497  this->wx = this->wy = this->wz = T(0);
1498  this->ww = T(1);
1499 }
1500 
1501 template <typename T, Size ROW, Size COL>
1502 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float64,3>& axis_v, Float64 angle)
1503 {
1504  mi_static_assert( COL == 4 && ROW == 4);
1505  Vector<T,3> axis( axis_v);
1506  const T min_angle = T(0.00024f);
1507 
1508  if( abs(T(angle)) < min_angle) {
1509  T xa = axis.x * T(angle);
1510  T ya = axis.y * T(angle);
1511  T za = axis.z * T(angle);
1512 
1513  this->xx = T(1);
1514  this->xy = za;
1515  this->xz = -ya;
1516  this->xw = T(0);
1517 
1518  this->yx = za;
1519  this->yy = T(1);
1520  this->yz = xa;
1521  this->yw = T(0);
1522 
1523  this->zx = -ya;
1524  this->zy = -xa;
1525  this->zz = T(1);
1526  this->zw = T(0);
1527  } else {
1528  T s = sin( T(angle));
1529  T c = cos( T(angle));
1530  T t = T(1) - c;
1531  T tmp;
1532 
1533  tmp = t * T(axis.x);
1534  this->xx = tmp * T(axis.x) + c;
1535  this->xy = tmp * T(axis.y) + s * T(axis.z);
1536  this->xz = tmp * T(axis.z) - s * T(axis.y);
1537  this->xw = T(0);
1538 
1539  tmp = t * T(axis.y);
1540  this->yx = tmp * T(axis.x) - s * T(axis.z);
1541  this->yy = tmp * T(axis.y) + c;
1542  this->yz = tmp * T(axis.z) + s * T(axis.x);
1543  this->yw = T(0);
1544 
1545  tmp = t * T(axis.z);
1546  this->zx = tmp * T(axis.x) + s * T(axis.y);
1547  this->zy = tmp * T(axis.y) - s * T(axis.x);
1548  this->zz = tmp * T(axis.z) + c;
1549  this->zw = T(0);
1550  }
1551  this->wx = this->wy = this->wz = T(0);
1552  this->ww = T(1);
1553 }
1554 
1555 template <typename T, Size ROW, Size COL>
1556 inline void Matrix<T,ROW,COL>::lookat(
1557  const Vector<Float32,3>& position,
1558  const Vector<Float32,3>& target,
1559  const Vector<Float32,3>& up)
1560 {
1561  mi_static_assert( COL == 4 && ROW == 4);
1562  Vector<Float32,3> xaxis, yaxis, zaxis;
1563 
1564  // Z vector
1565  zaxis = position - target;
1566  zaxis.normalize();
1567 
1568  // X vector = up cross Z
1569  xaxis = cross( up, zaxis);
1570  xaxis.normalize();
1571 
1572  // Recompute Y = Z cross X
1573  yaxis = cross( zaxis, xaxis);
1574  yaxis.normalize();
1575 
1576  // Build rotation matrix
1577  Matrix<T,4,4> rot(
1578  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1579  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1580  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1581  T(0), T(0), T(0), T(1));
1582 
1583  // Compute the new position by multiplying the inverse position with the rotation matrix
1584  Matrix<T,4,4> trans(
1585  T(1), T(0), T(0), T(0),
1586  T(0), T(1), T(0), T(0),
1587  T(0), T(0), T(1), T(0),
1588  T(-position.x), T(-position.y), T(-position.z), T(1));
1589 
1590  *this = trans * rot;
1591 }
1592 
1593 template <typename T, Size ROW, Size COL>
1594 inline void Matrix<T,ROW,COL>::lookat(
1595  const Vector<Float64,3>& position,
1596  const Vector<Float64,3>& target,
1597  const Vector<Float64,3>& up)
1598 {
1599  mi_static_assert( COL == 4 && ROW == 4);
1600  Vector<Float64,3> xaxis, yaxis, zaxis;
1601 
1602  // Z vector
1603  zaxis = position - target;
1604  zaxis.normalize();
1605 
1606  // X vector = up cross Z
1607  xaxis = cross( up, zaxis);
1608  xaxis.normalize();
1609 
1610  // Recompute Y = Z cross X
1611  yaxis = cross( zaxis, xaxis);
1612  yaxis.normalize();
1613 
1614  // Build rotation matrix
1615  Matrix<T,4,4> rot(
1616  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1617  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1618  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1619  T(0), T(0), T(0), T(1));
1620 
1621  // Compute the new position by multiplying the inverse position with the rotation matrix
1622  Matrix<T,4,4> trans(
1623  T(1), T(0), T(0), T(0),
1624  T(0), T(1), T(0), T(0),
1625  T(0), T(0), T(1), T(0),
1626  T(-position.x), T(-position.y), T(-position.z), T(1));
1627 
1628  *this = trans * rot;
1629 }
1630 
1631 
1632 //------ Definition and helper for matrix inversion ---------------------------
1633 
1634 // Matrix inversion class, specialized for symmetric matrices and low dimensions
1635 template <class T, Size ROW, Size COL>
1636 class Matrix_inverter
1637 {
1638 public:
1639  typedef math::Matrix<T,ROW,COL> Matrix;
1640 
1641  // Inverts the matrix \c M if possible and returns \c true.
1642  //
1643  // If the matrix cannot be inverted or if \c ROW != \c COL, this function returns \c false.
1644  static inline bool invert( Matrix& /* M */) { return false; }
1645 };
1646 
1647 // Matrix inversion class, specialization for symmetric matrices
1648 template <class T, Size DIM>
1649 class Matrix_inverter<T,DIM,DIM>
1650 {
1651 public:
1652  typedef math::Matrix<T,DIM,DIM> Matrix;
1653  typedef math::Vector<T,DIM> Vector;
1654  typedef math::Vector<Size,DIM> Index_vector;
1655 
1656  // LU decomposition of matrix lu in place.
1657  //
1658  // Returns \c false if matrix cannot be decomposed, for example, if it is not invertible, and \c
1659  // true otherwise. Returns also a row permutation in indx.
1660  static bool lu_decomposition(
1661  Matrix& lu, // matrix to decompose, modified in place
1662  Index_vector& indx); // row permutation info indx[4]
1663 
1664  // Solves the equation lu * x = b for x. lu and indx are the results of lu_decomposition. The
1665  // solution is returned in b.
1666  static void lu_backsubstitution(
1667  const Matrix& lu, // LU decomposed matrix
1668  const Index_vector& indx, // permutation vector
1669  Vector& b); // right hand side vector b, modified in place
1670 
1671  static bool invert( Matrix& mat);
1672 };
1673 
1674 template <class T, Size DIM>
1675 bool Matrix_inverter<T,DIM,DIM>::lu_decomposition(
1676  Matrix& lu,
1677  Index_vector& indx)
1678 {
1679  Vector vv; // implicit scaling of each row
1680 
1681  for( Size i = 0; i < DIM; i++) { // get implicit scaling
1682  T big = T(0);
1683  for( Size j = 0; j < DIM; j++) {
1684  T temp = abs(lu.get(i,j));
1685  if( temp > big)
1686  big = temp;
1687  }
1688  if( big == T(0))
1689  return false;
1690  vv[i] = T(1) / big; // save scaling
1691  }
1692  //
1693  // loop over columns of Crout's method
1694  //
1695  Size imax = 0;
1696  for( Size j = 0; j < DIM; j++) {
1697  for( Size i = 0; i < j; i++) {
1698  T sum = lu.get(i,j);
1699  for( Size k = 0; k < i; k++)
1700  sum -= lu.get(i,k) * lu.get(k,j);
1701  lu.set(i, j, sum);
1702  }
1703  T big = 0; // init for pivot search
1704  for( Size i = j; i < DIM; i++) {
1705  T sum = lu.get(i,j);
1706  for( Size k = 0; k < j; k++)
1707  sum -= lu.get(i,k) * lu.get(k,j);
1708  lu.set(i, j, sum);
1709  T dum = vv[i] * abs(sum);
1710  if( dum >= big) { // pivot good?
1711  big = dum;
1712  imax = i;
1713  }
1714  }
1715  if( j != imax) { // interchange rows
1716  for( Size k = 0; k < DIM; k++) {
1717  T dum = lu.get(imax,k);
1718  lu.set(imax, k, lu.get(j,k));
1719  lu.set(j, k, dum);
1720  }
1721  vv[imax] = vv[j]; // interchange scale factor
1722  }
1723  indx[j] = imax;
1724  if( lu.get(j,j) == 0)
1725  return false;
1726  if( j != DIM-1) { // divide by pivot element
1727  T dum = T(1) / lu.get(j,j);
1728  for( Size i = j + 1; i < DIM; i++)
1729  lu.set(i, j, lu.get(i,j) * dum);
1730  }
1731  }
1732  return true;
1733 }
1734 
1735 template <class T, Size DIM>
1736 void Matrix_inverter<T,DIM,DIM>::lu_backsubstitution(
1737  const Matrix& lu,
1738  const Index_vector& indx,
1739  Vector& b)
1740 {
1741  // when ii != DIM+1, ii is index of first non-vanishing element of b
1742  Size ii = DIM+1;
1743  for( Size i = 0; i < DIM; i++) {
1744  Size ip = indx[i];
1745  T sum = b[ip];
1746  b[ip] = b[i];
1747  if( ii != DIM+1) {
1748  for( Size j = ii; j < i; j++) {
1749  sum -= lu.get(i,j) * b[j];
1750  }
1751  } else {
1752  if( sum != T(0)) { // non-zero element
1753  ii = i;
1754  }
1755  }
1756  b[i] = sum;
1757  }
1758  for( Size i2 = DIM; i2 > 0;) { // backsubstitution
1759  --i2;
1760  T sum = b[i2];
1761  for( Size j = i2+1; j < DIM; j++)
1762  sum -= lu.get(i2,j) * b[j];
1763  b[i2] = sum / lu.get(i2,i2); // store comp. of sol. vector
1764  }
1765 }
1766 
1767 template <class T, Size DIM>
1768 bool Matrix_inverter<T,DIM,DIM>::invert( Matrix& mat)
1769 {
1770  Matrix lu(mat); // working copy of matrix to invert
1771  Index_vector indx; // row permutation info
1772 
1773  // LU decompose, return if it fails
1774  if( !lu_decomposition(lu, indx))
1775  return false;
1776 
1777  // solve for each column vector of the I matrix
1778  for( Size j = 0; j < DIM; ++j) {
1779  Vector col(T(0)); // TODO: do that directly in the mat matrix
1780  col[j] = T(1);
1781  lu_backsubstitution( lu, indx, col);
1782  for( Size i = 0; i < DIM; ++i) {
1783  mat.set( i, j, col[i]);
1784  }
1785  }
1786  return true;
1787 }
1788 
1789 // Matrix inversion class, specialization for 1x1 matrix
1790 template <class T>
1791 class Matrix_inverter<T,1,1>
1792 {
1793 public:
1794  typedef math::Matrix<T,1,1> Matrix;
1795 
1796  static inline bool invert( Matrix& mat)
1797  {
1798  T s = mat.get( 0, 0);
1799  if( s == T(0))
1800  return false;
1801  mat.set( 0, 0, T(1) / s);
1802  return true;
1803  }
1804 };
1805 
1806 // Matrix inversion class, specialization for 2x2 matrix
1807 template <class T>
1808 class Matrix_inverter<T,2,2>
1809 {
1810 public:
1811  typedef math::Matrix<T,2,2> Matrix;
1812 
1813  static inline bool invert( Matrix& mat)
1814  {
1815  T a = mat.get( 0, 0);
1816  T b = mat.get( 0, 1);
1817  T c = mat.get( 1, 0);
1818  T d = mat.get( 1, 1);
1819  T det = a*d - b*c;
1820  if( det == T( 0))
1821  return false;
1822  T rdet = T(1) / det;
1823  mat.set( 0, 0, d * rdet);
1824  mat.set( 0, 1,-b * rdet);
1825  mat.set( 1, 0,-c * rdet);
1826  mat.set( 1, 1, a * rdet);
1827  return true;
1828  }
1829 };
1830 
1831 template <typename T, Size ROW, Size COL>
1832 inline bool Matrix<T,ROW,COL>::invert()
1833 {
1834  return Matrix_inverter<T,ROW,COL>::invert( *this);
1835 }
1836 
1837 
1838 
1839 //------ Specializations and (specialized) overloads --------------------------
1840 
1841 // Specialization of common matrix multiplication for 4x4 matrices.
1842 template <typename T>
1843 inline Matrix<T,4,4>& operator*=(
1844  Matrix<T,4,4>& lhs,
1845  const Matrix<T,4,4>& rhs)
1846 {
1847  Matrix<T,4,4> old( lhs);
1848 
1849  lhs.xx = old.xx * rhs.xx + old.xy * rhs.yx + old.xz * rhs.zx + old.xw * rhs.wx;
1850  lhs.xy = old.xx * rhs.xy + old.xy * rhs.yy + old.xz * rhs.zy + old.xw * rhs.wy;
1851  lhs.xz = old.xx * rhs.xz + old.xy * rhs.yz + old.xz * rhs.zz + old.xw * rhs.wz;
1852  lhs.xw = old.xx * rhs.xw + old.xy * rhs.yw + old.xz * rhs.zw + old.xw * rhs.ww;
1853 
1854  lhs.yx = old.yx * rhs.xx + old.yy * rhs.yx + old.yz * rhs.zx + old.yw * rhs.wx;
1855  lhs.yy = old.yx * rhs.xy + old.yy * rhs.yy + old.yz * rhs.zy + old.yw * rhs.wy;
1856  lhs.yz = old.yx * rhs.xz + old.yy * rhs.yz + old.yz * rhs.zz + old.yw * rhs.wz;
1857  lhs.yw = old.yx * rhs.xw + old.yy * rhs.yw + old.yz * rhs.zw + old.yw * rhs.ww;
1858 
1859  lhs.zx = old.zx * rhs.xx + old.zy * rhs.yx + old.zz * rhs.zx + old.zw * rhs.wx;
1860  lhs.zy = old.zx * rhs.xy + old.zy * rhs.yy + old.zz * rhs.zy + old.zw * rhs.wy;
1861  lhs.zz = old.zx * rhs.xz + old.zy * rhs.yz + old.zz * rhs.zz + old.zw * rhs.wz;
1862  lhs.zw = old.zx * rhs.xw + old.zy * rhs.yw + old.zz * rhs.zw + old.zw * rhs.ww;
1863 
1864  lhs.wx = old.wx * rhs.xx + old.wy * rhs.yx + old.wz * rhs.zx + old.ww * rhs.wx;
1865  lhs.wy = old.wx * rhs.xy + old.wy * rhs.yy + old.wz * rhs.zy + old.ww * rhs.wy;
1866  lhs.wz = old.wx * rhs.xz + old.wy * rhs.yz + old.wz * rhs.zz + old.ww * rhs.wz;
1867  lhs.ww = old.wx * rhs.xw + old.wy * rhs.yw + old.wz * rhs.zw + old.ww * rhs.ww;
1868 
1869  return lhs;
1870 }
1871 
1872 // Specialization of common matrix multiplication for 4x4 matrices.
1873 template <typename T>
1874 inline Matrix<T,4,4> operator*(
1875  const Matrix<T,4,4>& lhs,
1876  const Matrix<T,4,4>& rhs)
1877 {
1878  Matrix<T,4,4> temp( lhs);
1879  temp *= rhs;
1880  return temp;
1881 }
1882 
1883 #endif // MI_FOR_DOXYGEN_ONLY
1884  // end group mi_math_matrix
1886 
1887 } // namespace math
1888 
1889 } // namespace mi
1890 
1891 #endif // MI_MATH_MATRIX_H