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