Mitsuba Renderer  0.5.0
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
matrix.h
Go to the documentation of this file.
1 /*
2  This file is part of Mitsuba, a physically based rendering system.
3 
4  Copyright (c) 2007-2014 by Wenzel Jakob and others.
5 
6  Mitsuba is free software; you can redistribute it and/or modify
7  it under the terms of the GNU General Public License Version 3
8  as published by the Free Software Foundation.
9 
10  Mitsuba is distributed in the hope that it will be useful,
11  but WITHOUT ANY WARRANTY; without even the implied warranty of
12  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13  GNU General Public License for more details.
14 
15  You should have received a copy of the GNU General Public License
16  along with this program. If not, see <http://www.gnu.org/licenses/>.
17 */
18 
19 #pragma once
20 #if !defined(__MITSUBA_CORE_MATRIX_H_)
21 #define __MITSUBA_CORE_MATRIX_H_
22 
23 #include <mitsuba/mitsuba.h>
24 #include <boost/static_assert.hpp>
25 
27 
28 /**
29  * \brief Generic fixed-size dense matrix class using a row-major storage format
30  *
31  * \ingroup libcore
32  */
33 template <int M, int N, typename T> struct Matrix {
34 public:
35  T m[M][N];
36 
37  /**
38  * \brief Construct a new MxN matrix without initializing it.
39  *
40  * This constructor is useful when the matrix will either not
41  * be used at all (it might be part of a larger data structure)
42  * or initialized at a later point in time. Always make sure
43  * that one of the two is the case! Otherwise your program will do
44  * computations involving uninitialized memory, which will probably
45  * lead to a difficult-to-find bug.
46  */
47 #if !defined(MTS_DEBUG_UNINITIALIZED)
48  Matrix() { }
49 #else
50  Matrix() {
51  for (int i=0; i<M; ++i)
52  for (int j=0; j<N; ++j)
53  m[i][j] = std::numeric_limits<double>::quiet_NaN();
54  }
55 #endif
56 
57  /// Initialize the matrix with constant entries
58  explicit inline Matrix(T value) {
59  for (int i=0; i<M; ++i)
60  for (int j=0; j<N; ++j)
61  m[i][j] = value;
62  }
63 
64  /// Initialize the matrix from a given MxN array
65  explicit inline Matrix(const T _m[M][N]) {
66  memcpy(m, _m, sizeof(T) * M * N);
67  }
68 
69  /// Initialize the matrix from a given (flat) MxN array in row-major order
70  explicit inline Matrix(const T _m[M*N]) {
71  memcpy(m, _m, sizeof(T) * M * N);
72  }
73 
74  /// Unserialize a matrix from a stream
75  explicit inline Matrix(Stream *stream) {
76  stream->readArray(&m[0][0], M * N);
77  }
78 
79  /// Copy constructor
80  inline Matrix(const Matrix &mtx) {
81  memcpy(m, mtx.m, sizeof(T) * M * N);
82  }
83 
84  /// Initialize with the identity matrix
85  void setIdentity() {
86  for (int i=0; i<M; ++i)
87  for (int j=0; j<N; ++j)
88  m[i][j] = (i == j) ? 1.0f : 0.0f;
89  }
90 
91  /// Initialize with zeroes
92  void setZero() {
93  memset(m, 0, sizeof(T) * M * N);
94  }
95 
96 
97  /// Indexing operator
98  inline T &operator()(int i, int j) { return m[i][j]; }
99 
100  /// Indexing operator (const verions)
101  inline const T & operator()(int i, int j) const { return m[i][j]; }
102 
103  /// Equality operator
104  inline bool operator==(const Matrix &mat) const {
105  for (int i=0; i<M; ++i)
106  for (int j=0; j<N; ++j)
107  if (m[i][j] != mat.m[i][j])
108  return false;
109  return true;
110  }
111 
112  /// Inequality operator
113  inline bool operator!=(const Matrix &mat) const {
114  for (int i=0; i<M; ++i)
115  for (int j=0; j<N; ++j)
116  if (m[i][j] != mat.m[i][j])
117  return true;
118  return false;
119  }
120 
121  /// Assignment operator
122  inline Matrix &operator=(const Matrix &mat) {
123  for (int i=0; i<M; ++i)
124  for (int j=0; j<N; ++j)
125  m[i][j] = mat.m[i][j];
126  return *this;
127  }
128 
129  /// Matrix addition (returns a temporary)
130  inline Matrix operator+(const Matrix &mat) const {
131  Matrix result;
132  for (int i=0; i<M; ++i)
133  for (int j=0; j<N; ++j)
134  result.m[i][j] = m[i][j] + mat.m[i][j];
135  return result;
136  }
137 
138  /// Matrix-scalar addition (returns a temporary)
139  inline Matrix operator+(T value) const {
140  Matrix result;
141  for (int i=0; i<M; ++i)
142  for (int j=0; j<N; ++j)
143  result.m[i][j] = m[i][j] + value;
144  return result;
145  }
146 
147  /// Matrix addition
148  inline const Matrix &operator+=(const Matrix &mat) {
149  for (int i=0; i<M; ++i)
150  for (int j=0; j<N; ++j)
151  m[i][j] += mat.m[i][j];
152  return *this;
153  }
154 
155  /// Matrix-scalar addition
156  inline const Matrix &operator+=(T value) {
157  for (int i=0; i<M; ++i)
158  for (int j=0; j<N; ++j)
159  m[i][j] += value;
160  return *this;
161  }
162 
163  /// Matrix subtraction (returns a temporary)
164  inline Matrix operator-(const Matrix &mat) const {
165  Matrix result;
166  for (int i=0; i<M; ++i)
167  for (int j=0; j<N; ++j)
168  result.m[i][j] = m[i][j] - mat.m[i][j];
169  return result;
170  }
171 
172  /// Matrix-scalar subtraction (returns a temporary)
173  inline Matrix operator-(T value) const {
174  Matrix result;
175  for (int i=0; i<M; ++i)
176  for (int j=0; j<N; ++j)
177  result.m[i][j] = m[i][j] - value;
178  return result;
179  }
180 
181  /// Matrix subtraction
182  inline const Matrix &operator-=(const Matrix &mat) {
183  for (int i=0; i<M; ++i)
184  for (int j=0; j<N; ++j)
185  m[i][j] -= mat.m[i][j];
186  return *this;
187  }
188 
189  /// Matrix-scalar subtraction
190  inline const Matrix &operator-(T value) {
191  for (int i=0; i<M; ++i)
192  for (int j=0; j<N; ++j)
193  m[i][j] -= value;
194  return *this;
195  }
196 
197  /// Matrix-scalar addition
198  inline const Matrix &operator-=(T value) {
199  for (int i=0; i<M; ++i)
200  for (int j=0; j<N; ++j)
201  m[i][j] -= value;
202  return *this;
203  }
204 
205  /// Component-wise negation
206  inline Matrix operator-() const {
207  Matrix result;
208  for (int i=0; i<M; ++i)
209  for (int j=0; j<N; ++j)
210  result.m[i][j] = -m[i][j];
211  return result;
212  }
213 
214  /// Scalar multiplication (creates a temporary)
215  inline Matrix operator*(T value) const {
216  Matrix result;
217  for (int i=0; i<M; ++i)
218  for (int j=0; j<N; ++j)
219  result.m[i][j] = m[i][j]*value;
220  return result;
221  }
222 
223  /// Scalar multiplication
224  inline const Matrix& operator*=(T value) {
225  for (int i=0; i<M; ++i)
226  for (int j=0; j<N; ++j)
227  m[i][j] *= value;
228  return *this;
229  }
230 
231  /// Scalar division (creates a temporary)
232  inline Matrix operator/(T value) const {
233  Matrix result;
234 #ifdef MTS_DEBUG
235  if (value == 0)
236  SLog(EWarn, "Matrix: Division by zero!");
237 #endif
238  Float recip = 1/value;
239  for (int i=0; i<M; ++i)
240  for (int j=0; j<N; ++j)
241  result.m[i][j] = m[i][j]*recip;
242  return result;
243  }
244 
245  /// Scalar division
246  inline const Matrix& operator/=(T value) {
247 #ifdef MTS_DEBUG
248  if (value == 0)
249  SLog(EWarn, "Matrix: Division by zero!");
250 #endif
251  Float recip = 1/value;
252  for (int i=0; i<M; ++i)
253  for (int j=0; j<N; ++j)
254  m[i][j] *= recip;
255  return *this;
256  }
257 
258  /// Matrix multiplication (for square matrices)
259  inline const Matrix &operator*=(const Matrix &mat) {
260  BOOST_STATIC_ASSERT(M == N);
261  Matrix temp = *this * mat;
262  *this = temp;
263  return *this;
264  }
265 
266  /// Compute the trace of a square matrix
267  inline Float trace() const {
268  BOOST_STATIC_ASSERT(M == N);
269  Float sum = 0;
270  for (int i=0; i<M; ++i)
271  sum += m[i][i];
272  return sum;
273  }
274 
275  /// Compute the Frobenius norm
276  inline Float frob() const {
277  Float val = 0;
278  for (int i=0; i<M; ++i)
279  for (int j=0; j<N; ++j)
280  val += m[i][j] * m[i][j];
281  return std::sqrt(val);
282  }
283 
284  /**
285  * \brief Compute the LU decomposition of a matrix
286  *
287  * For an m-by-n matrix A with m >= n, the LU decomposition is an
288  * m-by-n unit lower triangular matrix L, an n-by-n upper triangular
289  * matrix U,
290  *
291  * and a permutation vector piv of length m so that A(piv,:) = L*U.
292  * If m < n, then L is m-by-m and U is m-by-n.
293  *
294  * The LU decomposition with pivoting always exists, even if the matrix is
295  * singular, so the constructor will never fail.
296  * The primary use of the
297  *
298  * LU decomposition is in the solution of square systems of simultaneous
299  * linear equations.
300  *
301  * \param Target matrix (the L and U parts will be stored
302  * together in a packed format)
303  * \param piv Storage for the permutation created by the pivoting
304  * \param pivsign Sign of the permutation
305  * \return \c true if the matrix was nonsingular.
306  *
307  * Based on the implementation in JAMA.
308  */
309  bool lu(Matrix &LU, int piv[M], int &pivsign) const;
310 
311  /**
312  * Compute the Cholesky decomposition of a symmetric
313  * positive definite matrix
314  *
315  * \param L Target matrix (a lower triangular matrix such
316  * that A=L*L')
317  * \return \c false If the matrix is not symmetric positive
318  * definite.
319  * Based on the implementation in JAMA.
320  */
321  bool chol(Matrix &L) const;
322 
323  /**
324  * Solve A*X==B, where \a this is a Cholesky decomposition of \a A created by \ref chol()
325  *
326  * \param B A matrix with as many rows as \a A and any number of columns
327  * \param X A matrix such that L*L'*X == B
328  *
329  * Based on the implementation in JAMA.
330  */
331  template <int K> void cholSolve(const Matrix<M, K, T> &B,
332  Matrix<M, K, T> &X) const;
333 
334  /**
335  * Solve A*X==B, where \a this is a LU decomposition of \a A created by \ref lu()
336  *
337  * \param B A matrix with as many rows as \a A and any number of columns
338  * \param X A matrix such that L*U*X == B(piv, :)
339  * \param piv Pivot vector returned by \ref lu()
340  *
341  * Based on the implementation in JAMA.
342  */
343  template <int K> void luSolve(const Matrix<M, K, T> &B,
344  Matrix<M, K, T> &X, int piv[M]) const;
345 
346  /**
347  * \brief Compute the determinant of a decomposed matrix
348  * created by \ref lu()
349  *
350  * \param pivsign The sign of the pivoting permutation returned
351  * by \ref lu()
352  *
353  * Based on the implementation in JAMA.
354  */
355  T luDet(int pivsign) const;
356 
357  /**
358  * \brief Compute the determinant of a decomposed matrix
359  * created by \ref chol()
360  */
361 
362  T cholDet() const;
363 
364  /// Check if the matrix is identically zero
365  inline bool isZero() const {
366  for (int i=0; i<M; ++i)
367  for (int j=0; j<N; ++j)
368  if (m[i][j] != 0)
369  return false;
370  return true;
371  }
372 
373  /// Test if this is the identity matrix
374  inline bool isIdentity() const {
375  for (int i=0; i<M; ++i) {
376  for (int j=0; j<N; ++j) {
377  if (m[i][j] != ((i==j) ? 1 : 0))
378  return false;
379  }
380  }
381  return true;
382  }
383 
384  /**
385  * \brief Compute the determinant of a square matrix (internally
386  * creates a LU decomposition)
387  */
388  inline T det() const {
389  Matrix LU;
390  int piv[M], pivsign;
391  if (!lu(LU, piv, pivsign))
392  return 0.0f;
393  return LU.luDet(pivsign);
394  }
395 
396  /// Compute the inverse of a square matrix using the Gauss-Jordan algorithm
397  bool invert(Matrix &target) const;
398 
399  /**
400  * \brief Perform a symmetric eigendecomposition of a square matrix
401  * into Q and D.
402  *
403  * Based on the implementation in JAMA.
404  */
405  inline void symEig(Matrix &Q, T d[M]) const {
406  BOOST_STATIC_ASSERT(M == N);
407  T e[M];
408  Q = *this;
409  tred2(Q.m, d, e);
410  tql2(Q.m, d, e);
411  }
412 
413  /// Compute the transpose of this matrix
414  inline void transpose(Matrix<N, M, T> &target) const {
415  for (int i=0; i<M; ++i)
416  for (int j=0; j<N; ++j)
417  target.m[i][j] = m[j][i];
418  }
419 
420  /// Serialize the matrix to a stream
421  inline void serialize(Stream *stream) const {
422  stream->writeArray(&m[0][0], M*N);
423  }
424 
425  /// Return a string representation
426  std::string toString() const {
427  std::ostringstream oss;
428  oss << "Matrix" << M << "x" << N << "["<< std::endl;
429  for (int i=0; i<M; ++i) {
430  oss << " ";
431  for (int j=0; j<N; ++j) {
432  oss << m[i][j];
433  if (j != N-1)
434  oss << ", ";
435  }
436  if (i != M-1)
437  oss << ";";
438  oss << std::endl;
439  }
440  oss << "]";
441  return oss.str();
442  }
443 protected:
444  /// Symmetric Householder reduction to tridiagonal form.
445  static void tred2(T V[M][N], T d[N], T e[N]);
446 
447  /// Symmetric tridiagonal QL algorithm.
448  static void tql2(T V[M][N], T d[N], T e[N]);
449 };
450 
451 /**
452  * \brief Basic 2x2 matrix data type
453  * \ingroup libcore
454  */
455 struct MTS_EXPORT_CORE Matrix2x2 : public Matrix<2, 2, Float> {
456 public:
457  inline Matrix2x2() { }
458 
459  /// Initialize the matrix with constant entries
460  explicit inline Matrix2x2(Float value) : Matrix<2, 2, Float>(value) { }
461 
462  /// Initialize the matrix from a given 2x2 array
463  explicit inline Matrix2x2(const Float _m[2][2]) : Matrix<2, 2, Float>(_m) { }
464 
465  /// Initialize the matrix from a given (float) 2x2 array in row-major order
466  explicit inline Matrix2x2(const Float _m[4]) : Matrix<2, 2, Float>(_m) { }
467 
468  /// Initialize the matrix from two 2D column vectors
469  explicit inline Matrix2x2(const Vector2 &v1, const Vector2 &v2) {
470  m[0][0] = v1.x; m[0][1] = v2.x;
471  m[1][0] = v1.y; m[1][1] = v2.y;
472  }
473 
474  /// Unserialize a matrix from a stream
475  explicit inline Matrix2x2(Stream *stream) : Matrix<2, 2, Float>(stream) { }
476 
477  /// Copy constructor
478  inline Matrix2x2(const Matrix<2, 2, Float> &mtx) : Matrix<2, 2, Float>(mtx) { }
479 
480  /// Initialize with the given values
481  inline Matrix2x2(Float a00, Float a01, Float a10, Float a11) {
482  m[0][0] = a00; m[0][1] = a01;
483  m[1][0] = a10; m[1][1] = a11;
484  }
485 
486  /// Return the determinant (Faster than Matrix::det)
487  inline Float det() const {
488  return m[0][0]*m[1][1] - m[0][1]*m[1][0];
489  }
490 
491  /// Compute the inverse (Faster than Matrix::invert)
492  FINLINE bool invert2x2(Matrix2x2 &target) const {
493  Float det = m[0][0]*m[1][1] - m[0][1]*m[1][0];
494  if (std::abs(det) <= RCPOVERFLOW)
495  return false;
496  Float invDet = 1/det;
497  target.m[0][0] = m[1][1] * invDet;
498  target.m[0][1] = -m[0][1] * invDet;
499  target.m[1][1] = m[0][0] * invDet;
500  target.m[1][0] = -m[1][0] * invDet;
501  return true;
502  }
503  FINLINE bool invert(Matrix2x2 &target) const {
504  return invert2x2(target);
505  }
506 
507  /// Compute the inverse with det (Faster than Matrix::invert)
508  FINLINE bool invert2x2(Matrix2x2 &target, Float& det) const {
509  det = m[0][0]*m[1][1] - m[0][1]*m[1][0];
510  if(std::abs(det) <= RCPOVERFLOW)
511  return false;
512  Float invDet = 1/det;
513  target.m[0][0] = m[1][1] * invDet;
514  target.m[0][1] = -m[0][1] * invDet;
515  target.m[1][1] = m[0][0] * invDet;
516  target.m[1][0] = -m[1][0] * invDet;
517  return true;
518  }
519 
520  /// Matrix-vector multiplication
521  inline Vector2 operator*(const Vector2 &v) const {
522  return Vector2(
523  m[0][0] * v.x + m[0][1] * v.y,
524  m[1][0] * v.x + m[1][1] * v.y
525  );
526  }
527 
528  /// Scalar multiplication (creates a temporary)
529  inline Matrix2x2 operator*(Float value) const {
530  Matrix2x2 result;
531  for (int i=0; i<2; ++i)
532  for (int j=0; j<2; ++j)
533  result.m[i][j] = m[i][j]*value;
534  return result;
535  }
536 
537  /// Assignment operator
539  for (int i=0; i<2; ++i)
540  for (int j=0; j<2; ++j)
541  m[i][j] = mat.m[i][j];
542  return *this;
543  }
544 
545  /// Return a row by index
546  inline Vector2 row(int i) const {
547  return Vector2(m[i][0], m[i][1]);
548  }
549 
550  /// Return a column by index
551  inline Vector2 col(int i) const {
552  return Vector2(m[0][i], m[1][i]);
553  }
554 };
555 
556 /**
557  * \brief Basic 3x3 matrix data type
558  * \ingroup libcore
559  */
560 struct MTS_EXPORT_CORE Matrix3x3 : public Matrix<3, 3, Float> {
561 public:
562  inline Matrix3x3() { }
563 
564  /// Initialize the matrix with constant entries
565  explicit inline Matrix3x3(Float value) : Matrix<3, 3, Float>(value) { }
566 
567  /// Initialize the matrix from a given 3x3 array
568  explicit inline Matrix3x3(const Float _m[3][3]) : Matrix<3, 3, Float>(_m) { }
569 
570  /// Initialize the matrix from a given (float) 3x3 array in row-major order
571  explicit inline Matrix3x3(const Float _m[9]) : Matrix<3, 3, Float>(_m) { }
572 
573  /// Initialize the matrix from three 3D column vectors
574  explicit inline Matrix3x3(const Vector &v1, const Vector &v2, const Vector &v3) {
575  m[0][0] = v1.x; m[0][1] = v2.x; m[0][2] = v3.x;
576  m[1][0] = v1.y; m[1][1] = v2.y; m[1][2] = v3.y;
577  m[2][0] = v1.z; m[2][1] = v2.z; m[2][2] = v3.z;
578  }
579 
580  /// Unserialize a matrix from a stream
581  explicit inline Matrix3x3(Stream *stream) : Matrix<3, 3, Float>(stream) { }
582 
583  /// Copy constructor
584  inline Matrix3x3(const Matrix<3, 3, Float> &mtx) : Matrix<3, 3, Float>(mtx) { }
585 
586  /// Initialize with the given values
587  inline Matrix3x3(Float a00, Float a01, Float a02,
588  Float a10, Float a11, Float a12,
589  Float a20, Float a21, Float a22) {
590  m[0][0] = a00; m[0][1] = a01; m[0][2] = a02;
591  m[1][0] = a10; m[1][1] = a11; m[1][2] = a12;
592  m[2][0] = a20; m[2][1] = a21; m[2][2] = a22;
593  }
594 
595  /// Return the determinant (Faster than Matrix::det())
596  inline Float det() const {
597  return ((m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]))
598  - (m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]))
599  + (m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])));
600  }
601 
602  /// Matrix-vector multiplication
603  inline Vector operator*(const Vector &v) const {
604  return Vector(
605  m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z,
606  m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z,
607  m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z);
608  }
609 
610  /// multiply the matrix by transpose of v, returns the transpose of the line vector.
611  inline Vector preMult(const Vector &v) const {
612  return Vector(v.x * m[0][0] + v.y * m[1][0] + v.z * m[2][0],
613  v.x * m[0][1] + v.y * m[1][1] + v.z * m[2][1],
614  v.x * m[0][2] + v.y * m[1][2] + v.z * m[2][2]);
615  }
616 
617  /// Scalar multiplication (creates a temporary)
618  inline Matrix3x3 operator*(Float value) const {
619  Matrix3x3 result;
620  for (int i=0; i<3; ++i)
621  for (int j=0; j<3; ++j)
622  result.m[i][j] = m[i][j]*value;
623  return result;
624  }
625 
626 
627  /// Assignment operator
629  for (int i=0; i<3; ++i)
630  for (int j=0; j<3; ++j)
631  m[i][j] = mat.m[i][j];
632  return *this;
633  }
634 
635  /// Return a row by index
636  inline Vector3 row(int i) const {
637  return Vector3(
638  m[i][0], m[i][1], m[i][2]
639  );
640  }
641 
642  /// Return a column by index
643  inline Vector3 col(int i) const {
644  return Vector3(
645  m[0][i], m[1][i], m[2][i]
646  );
647  }
648 };
649 
650 
651 /**
652  * \brief Basic 4x4 matrix data type
653  * \ingroup libcore
654  * \ingroup libpython
655  */
656 struct MTS_EXPORT_CORE Matrix4x4 : public Matrix<4, 4, Float> {
657  inline Matrix4x4() { }
658 
659  /// Initialize the matrix with constant entries
660  explicit inline Matrix4x4(Float value) : Matrix<4, 4, Float>(value) { }
661 
662  /// Initialize the matrix from a given 4x4 array
663  explicit inline Matrix4x4(const Float _m[4][4]) : Matrix<4, 4, Float>(_m) { }
664 
665  /// Initialize the matrix from a given (float) 4x4 array in row-major order
666  explicit inline Matrix4x4(const Float _m[16]) : Matrix<4, 4, Float>(_m) { }
667 
668  /// Initialize the matrix from four 4D column vectors
669  explicit inline Matrix4x4(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3, const Vector4 &v4) {
670  m[0][0] = v1.x; m[0][1] = v2.x; m[0][2] = v3.x; m[0][3] = v4.x;
671  m[1][0] = v1.y; m[1][1] = v2.y; m[1][2] = v3.y; m[1][3] = v4.y;
672  m[2][0] = v1.z; m[2][1] = v2.z; m[2][2] = v3.z; m[2][3] = v4.z;
673  m[3][0] = v1.w; m[3][1] = v2.w; m[3][2] = v3.w; m[3][3] = v4.w;
674  }
675 
676  /// Unserialize a matrix from a stream
677  explicit inline Matrix4x4(Stream *stream) : Matrix<4, 4, Float>(stream) { }
678 
679  /// Copy constructor
680  inline Matrix4x4(const Matrix<4, 4, Float> &mtx) : Matrix<4, 4, Float>(mtx) { }
681 
682  /// Initialize with the given values
683  inline Matrix4x4(
684  Float a00, Float a01, Float a02, Float a03,
685  Float a10, Float a11, Float a12, Float a13,
686  Float a20, Float a21, Float a22, Float a23,
687  Float a30, Float a31, Float a32, Float a33) {
688  m[0][0] = a00; m[0][1] = a01; m[0][2] = a02; m[0][3] = a03;
689  m[1][0] = a10; m[1][1] = a11; m[1][2] = a12; m[1][3] = a13;
690  m[2][0] = a20; m[2][1] = a21; m[2][2] = a22; m[2][3] = a23;
691  m[3][0] = a30; m[3][1] = a31; m[3][2] = a32; m[3][3] = a33;
692  }
693 
694  /// Return the determinant of the upper left 3x3 sub-matrix
695  inline Float det3x3() const {
696  return ((m[0][0] * (m[1][1] * m[2][2] - m[1][2] * m[2][1]))
697  - (m[0][1] * (m[1][0] * m[2][2] - m[1][2] * m[2][0]))
698  + (m[0][2] * (m[1][0] * m[2][1] - m[1][1] * m[2][0])));
699  }
700 
701  /// Matrix-vector multiplication
702  inline Vector4 operator*(const Vector4 &v) const {
703  return Vector4(
704  m[0][0] * v.x + m[0][1] * v.y + m[0][2] * v.z + m[0][3] * v.w,
705  m[1][0] * v.x + m[1][1] * v.y + m[1][2] * v.z + m[1][3] * v.w,
706  m[2][0] * v.x + m[2][1] * v.y + m[2][2] * v.z + m[2][3] * v.w,
707  m[3][0] * v.x + m[3][1] * v.y + m[3][2] * v.z + m[3][3] * v.w
708  );
709  }
710 
711  /// Scalar multiplication (creates a temporary)
712  inline Matrix4x4 operator*(Float value) const {
713  Matrix4x4 result;
714  for (int i=0; i<4; ++i)
715  for (int j=0; j<4; ++j)
716  result.m[i][j] = m[i][j]*value;
717  return result;
718  }
719 
720  /// Assignment operator
722  for (int i=0; i<4; ++i)
723  for (int j=0; j<4; ++j)
724  m[i][j] = mat.m[i][j];
725  return *this;
726  }
727 
728  /// Return a row by index
729  inline Vector4 row(int i) const {
730  return Vector4(
731  m[i][0], m[i][1], m[i][2], m[i][3]
732  );
733  }
734 
735  /// Return a column by index
736  inline Vector4 col(int i) const {
737  return Vector4(
738  m[0][i], m[1][i], m[2][i], m[3][i]
739  );
740  }
741 };
742 
743 /// Matrix multiplication (creates a temporary)
744 template <typename T, int M1, int N1, int M2, int N2> inline Matrix<M1, N2, T>
745  operator*(const Matrix<M1, N1, T> &mat1, const Matrix<M2, N2, T> &mat2) {
746  BOOST_STATIC_ASSERT(N1 == M2);
747  Matrix<M1, N2, T> result;
748  for (int i=0; i<M1; ++i) {
749  for (int j=0; j<N2; ++j) {
750  T sum = 0;
751  for (int k=0; k<N1; ++k)
752  sum += mat1.m[i][k] * mat2.m[k][j];
753  result.m[i][j] = sum;
754  }
755  }
756  return result;
757 }
758 
759 template <typename T, int M, int N> inline Matrix<M, N, T> operator*(T f, const Matrix<M, N, T> &m) {
760  return m*f;
761 }
762 
763 /**
764  * \brief Fast 3x3 eigenvalue decomposition
765  *
766  * \param m
767  * Matrix in question -- will be replaced with the eigenvectors
768  * \param lambda
769  * Parameter used to returns the eigenvalues
770  * \return
771  * \c true upon success.
772  */
773 extern MTS_EXPORT_CORE bool eig3(Matrix3x3 &m, Float lambda[3]);
774 
775 /**
776  * \brief Fast non-iterative 3x3 eigenvalue decomposition
777  *
778  * \param m
779  * Matrix in question -- will be replaced with the eigenvectors
780  * \param lambda
781  * Parameter used to returns the eigenvalues
782  */
783 extern MTS_EXPORT_CORE void eig3_noniter(Matrix3x3 &m, Float lambda[3]);
784 
786 
787 #include "matrix.inl"
788 
789 #endif /* __MITSUBA_CORE_MATRIX_H_ */
void serialize(Stream *stream) const
Serialize the matrix to a stream.
Definition: matrix.h:421
Vector operator*(const Vector &v) const
Matrix-vector multiplication.
Definition: matrix.h:603
std::string toString() const
Return a string representation.
Definition: matrix.h:426
Matrix< M, N, T > operator*(T f, const Matrix< M, N, T > &m)
Definition: matrix.h:759
Float trace() const
Compute the trace of a square matrix.
Definition: matrix.h:267
Matrix operator+(const Matrix &mat) const
Matrix addition (returns a temporary)
Definition: matrix.h:130
bool operator!=(const Matrix &mat) const
Inequality operator.
Definition: matrix.h:113
Matrix3x3(const Float _m[3][3])
Initialize the matrix from a given 3x3 array.
Definition: matrix.h:568
Vector3 row(int i) const
Return a row by index.
Definition: matrix.h:636
Matrix3x3 operator*(Float value) const
Scalar multiplication (creates a temporary)
Definition: matrix.h:618
Matrix2x2(const Matrix< 2, 2, Float > &mtx)
Copy constructor.
Definition: matrix.h:478
Matrix operator*(T value) const
Scalar multiplication (creates a temporary)
Definition: matrix.h:215
FINLINE bool invert(Matrix2x2 &target) const
Definition: matrix.h:503
FINLINE bool invert2x2(Matrix2x2 &target) const
Compute the inverse (Faster than Matrix::invert)
Definition: matrix.h:492
Vector4 row(int i) const
Return a row by index.
Definition: matrix.h:729
Definition: fwd.h:48
Matrix2x2(Float value)
Initialize the matrix with constant entries.
Definition: matrix.h:460
Matrix4x4(Float value)
Initialize the matrix with constant entries.
Definition: matrix.h:660
Generic fixed-size dense matrix class using a row-major storage format.
Definition: matrix.h:33
Matrix4x4 operator*(Float value) const
Scalar multiplication (creates a temporary)
Definition: matrix.h:712
Matrix4x4(Stream *stream)
Unserialize a matrix from a stream.
Definition: matrix.h:677
Matrix4x4(const Vector4 &v1, const Vector4 &v2, const Vector4 &v3, const Vector4 &v4)
Initialize the matrix from four 4D column vectors.
Definition: matrix.h:669
Matrix4x4(Float a00, Float a01, Float a02, Float a03, Float a10, Float a11, Float a12, Float a13, Float a20, Float a21, Float a22, Float a23, Float a30, Float a31, Float a32, Float a33)
Initialize with the given values.
Definition: matrix.h:683
Matrix(T value)
Initialize the matrix with constant entries.
Definition: matrix.h:58
Vector3 col(int i) const
Return a column by index.
Definition: matrix.h:643
T & operator()(int i, int j)
Indexing operator.
Definition: matrix.h:98
Matrix3x3 & operator=(const Matrix< 3, 3, Float > &mat)
Assignment operator.
Definition: matrix.h:628
TVector3< Float > Vector
Definition: fwd.h:113
void eig3_noniter(Matrix3x3 &m, Float lambda[3])
Fast non-iterative 3x3 eigenvalue decomposition.
Matrix3x3(Stream *stream)
Unserialize a matrix from a stream.
Definition: matrix.h:581
Vector preMult(const Vector &v) const
multiply the matrix by transpose of v, returns the transpose of the line vector.
Definition: matrix.h:611
Matrix3x3(Float a00, Float a01, Float a02, Float a10, Float a11, Float a12, Float a20, Float a21, Float a22)
Initialize with the given values.
Definition: matrix.h:587
Basic 2x2 matrix data type.
Definition: matrix.h:455
Matrix operator-(const Matrix &mat) const
Matrix subtraction (returns a temporary)
Definition: matrix.h:164
#define MTS_EXPORT_CORE
Definition: getopt.h:29
Matrix2x2(const Float _m[4])
Initialize the matrix from a given (float) 2x2 array in row-major order.
Definition: matrix.h:466
#define SLog(level, fmt,...)
Write a Log message to the console (static version - to be used outside of classes that derive from O...
Definition: logger.h:49
bool operator==(const Matrix &mat) const
Equality operator.
Definition: matrix.h:104
Vector2 operator*(const Vector2 &v) const
Matrix-vector multiplication.
Definition: matrix.h:521
const T & operator()(int i, int j) const
Indexing operator (const verions)
Definition: matrix.h:101
Matrix2x2(Stream *stream)
Unserialize a matrix from a stream.
Definition: matrix.h:475
T det(const TVector2< T > &v1, const TVector2< T > &v2)
Definition: vector.h:411
Vector4 col(int i) const
Return a column by index.
Definition: matrix.h:736
Float det() const
Return the determinant (Faster than Matrix::det)
Definition: matrix.h:487
#define MTS_NAMESPACE_BEGIN
Definition: platform.h:137
void writeArray(const T *array, size_t count)
Write an array to the stream (uses partial template specialization to select a method appropriate to ...
Definition: stream.h:520
bool isZero() const
Check if the matrix is identically zero.
Definition: matrix.h:365
Matrix3x3(const Float _m[9])
Initialize the matrix from a given (float) 3x3 array in row-major order.
Definition: matrix.h:571
Matrix(const T _m[M *N])
Initialize the matrix from a given (flat) MxN array in row-major order.
Definition: matrix.h:70
Matrix4x4()
Definition: matrix.h:657
Matrix3x3(const Matrix< 3, 3, Float > &mtx)
Copy constructor.
Definition: matrix.h:584
Matrix4x4 & operator=(const Matrix< 4, 4, Float > &mat)
Assignment operator.
Definition: matrix.h:721
Matrix4x4(const Float _m[16])
Initialize the matrix from a given (float) 4x4 array in row-major order.
Definition: matrix.h:666
Matrix3x3(const Vector &v1, const Vector &v2, const Vector &v3)
Initialize the matrix from three 3D column vectors.
Definition: matrix.h:574
Vector4 operator*(const Vector4 &v) const
Matrix-vector multiplication.
Definition: matrix.h:702
Matrix()
Construct a new MxN matrix without initializing it.
Definition: matrix.h:48
Vector2 row(int i) const
Return a row by index.
Definition: matrix.h:546
Matrix(const T _m[M][N])
Initialize the matrix from a given MxN array.
Definition: matrix.h:65
void transpose(Matrix< N, M, T > &target) const
Compute the transpose of this matrix.
Definition: matrix.h:414
Matrix(Stream *stream)
Unserialize a matrix from a stream.
Definition: matrix.h:75
Matrix2x2 operator*(Float value) const
Scalar multiplication (creates a temporary)
Definition: matrix.h:529
Matrix4x4(const Float _m[4][4])
Initialize the matrix from a given 4x4 array.
Definition: matrix.h:663
TVector2< Float > Vector2
Definition: fwd.h:106
const Matrix & operator-=(const Matrix &mat)
Matrix subtraction.
Definition: matrix.h:182
TVector4< Float > Vector4
Definition: fwd.h:122
T det() const
Compute the determinant of a square matrix (internally creates a LU decomposition) ...
Definition: matrix.h:388
Abstract seekable stream class.
Definition: stream.h:58
Matrix3x3()
Definition: matrix.h:562
T m[M][N]
Definition: matrix.h:35
Warning message.
Definition: formatter.h:32
const Matrix & operator*=(const Matrix &mat)
Matrix multiplication (for square matrices)
Definition: matrix.h:259
const Matrix & operator-=(T value)
Matrix-scalar addition.
Definition: matrix.h:198
Matrix2x2()
Definition: matrix.h:457
Definition: fwd.h:96
TVector3< Float > Vector3
Definition: fwd.h:115
const Matrix & operator+=(T value)
Matrix-scalar addition.
Definition: matrix.h:156
const Matrix & operator/=(T value)
Scalar division.
Definition: matrix.h:246
Basic 4x4 matrix data type.
Definition: matrix.h:656
Matrix2x2(Float a00, Float a01, Float a10, Float a11)
Initialize with the given values.
Definition: matrix.h:481
Matrix2x2(const Vector2 &v1, const Vector2 &v2)
Initialize the matrix from two 2D column vectors.
Definition: matrix.h:469
Matrix operator-() const
Component-wise negation.
Definition: matrix.h:206
const Matrix & operator*=(T value)
Scalar multiplication.
Definition: matrix.h:224
Matrix2x2(const Float _m[2][2])
Initialize the matrix from a given 2x2 array.
Definition: matrix.h:463
Float det3x3() const
Return the determinant of the upper left 3x3 sub-matrix.
Definition: matrix.h:695
Matrix operator/(T value) const
Scalar division (creates a temporary)
Definition: matrix.h:232
void readArray(T *array, size_t count)
Read an array from the stream (uses partial template specialization to select a method appropriate to...
Definition: stream.h:516
Matrix operator+(T value) const
Matrix-scalar addition (returns a temporary)
Definition: matrix.h:139
Matrix3x3(Float value)
Initialize the matrix with constant entries.
Definition: matrix.h:565
const Matrix & operator-(T value)
Matrix-scalar subtraction.
Definition: matrix.h:190
Matrix operator-(T value) const
Matrix-scalar subtraction (returns a temporary)
Definition: matrix.h:173
FINLINE bool invert2x2(Matrix2x2 &target, Float &det) const
Compute the inverse with det (Faster than Matrix::invert)
Definition: matrix.h:508
Matrix & operator=(const Matrix &mat)
Assignment operator.
Definition: matrix.h:122
#define RCPOVERFLOW
Definition: constants.h:97
Matrix(const Matrix &mtx)
Copy constructor.
Definition: matrix.h:80
Float det() const
Return the determinant (Faster than Matrix::det())
Definition: matrix.h:596
const Matrix & operator+=(const Matrix &mat)
Matrix addition.
Definition: matrix.h:148
Definition: fwd.h:97
void setIdentity()
Initialize with the identity matrix.
Definition: matrix.h:85
Float frob() const
Compute the Frobenius norm.
Definition: matrix.h:276
Definition: fwd.h:95
#define MTS_NAMESPACE_END
Definition: platform.h:138
bool isIdentity() const
Test if this is the identity matrix.
Definition: matrix.h:374
Matrix4x4(const Matrix< 4, 4, Float > &mtx)
Copy constructor.
Definition: matrix.h:680
Basic 3x3 matrix data type.
Definition: matrix.h:560
Matrix2x2 & operator=(const Matrix< 2, 2, Float > &mat)
Assignment operator.
Definition: matrix.h:538
T luDet(int pivsign) const
Compute the determinant of a decomposed matrix created by lu()
Definition: matrix.inl:122
void symEig(Matrix &Q, T d[M]) const
Perform a symmetric eigendecomposition of a square matrix into Q and D.
Definition: matrix.h:405
Vector2 col(int i) const
Return a column by index.
Definition: matrix.h:551
bool eig3(Matrix3x3 &m, Float lambda[3])
Fast 3x3 eigenvalue decomposition.
void setZero()
Initialize with zeroes.
Definition: matrix.h:92