NVIDIA Iray: Math API
 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, 2014 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 #ifdef MI_NEURAYLIB_DEPRECATED_5_0
929  inline void rotateaxis( const Vector<Float32,3>& axis, Float64 angle)
930  {
931  set_rotation( axis, angle);
932  }
933 
934  inline void rotateaxis( const Vector<Float64,3>& axis, Float64 angle)
935  {
936  set_rotation( axis, angle);
937  }
938 
939  inline void multiply( const Matrix<Float32,4,4>& matrix);
940 
941  inline void multiply( const Matrix<Float64,4,4>& matrix);
942 
943  inline Vector<Float32,3> transform( const Vector<Float32,3> vector) const;
944 
945  inline Vector<Float64,3> transform( const Vector<Float64,3> vector) const;
946 
947  inline Vector<Float32,3> transform33( const Vector<Float32,3> vector) const;
948 
949  inline Vector<Float64,3> transform33( const Vector<Float64,3> vector) const;
950 
951  inline Vector<Float32,3> transform33t( const Vector<Float32,3> vector) const;
952 
953  inline Vector<Float64,3> transform33t( const Vector<Float64,3> vector) const;
954 #endif // MI_NEURAYLIB_DEPRECATED_5_0
955 };
956 
957 
958 //------ Free comparison operators ==, !=, <, <=, >, >= for matrices -----------
959 
961 template <typename T, Size ROW, Size COL>
962 inline bool operator==(
963  const Matrix<T,ROW,COL>& lhs,
964  const Matrix<T,ROW,COL>& rhs)
965 {
966  return EQUAL == lexicographically_compare( lhs, rhs);
967 }
968 
970 template <typename T, Size ROW, Size COL>
971 inline bool operator!=(
972  const Matrix<T,ROW,COL>& lhs,
973  const Matrix<T,ROW,COL>& rhs)
974 {
975  return EQUAL != lexicographically_compare( lhs, rhs);
976 }
977 
981 template <typename T, Size ROW, Size COL>
982 inline bool operator<(
983  const Matrix<T,ROW,COL>& lhs,
984  const Matrix<T,ROW,COL>& rhs)
985 {
986  return lexicographically_less( lhs, rhs);
987 }
988 
992 template <typename T, Size ROW, Size COL>
993 inline bool operator<=(
994  const Matrix<T,ROW,COL>& lhs,
995  const Matrix<T,ROW,COL>& rhs)
996 {
997  return ! lexicographically_less( rhs, lhs);
998 }
999 
1003 template <typename T, Size ROW, Size COL>
1004 inline bool operator>(
1005  const Matrix<T,ROW,COL>& lhs,
1006  const Matrix<T,ROW,COL>& rhs)
1007 {
1008  return lexicographically_less( rhs, lhs);
1009 }
1010 
1014 template <typename T, Size ROW, Size COL>
1015 inline bool operator>=(
1016  const Matrix<T,ROW,COL>& lhs,
1017  const Matrix<T,ROW,COL>& rhs)
1018 {
1019  return ! lexicographically_less( lhs, rhs);
1020 }
1021 
1022 //------ Operator declarations for Matrix -------------------------------------
1023 
1025 template <typename T, Size ROW, Size COL>
1026 Matrix<T,ROW,COL>& operator+=( Matrix<T,ROW,COL>& lhs, const Matrix<T,ROW,COL>& rhs);
1027 
1029 template <typename T, Size ROW, Size COL>
1030 Matrix<T,ROW,COL>& operator-=( Matrix<T,ROW,COL>& lhs, const Matrix<T,ROW,COL>& rhs);
1031 
1033 template <typename T, Size ROW, Size COL>
1034 inline Matrix<T,ROW,COL> operator+(
1035  const Matrix<T,ROW,COL>& lhs,
1036  const Matrix<T,ROW,COL>& rhs)
1037 {
1038  Matrix<T,ROW,COL> temp( lhs);
1039  temp += rhs;
1040  return temp;
1041 }
1042 
1044 template <typename T, Size ROW, Size COL>
1045 inline Matrix<T,ROW,COL> operator-(
1046  const Matrix<T,ROW,COL>& lhs,
1047  const Matrix<T,ROW,COL>& rhs)
1048 {
1049  Matrix<T,ROW,COL> temp( lhs);
1050  temp -= rhs;
1051  return temp;
1052 }
1053 
1055 template <typename T, Size ROW, Size COL>
1056 inline Matrix<T,ROW,COL> operator-(
1057  const Matrix<T,ROW,COL>& mat)
1058 {
1059  Matrix<T,ROW,COL> temp;
1060  for( Size i(0u); i < (ROW*COL); ++i)
1061  temp.begin()[i] = -mat.get(i);
1062  return temp;
1063 }
1064 
1071 template <typename T, Size ROW, Size COL>
1072 inline Matrix<T,ROW,COL>& operator*=(
1073  Matrix<T,ROW,COL>& lhs,
1074  const Matrix<T,COL,COL>& rhs)
1075 {
1076  // There are more efficient ways of implementing this. Its a default solution. Specializations
1077  // exist below.
1078  Matrix<T,ROW,COL> old( lhs);
1079  for( Size rrow = 0; rrow < ROW; ++rrow) {
1080  for( Size rcol = 0; rcol < COL; ++rcol) {
1081  lhs( rrow, rcol) = T(0);
1082  for( Size k = 0; k < COL; ++k)
1083  lhs( rrow, rcol) += old( rrow, k) * rhs( k, rcol);
1084  }
1085  }
1086  return lhs;
1087 }
1088 
1094 template <typename T, Size ROW1, Size COL1, Size ROW2, Size COL2>
1095 inline Matrix<T,ROW1,COL2> operator*(
1096  const Matrix<T,ROW1,COL1>& lhs,
1097  const Matrix<T,ROW2,COL2>& rhs)
1098 {
1099  // There are more efficient ways of implementing this. Its a default solution. Specializations
1100  // exist below.
1101  mi_static_assert( COL1 == ROW2);
1102  Matrix<T,ROW1,COL2> result;
1103  for( Size rrow = 0; rrow < ROW1; ++rrow) {
1104  for( Size rcol = 0; rcol < COL2; ++rcol) {
1105  result( rrow, rcol) = T(0);
1106  for( Size k = 0; k < COL1; ++k)
1107  result( rrow, rcol) += lhs( rrow, k) * rhs( k, rcol);
1108  }
1109  }
1110  return result;
1111 }
1112 
1120 template <typename T, Size ROW, Size COL, Size DIM>
1121 inline Vector<T,ROW> operator*(
1122  const Matrix<T,ROW,COL>& mat,
1123  const Vector<T,DIM>& vec)
1124 {
1125  mi_static_assert( COL == DIM);
1126  Vector<T,ROW> result;
1127  for( Size row = 0; row < ROW; ++row) {
1128  result[row] = T(0);
1129  for( Size col = 0; col < COL; ++col)
1130  result[row] += mat( row, col) * vec[col];
1131  }
1132  return result;
1133 }
1134 
1141 template <Size DIM, typename T, Size ROW, Size COL>
1142 inline Vector<T,COL> operator*(
1143  const Vector<T,DIM>& vec,
1144  const Matrix<T,ROW,COL>& mat)
1145 {
1146  mi_static_assert( DIM == ROW);
1147  Vector<T,COL> result;
1148  for( Size col = 0; col < COL; ++col) {
1149  result[col] = T(0);
1150  for( Size row = 0; row < ROW; ++row)
1151  result[col] += mat( row, col) * vec[row];
1152  }
1153  return result;
1154 }
1155 
1158 template <typename T, Size ROW, Size COL>
1159 inline Matrix<T,ROW,COL>& operator*=( Matrix<T,ROW,COL>& mat, T factor)
1160 {
1161  for( Size i=0; i < (ROW*COL); ++i)
1162  mat.begin()[i] *= factor;
1163  return mat;
1164 }
1165 
1167 template <typename T, Size ROW, Size COL>
1168 inline Matrix<T,ROW,COL> operator*( const Matrix<T,ROW,COL>& mat, T factor)
1169 {
1170  Matrix<T,ROW,COL> temp( mat);
1171  temp *= factor;
1172  return temp;
1173 }
1174 
1176 template <typename T, Size ROW, Size COL>
1177 inline Matrix<T,ROW,COL> operator*( T factor, const Matrix<T,ROW,COL>& mat)
1178 {
1179  Matrix<T,ROW,COL> temp( mat);
1180  temp *= factor; // * on T should be commutative (IEEE-754)
1181  return temp;
1182 }
1183 
1184 
1185 //------ Free Functions for Matrix --------------------------------------------
1186 
1190 template <Size NEW_ROW, Size NEW_COL, typename T, Size ROW, Size COL>
1191 inline Matrix<T,NEW_ROW,NEW_COL> sub_matrix(
1192  const Matrix<T,ROW,COL>& mat)
1193 {
1194  mi_static_assert( NEW_ROW <= ROW);
1195  mi_static_assert( NEW_COL <= COL);
1196  Matrix<T,NEW_ROW,NEW_COL> result;
1197  for( Size i=0; i < NEW_ROW; ++i)
1198  for( Size j=0; j < NEW_COL; ++j)
1199  result( i, j) = mat( i, j);
1200  return result;
1201 }
1202 
1203 
1204 
1206 template <typename T, Size ROW, Size COL>
1207 inline Matrix<T,COL,ROW> transpose(
1208  const Matrix<T,ROW,COL>& mat)
1209 {
1210  Matrix<T,COL,ROW> result;
1211  for( Size i=0; i < ROW; ++i)
1212  for( Size j=0; j < COL; ++j)
1213  result( j, i) = mat( i, j);
1214  return result;
1215 }
1216 
1222 template <typename T, typename U>
1223 inline U transform_point(
1224  const Matrix<T,4,4>& mat,
1225  const U& point)
1226 {
1227  const T w = T(mat.xw * point + mat.ww);
1228  if( w == T(0) || w == T(1))
1229  return U(mat.xx * point + mat.wx);
1230  else
1231  return U((mat.xx * point + mat.wx) / w);
1232 }
1233 
1239 template <typename T, typename U>
1240 inline Vector<U,2> transform_point(
1241  const Matrix<T,4,4>& mat,
1242  const Vector<U,2>& point)
1243 {
1244  T w = T(mat.xw * point.x + mat.yw * point.y + mat.ww);
1245  if( w == T(0) || w == T(1))
1246  return Vector<U, 2>(
1247  U(mat.xx * point.x + mat.yx * point.y + mat.wx),
1248  U(mat.xy * point.x + mat.yy * point.y + mat.wy));
1249  else {
1250  w = T(1)/w; // optimization
1251  return Vector<U, 2>(
1252  U((mat.xx * point.x + mat.yx * point.y + mat.wx) * w),
1253  U((mat.xy * point.x + mat.yy * point.y + mat.wy) * w));
1254  }
1255 }
1256 
1262 template <typename T, typename U>
1263 inline Vector<U,3> transform_point(
1264  const Matrix<T,4,4>& mat,
1265  const Vector<U,3>& point)
1266 {
1267  T w = T(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww);
1268  if( w == T(0) || w == T(1)) // avoids temporary and division
1269  return Vector<U,3>(
1270  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx),
1271  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy),
1272  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz));
1273  else {
1274  w = T(1)/w; // optimization
1275  return Vector<U,3>(
1276  U((mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx) * w),
1277  U((mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy) * w),
1278  U((mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz) * w));
1279  }
1280 }
1281 
1287 template <typename T, typename U>
1288 inline Vector<U,4> transform_point(
1289  const Matrix<T,4,4>& mat,
1290  const Vector<U,4>& point)
1291 {
1292  return Vector<U, 4>(
1293  U(mat.xx * point.x + mat.yx * point.y + mat.zx * point.z + mat.wx * point.w),
1294  U(mat.xy * point.x + mat.yy * point.y + mat.zy * point.z + mat.wy * point.w),
1295  U(mat.xz * point.x + mat.yz * point.y + mat.zz * point.z + mat.wz * point.w),
1296  U(mat.xw * point.x + mat.yw * point.y + mat.zw * point.z + mat.ww * point.w));
1297 }
1298 
1304 template <typename T, typename U>
1305 inline U transform_vector(
1306  const Matrix<T,4,4>& mat,
1307  const U& vector)
1308 {
1309  return U(mat.xx * vector);
1310 }
1311 
1317 template <typename T, typename U>
1318 inline Vector<U,2> transform_vector(
1319  const Matrix<T,4,4>& mat,
1320  const Vector<U,2>& vector)
1321 {
1322  return Vector<U,2>(
1323  U(mat.xx * vector.x + mat.yx * vector.y),
1324  U(mat.xy * vector.x + mat.yy * vector.y));
1325 }
1326 
1332 template <typename T, typename U>
1333 inline Vector<U,3> transform_vector(
1334  const Matrix<T,4,4>& mat,
1335  const Vector<U,3>& vector)
1336 {
1337  return Vector<U,3>(
1338  U(mat.xx * vector.x + mat.yx * vector.y + mat.zx * vector.z),
1339  U(mat.xy * vector.x + mat.yy * vector.y + mat.zy * vector.z),
1340  U(mat.xz * vector.x + mat.yz * vector.y + mat.zz * vector.z));
1341 }
1342 
1354 template <typename T, typename U>
1355 inline Vector<U,3> transform_normal_inv(
1356  const Matrix<T,4,4>& inv_mat,
1357  const Vector<U,3>& normal)
1358 {
1359  return Vector<U,3>(
1360  U(inv_mat.xx * normal.x + inv_mat.xy * normal.y + inv_mat.xz * normal.z),
1361  U(inv_mat.yx * normal.x + inv_mat.yy * normal.y + inv_mat.yz * normal.z),
1362  U(inv_mat.zx * normal.x + inv_mat.zy * normal.y + inv_mat.zz * normal.z));
1363 }
1364 
1378 template <typename T, typename U>
1379 inline Vector<U,3> transform_normal(
1380  const Matrix<T,4,4>& mat,
1381  const Vector<U,3>& normal)
1382 {
1383  Matrix<T,3,3> sub_mat( mat.xx, mat.xy, mat.xz,
1384  mat.yx, mat.yy, mat.yz,
1385  mat.zx, mat.zy, mat.zz);
1386  bool inverted = sub_mat.invert();
1387  if( !inverted)
1388  return normal;
1389  return Vector<U,3>(
1390  U(sub_mat.xx * normal.x + sub_mat.xy * normal.y + sub_mat.xz * normal.z),
1391  U(sub_mat.yx * normal.x + sub_mat.yy * normal.y + sub_mat.yz * normal.z),
1392  U(sub_mat.zx * normal.x + sub_mat.zy * normal.y + sub_mat.zz * normal.z));
1393 }
1394 
1395 
1396 //------ Definitions of non-inline function -----------------------------------
1397 
1398 template <typename T, Size ROW, Size COL>
1399 Matrix<T,ROW,COL>& operator+=(
1400  Matrix<T,ROW,COL>& lhs,
1401  const Matrix<T,ROW,COL>& rhs)
1402 {
1403  for( Size i=0; i < ROW*COL; ++i)
1404  lhs.begin()[i] += rhs.get(i);
1405  return lhs;
1406 }
1407 
1408 template <typename T, Size ROW, Size COL>
1409 Matrix<T,ROW,COL>& operator-=(
1410  Matrix<T,ROW,COL>& lhs,
1411  const Matrix<T,ROW,COL>& rhs)
1412 {
1413  for( Size i=0; i < ROW*COL; ++i)
1414  lhs.begin()[i] -= rhs.get(i);
1415  return lhs;
1416 }
1417 
1418 #ifndef MI_FOR_DOXYGEN_ONLY
1419 
1420 template <typename T, Size ROW, Size COL>
1421 inline void Matrix<T,ROW,COL>::set_rotation( T xangle, T yangle, T zangle)
1422 {
1423  mi_static_assert( COL == 4 && ROW == 4);
1424  T tsx, tsy, tsz; // sine of [xyz]_angle
1425  T tcx, tcy, tcz; // cosine of [xyz]_angle
1426  T tmp;
1427  const T min_angle = T(0.00024f);
1428 
1429  if( abs( xangle) > min_angle) {
1430  tsx = sin( xangle);
1431  tcx = cos( xangle);
1432  } else {
1433  tsx = xangle;
1434  tcx = T(1);
1435  }
1436  if( abs( yangle) > min_angle) {
1437  tsy = sin( yangle);
1438  tcy = cos( yangle);
1439  } else {
1440  tsy = yangle;
1441  tcy = T(1);
1442  }
1443  if( abs(zangle) > min_angle) {
1444  tsz = sin( zangle);
1445  tcz = cos( zangle);
1446  } else {
1447  tsz = T(zangle);
1448  tcz = T(1);
1449  }
1450  this->xx = tcy * tcz;
1451  this->xy = tcy * tsz;
1452  this->xz = -tsy;
1453 
1454  tmp = tsx * tsy;
1455  this->yx = tmp * tcz - tcx * tsz;
1456  this->yy = tmp * tsz + tcx * tcz;
1457  this->yz = tsx * tcy;
1458 
1459  tmp = tcx * tsy;
1460  this->zx = tmp * tcz + tsx * tsz;
1461  this->zy = tmp * tsz - tsx * tcz;
1462  this->zz = tcx * tcy;
1463 }
1464 
1465 template <typename T, Size ROW, Size COL>
1466 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float32,3>& axis_v, Float64 angle)
1467 {
1468  mi_static_assert( COL == 4 && ROW == 4);
1469  Vector<T,3> axis( axis_v);
1470  const T min_angle = T(0.00024f);
1471 
1472  if( abs( T(angle)) < min_angle) {
1473  T xa = axis.x * T(angle);
1474  T ya = axis.y * T(angle);
1475  T za = axis.z * T(angle);
1476 
1477  this->xx = T(1);
1478  this->xy = za;
1479  this->xz = -ya;
1480  this->xw = T(0);
1481 
1482  this->yx = za;
1483  this->yy = T(1);
1484  this->yz = xa;
1485  this->yw = T(0);
1486 
1487  this->zx = -ya;
1488  this->zy = -xa;
1489  this->zz = T(1);
1490  this->zw = T(0);
1491  } else {
1492  T s = sin( T(angle));
1493  T c = cos( T(angle));
1494  T t = T(1) - c;
1495  T tmp;
1496 
1497  tmp = t * T(axis.x);
1498  this->xx = tmp * T(axis.x) + c;
1499  this->xy = tmp * T(axis.y) + s * T(axis.z);
1500  this->xz = tmp * T(axis.z) - s * T(axis.y);
1501  this->xw = T(0);
1502 
1503  tmp = t * T(axis.y);
1504  this->yx = tmp * T(axis.x) - s * T(axis.z);
1505  this->yy = tmp * T(axis.y) + c;
1506  this->yz = tmp * T(axis.z) + s * T(axis.x);
1507  this->yw = T(0);
1508 
1509  tmp = t * T(axis.z);
1510  this->zx = tmp * T(axis.x) + s * T(axis.y);
1511  this->zy = tmp * T(axis.y) - s * T(axis.x);
1512  this->zz = tmp * T(axis.z) + c;
1513  this->zw = T(0);
1514  }
1515  this->wx = this->wy = this->wz = T(0);
1516  this->ww = T(1);
1517 }
1518 
1519 template <typename T, Size ROW, Size COL>
1520 inline void Matrix<T,ROW,COL>::set_rotation( const Vector<Float64,3>& axis_v, Float64 angle)
1521 {
1522  mi_static_assert( COL == 4 && ROW == 4);
1523  Vector<T,3> axis( axis_v);
1524  const T min_angle = T(0.00024f);
1525 
1526  if( abs(T(angle)) < min_angle) {
1527  T xa = axis.x * T(angle);
1528  T ya = axis.y * T(angle);
1529  T za = axis.z * T(angle);
1530 
1531  this->xx = T(1);
1532  this->xy = za;
1533  this->xz = -ya;
1534  this->xw = T(0);
1535 
1536  this->yx = za;
1537  this->yy = T(1);
1538  this->yz = xa;
1539  this->yw = T(0);
1540 
1541  this->zx = -ya;
1542  this->zy = -xa;
1543  this->zz = T(1);
1544  this->zw = T(0);
1545  } else {
1546  T s = sin( T(angle));
1547  T c = cos( T(angle));
1548  T t = T(1) - c;
1549  T tmp;
1550 
1551  tmp = t * T(axis.x);
1552  this->xx = tmp * T(axis.x) + c;
1553  this->xy = tmp * T(axis.y) + s * T(axis.z);
1554  this->xz = tmp * T(axis.z) - s * T(axis.y);
1555  this->xw = T(0);
1556 
1557  tmp = t * T(axis.y);
1558  this->yx = tmp * T(axis.x) - s * T(axis.z);
1559  this->yy = tmp * T(axis.y) + c;
1560  this->yz = tmp * T(axis.z) + s * T(axis.x);
1561  this->yw = T(0);
1562 
1563  tmp = t * T(axis.z);
1564  this->zx = tmp * T(axis.x) + s * T(axis.y);
1565  this->zy = tmp * T(axis.y) - s * T(axis.x);
1566  this->zz = tmp * T(axis.z) + c;
1567  this->zw = T(0);
1568  }
1569  this->wx = this->wy = this->wz = T(0);
1570  this->ww = T(1);
1571 }
1572 
1573 template <typename T, Size ROW, Size COL>
1574 inline void Matrix<T,ROW,COL>::lookat(
1575  const Vector<Float32,3>& position,
1576  const Vector<Float32,3>& target,
1577  const Vector<Float32,3>& up)
1578 {
1579  mi_static_assert( COL == 4 && ROW == 4);
1580  Vector<Float32,3> xaxis, yaxis, zaxis;
1581 
1582  // Z vector
1583  zaxis = position - target;
1584  zaxis.normalize();
1585 
1586  // X vector = up cross Z
1587  xaxis = cross( up, zaxis);
1588  xaxis.normalize();
1589 
1590  // Recompute Y = Z cross X
1591  yaxis = cross( zaxis, xaxis);
1592  yaxis.normalize();
1593 
1594  // Build rotation matrix
1595  Matrix<T,4,4> rot(
1596  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1597  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1598  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1599  T(0), T(0), T(0), T(1));
1600 
1601  // Compute the new position by multiplying the inverse position with the rotation matrix
1602  Matrix<T,4,4> trans(
1603  T(1), T(0), T(0), T(0),
1604  T(0), T(1), T(0), T(0),
1605  T(0), T(0), T(1), T(0),
1606  T(-position.x), T(-position.y), T(-position.z), T(1));
1607 
1608  *this = trans * rot;
1609 }
1610 
1611 template <typename T, Size ROW, Size COL>
1612 inline void Matrix<T,ROW,COL>::lookat(
1613  const Vector<Float64,3>& position,
1614  const Vector<Float64,3>& target,
1615  const Vector<Float64,3>& up)
1616 {
1617  mi_static_assert( COL == 4 && ROW == 4);
1618  Vector<Float64,3> xaxis, yaxis, zaxis;
1619 
1620  // Z vector
1621  zaxis = position - target;
1622  zaxis.normalize();
1623 
1624  // X vector = up cross Z
1625  xaxis = cross( up, zaxis);
1626  xaxis.normalize();
1627 
1628  // Recompute Y = Z cross X
1629  yaxis = cross( zaxis, xaxis);
1630  yaxis.normalize();
1631 
1632  // Build rotation matrix
1633  Matrix<T,4,4> rot(
1634  T(xaxis.x), T(yaxis.x), T(zaxis.x), T(0),
1635  T(xaxis.y), T(yaxis.y), T(zaxis.y), T(0),
1636  T(xaxis.z), T(yaxis.z), T(zaxis.z), T(0),
1637  T(0), T(0), T(0), T(1));
1638 
1639  // Compute the new position by multiplying the inverse position with the rotation matrix
1640  Matrix<T,4,4> trans(
1641  T(1), T(0), T(0), T(0),
1642  T(0), T(1), T(0), T(0),
1643  T(0), T(0), T(1), T(0),
1644  T(-position.x), T(-position.y), T(-position.z), T(1));
1645 
1646  *this = trans * rot;
1647 }
1648 
1649 
1650 //------ Definition and helper for matrix inversion ---------------------------
1651 
1652 // Matrix inversion class, specialized for symmetric matrices and low dimensions
1653 template <class T, Size ROW, Size COL>
1654 class Matrix_inverter
1655 {
1656 public:
1657  typedef math::Matrix<T,ROW,COL> Matrix;
1658 
1659  // Inverts the matrix \c M if possible and returns \c true.
1660  //
1661  // If the matrix cannot be inverted or if \c ROW != \c COL, this function returns \c false.
1662  static inline bool invert( Matrix& /* M */) { return false; }
1663 };
1664 
1665 // Matrix inversion class, specialization for symmetric matrices
1666 template <class T, Size DIM>
1667 class Matrix_inverter<T,DIM,DIM>
1668 {
1669 public:
1670  typedef math::Matrix<T,DIM,DIM> Matrix;
1671  typedef math::Vector<T,DIM> Vector;
1672  typedef math::Vector<Size,DIM> Index_vector;
1673 
1674  // LU decomposition of matrix lu in place.
1675  //
1676  // Returns \c false if matrix cannot be decomposed, for example, if it is not invertible, and \c
1677  // true otherwise. Returns also a row permutation in indx.
1678  static bool lu_decomposition(
1679  Matrix& lu, // matrix to decompose, modified in place
1680  Index_vector& indx); // row permutation info indx[4]
1681 
1682  // Solves the equation lu * x = b for x. lu and indx are the results of lu_decomposition. The
1683  // solution is returned in b.
1684  static void lu_backsubstitution(
1685  const Matrix& lu, // LU decomposed matrix
1686  const Index_vector& indx, // permutation vector
1687  Vector& b); // right hand side vector b, modified in place
1688 
1689  static bool invert( Matrix& mat);
1690 };
1691 
1692 template <class T, Size DIM>
1693 bool Matrix_inverter<T,DIM,DIM>::lu_decomposition(
1694  Matrix& lu,
1695  Index_vector& indx)
1696 {
1697  Vector vv; // implicit scaling of each row
1698 
1699  for( Size i = 0; i < DIM; i++) { // get implicit scaling
1700  T big = T(0);
1701  for( Size j = 0; j < DIM; j++) {
1702  T temp = abs(lu.get(i,j));
1703  if( temp > big)
1704  big = temp;
1705  }
1706  if( big == T(0))
1707  return false;
1708  vv[i] = T(1) / big; // save scaling
1709  }
1710  //
1711  // loop over columns of Crout's method
1712  //
1713  Size imax = 0;
1714  for( Size j = 0; j < DIM; j++) {
1715  for( Size i = 0; i < j; i++) {
1716  T sum = lu.get(i,j);
1717  for( Size k = 0; k < i; k++)
1718  sum -= lu.get(i,k) * lu.get(k,j);
1719  lu.set(i, j, sum);
1720  }
1721  T big = 0; // init for pivot search
1722  for( Size i = j; i < DIM; i++) {
1723  T sum = lu.get(i,j);
1724  for( Size k = 0; k < j; k++)
1725  sum -= lu.get(i,k) * lu.get(k,j);
1726  lu.set(i, j, sum);
1727  T dum = vv[i] * abs(sum);
1728  if( dum >= big) { // pivot good?
1729  big = dum;
1730  imax = i;
1731  }
1732  }
1733  if( j != imax) { // interchange rows
1734  for( Size k = 0; k < DIM; k++) {
1735  T dum = lu.get(imax,k);
1736  lu.set(imax, k, lu.get(j,k));
1737  lu.set(j, k, dum);
1738  }
1739  vv[imax] = vv[j]; // interchange scale factor
1740  }
1741  indx[j] = imax;
1742  if( lu.get(j,j) == 0)
1743  return false;
1744  if( j != DIM-1) { // divide by pivot element
1745  T dum = T(1) / lu.get(j,j);
1746  for( Size i = j + 1; i < DIM; i++)
1747  lu.set(i, j, lu.get(i,j) * dum);
1748  }
1749  }
1750  return true;
1751 }
1752 
1753 template <class T, Size DIM>
1754 void Matrix_inverter<T,DIM,DIM>::lu_backsubstitution(
1755  const Matrix& lu,
1756  const Index_vector& indx,
1757  Vector& b)
1758 {
1759  // when ii != DIM+1, ii is index of first non-vanishing element of b
1760  Size ii = DIM+1;
1761  for( Size i = 0; i < DIM; i++) {
1762  Size ip = indx[i];
1763  T sum = b[ip];
1764  b[ip] = b[i];
1765  if( ii != DIM+1) {
1766  for( Size j = ii; j < i; j++) {
1767  sum -= lu.get(i,j) * b[j];
1768  }
1769  } else {
1770  if( sum != T(0)) { // non-zero element
1771  ii = i;
1772  }
1773  }
1774  b[i] = sum;
1775  }
1776  for( Size i2 = DIM; i2 > 0;) { // backsubstitution
1777  --i2;
1778  T sum = b[i2];
1779  for( Size j = i2+1; j < DIM; j++)
1780  sum -= lu.get(i2,j) * b[j];
1781  b[i2] = sum / lu.get(i2,i2); // store comp. of sol. vector
1782  }
1783 }
1784 
1785 template <class T, Size DIM>
1786 bool Matrix_inverter<T,DIM,DIM>::invert( Matrix& mat)
1787 {
1788  Matrix lu(mat); // working copy of matrix to invert
1789  Index_vector indx; // row permutation info
1790 
1791  // LU decompose, return if it fails
1792  if( !lu_decomposition(lu, indx))
1793  return false;
1794 
1795  // solve for each column vector of the I matrix
1796  for( Size j = 0; j < DIM; ++j) {
1797  Vector col(T(0)); // TODO: do that directly in the mat matrix
1798  col[j] = T(1);
1799  lu_backsubstitution( lu, indx, col);
1800  for( Size i = 0; i < DIM; ++i) {
1801  mat.set( i, j, col[i]);
1802  }
1803  }
1804  return true;
1805 }
1806 
1807 // Matrix inversion class, specialization for 1x1 matrix
1808 template <class T>
1809 class Matrix_inverter<T,1,1>
1810 {
1811 public:
1812  typedef math::Matrix<T,1,1> Matrix;
1813 
1814  static inline bool invert( Matrix& mat)
1815  {
1816  T s = mat.get( 0, 0);
1817  if( s == T(0))
1818  return false;
1819  mat.set( 0, 0, T(1) / s);
1820  return true;
1821  }
1822 };
1823 
1824 // Matrix inversion class, specialization for 2x2 matrix
1825 template <class T>
1826 class Matrix_inverter<T,2,2>
1827 {
1828 public:
1829  typedef math::Matrix<T,2,2> Matrix;
1830 
1831  static inline bool invert( Matrix& mat)
1832  {
1833  T a = mat.get( 0, 0);
1834  T b = mat.get( 0, 1);
1835  T c = mat.get( 1, 0);
1836  T d = mat.get( 1, 1);
1837  T det = a*d - b*c;
1838  if( det == T( 0))
1839  return false;
1840  T rdet = T(1) / det;
1841  mat.set( 0, 0, d * rdet);
1842  mat.set( 0, 1,-b * rdet);
1843  mat.set( 1, 0,-c * rdet);
1844  mat.set( 1, 1, a * rdet);
1845  return true;
1846  }
1847 };
1848 
1849 template <typename T, Size ROW, Size COL>
1850 inline bool Matrix<T,ROW,COL>::invert()
1851 {
1852  return Matrix_inverter<T,ROW,COL>::invert( *this);
1853 }
1854 
1855 
1856 
1857 //------ Specializations and (specialized) overloads --------------------------
1858 
1859 // Specialization of common matrix multiplication for 4x4 matrices.
1860 template <typename T>
1861 inline Matrix<T,4,4>& operator*=(
1862  Matrix<T,4,4>& lhs,
1863  const Matrix<T,4,4>& rhs)
1864 {
1865  Matrix<T,4,4> old( lhs);
1866 
1867  lhs.xx = old.xx * rhs.xx + old.xy * rhs.yx + old.xz * rhs.zx + old.xw * rhs.wx;
1868  lhs.xy = old.xx * rhs.xy + old.xy * rhs.yy + old.xz * rhs.zy + old.xw * rhs.wy;
1869  lhs.xz = old.xx * rhs.xz + old.xy * rhs.yz + old.xz * rhs.zz + old.xw * rhs.wz;
1870  lhs.xw = old.xx * rhs.xw + old.xy * rhs.yw + old.xz * rhs.zw + old.xw * rhs.ww;
1871 
1872  lhs.yx = old.yx * rhs.xx + old.yy * rhs.yx + old.yz * rhs.zx + old.yw * rhs.wx;
1873  lhs.yy = old.yx * rhs.xy + old.yy * rhs.yy + old.yz * rhs.zy + old.yw * rhs.wy;
1874  lhs.yz = old.yx * rhs.xz + old.yy * rhs.yz + old.yz * rhs.zz + old.yw * rhs.wz;
1875  lhs.yw = old.yx * rhs.xw + old.yy * rhs.yw + old.yz * rhs.zw + old.yw * rhs.ww;
1876 
1877  lhs.zx = old.zx * rhs.xx + old.zy * rhs.yx + old.zz * rhs.zx + old.zw * rhs.wx;
1878  lhs.zy = old.zx * rhs.xy + old.zy * rhs.yy + old.zz * rhs.zy + old.zw * rhs.wy;
1879  lhs.zz = old.zx * rhs.xz + old.zy * rhs.yz + old.zz * rhs.zz + old.zw * rhs.wz;
1880  lhs.zw = old.zx * rhs.xw + old.zy * rhs.yw + old.zz * rhs.zw + old.zw * rhs.ww;
1881 
1882  lhs.wx = old.wx * rhs.xx + old.wy * rhs.yx + old.wz * rhs.zx + old.ww * rhs.wx;
1883  lhs.wy = old.wx * rhs.xy + old.wy * rhs.yy + old.wz * rhs.zy + old.ww * rhs.wy;
1884  lhs.wz = old.wx * rhs.xz + old.wy * rhs.yz + old.wz * rhs.zz + old.ww * rhs.wz;
1885  lhs.ww = old.wx * rhs.xw + old.wy * rhs.yw + old.wz * rhs.zw + old.ww * rhs.ww;
1886 
1887  return lhs;
1888 }
1889 
1890 // Specialization of common matrix multiplication for 4x4 matrices.
1891 template <typename T>
1892 inline Matrix<T,4,4> operator*(
1893  const Matrix<T,4,4>& lhs,
1894  const Matrix<T,4,4>& rhs)
1895 {
1896  Matrix<T,4,4> temp( lhs);
1897  temp *= rhs;
1898  return temp;
1899 }
1900 
1901 #endif // MI_FOR_DOXYGEN_ONLY
1902 
1903 #ifdef MI_NEURAYLIB_DEPRECATED_5_0
1904 
1905 template <typename T, Size ROW, Size COL>
1906 inline void Matrix<T,ROW,COL>::multiply( const Matrix<Float32,4,4>& matrix)
1907 {
1908  mi_static_assert( COL == 4 && ROW == 4);
1909  (*this) *= Matrix<T,ROW,COL> (matrix);
1910 }
1911 
1912 template <typename T, Size ROW, Size COL>
1913 inline void Matrix<T,ROW,COL>::multiply( const Matrix<Float64,4,4>& matrix)
1914 {
1915  mi_static_assert( COL == 4 && ROW == 4);
1916  (*this) *= Matrix<T,ROW,COL>( matrix);
1917 }
1918 
1919 template <typename T, Size ROW, Size COL>
1920 inline Vector<Float32,3> Matrix<T,ROW,COL>::transform( const Vector<Float32,3> vector) const
1921 {
1922  mi_static_assert( COL == 4 && ROW == 4);
1923  return transform_point( *this, vector);
1924 }
1925 
1926 template <typename T, Size ROW, Size COL>
1927 inline Vector<Float64,3> Matrix<T,ROW,COL>::transform( const Vector<Float64,3> vector) const
1928 {
1929  mi_static_assert( COL == 4 && ROW == 4);
1930  return transform_point( *this, vector);
1931 }
1932 
1933 template <typename T, Size ROW, Size COL>
1934 inline Vector<Float32,3> Matrix<T,ROW,COL>::transform33( const Vector<Float32,3> vector) const
1935 {
1936  mi_static_assert( COL == 4 && ROW == 4);
1937  return transform_vector( *this, vector);
1938 }
1939 
1940 template <typename T, Size ROW, Size COL>
1941 inline Vector<Float64,3> Matrix<T,ROW,COL>::transform33( const Vector<Float64,3> vector) const
1942 {
1943  mi_static_assert( COL == 4 && ROW == 4);
1944  return transform_vector( *this, vector);
1945 }
1946 
1947 template <typename T, Size ROW, Size COL>
1948 inline Vector<Float32,3> Matrix<T,ROW,COL>::transform33t( const Vector<Float32,3> vector) const
1949 {
1950  mi_static_assert( COL == 4 && ROW == 4);
1951  return transform_normal_inv( *this, vector);
1952 }
1953 
1954 template <typename T, Size ROW, Size COL>
1955 inline Vector<Float64,3> Matrix<T,ROW,COL>::transform33t( const Vector<Float64,3> vector) const
1956 {
1957  mi_static_assert( COL == 4 && ROW == 4);
1958  return transform_normal_inv( *this, vector);
1959 }
1960 
1961 #endif // MI_NEURAYLIB_DEPRECATED_5_0
1962  // end group mi_math_matrix
1964 
1965 } // namespace math
1966 
1967 } // namespace mi
1968 
1969 #endif // MI_MATH_MATRIX_H