Hermes
matrix.h
Go to the documentation of this file.
1 
32 #ifndef HERMES_GEOMETRY_MATRIX_H
33 #define HERMES_GEOMETRY_MATRIX_H
34 
35 #include <hermes/geometry/vector.h>
36 #include <vector>
37 #include <hermes/numeric/math_element.h>
38 
39 namespace hermes {
40 
41 // *********************************************************************************************************************
42 // AUXILIARY FUNCTIONS
43 // *********************************************************************************************************************
50 template<typename T>
51 HERMES_DEVICE_CALLABLE bool gluInvertMatrix(const T m[16], T invOut[16]) {
52  T inv[16], det;
53  int i;
54 
55  inv[0] = m[5] * m[10] * m[15] - m[5] * m[11] * m[14] - m[9] * m[6] * m[15] +
56  m[9] * m[7] * m[14] + m[13] * m[6] * m[11] - m[13] * m[7] * m[10];
57 
58  inv[4] = -m[4] * m[10] * m[15] + m[4] * m[11] * m[14] + m[8] * m[6] * m[15] -
59  m[8] * m[7] * m[14] - m[12] * m[6] * m[11] + m[12] * m[7] * m[10];
60 
61  inv[8] = m[4] * m[9] * m[15] - m[4] * m[11] * m[13] - m[8] * m[5] * m[15] +
62  m[8] * m[7] * m[13] + m[12] * m[5] * m[11] - m[12] * m[7] * m[9];
63 
64  inv[12] = -m[4] * m[9] * m[14] + m[4] * m[10] * m[13] + m[8] * m[5] * m[14] -
65  m[8] * m[6] * m[13] - m[12] * m[5] * m[10] + m[12] * m[6] * m[9];
66 
67  inv[1] = -m[1] * m[10] * m[15] + m[1] * m[11] * m[14] + m[9] * m[2] * m[15] -
68  m[9] * m[3] * m[14] - m[13] * m[2] * m[11] + m[13] * m[3] * m[10];
69 
70  inv[5] = m[0] * m[10] * m[15] - m[0] * m[11] * m[14] - m[8] * m[2] * m[15] +
71  m[8] * m[3] * m[14] + m[12] * m[2] * m[11] - m[12] * m[3] * m[10];
72 
73  inv[9] = -m[0] * m[9] * m[15] + m[0] * m[11] * m[13] + m[8] * m[1] * m[15] -
74  m[8] * m[3] * m[13] - m[12] * m[1] * m[11] + m[12] * m[3] * m[9];
75 
76  inv[13] = m[0] * m[9] * m[14] - m[0] * m[10] * m[13] - m[8] * m[1] * m[14] +
77  m[8] * m[2] * m[13] + m[12] * m[1] * m[10] - m[12] * m[2] * m[9];
78 
79  inv[2] = m[1] * m[6] * m[15] - m[1] * m[7] * m[14] - m[5] * m[2] * m[15] +
80  m[5] * m[3] * m[14] + m[13] * m[2] * m[7] - m[13] * m[3] * m[6];
81 
82  inv[6] = -m[0] * m[6] * m[15] + m[0] * m[7] * m[14] + m[4] * m[2] * m[15] -
83  m[4] * m[3] * m[14] - m[12] * m[2] * m[7] + m[12] * m[3] * m[6];
84 
85  inv[10] = m[0] * m[5] * m[15] - m[0] * m[7] * m[13] - m[4] * m[1] * m[15] +
86  m[4] * m[3] * m[13] + m[12] * m[1] * m[7] - m[12] * m[3] * m[5];
87 
88  inv[14] = -m[0] * m[5] * m[14] + m[0] * m[6] * m[13] + m[4] * m[1] * m[14] -
89  m[4] * m[2] * m[13] - m[12] * m[1] * m[6] + m[12] * m[2] * m[5];
90 
91  inv[3] = -m[1] * m[6] * m[11] + m[1] * m[7] * m[10] + m[5] * m[2] * m[11] -
92  m[5] * m[3] * m[10] - m[9] * m[2] * m[7] + m[9] * m[3] * m[6];
93 
94  inv[7] = m[0] * m[6] * m[11] - m[0] * m[7] * m[10] - m[4] * m[2] * m[11] +
95  m[4] * m[3] * m[10] + m[8] * m[2] * m[7] - m[8] * m[3] * m[6];
96 
97  inv[11] = -m[0] * m[5] * m[11] + m[0] * m[7] * m[9] + m[4] * m[1] * m[11] -
98  m[4] * m[3] * m[9] - m[8] * m[1] * m[7] + m[8] * m[3] * m[5];
99 
100  inv[15] = m[0] * m[5] * m[10] - m[0] * m[6] * m[9] - m[4] * m[1] * m[10] +
101  m[4] * m[2] * m[9] + m[8] * m[1] * m[6] - m[8] * m[2] * m[5];
102 
103  det = m[0] * inv[0] + m[1] * inv[4] + m[2] * inv[8] + m[3] * inv[12];
104 
105  if (det == 0)
106  return false;
107 
108  det = 1.0f / det;
109 
110  for (i = 0; i < 16; i++)
111  invOut[i] = inv[i] * det;
112 
113  return true;
114 }
115 
116 // *********************************************************************************************************************
117 // Matrix4x4
118 // *********************************************************************************************************************
121 template<typename T> class Matrix4x4 : public MathElement<T, 16> {
122 public:
123  // *******************************************************************************************************************
124  // STATIC
125  // *******************************************************************************************************************
128  HERMES_DEVICE_CALLABLE static inline Matrix4x4 I() {
129  return {1, 0, 0, 0,
130  0, 1, 0, 0,
131  0, 0, 1, 0,
132  0, 0, 0, 1};
133  }
134  // *******************************************************************************************************************
135  // CONSTRUCTORS
136  // *******************************************************************************************************************
139  HERMES_DEVICE_CALLABLE explicit Matrix4x4(bool is_identity = false) {
140  std::memset(m_, 0, sizeof(m_));
141  if (is_identity)
142  for (int i = 0; i < 4; i++)
143  m_[i][i] = 1.f;
144  }
147  HERMES_DEVICE_CALLABLE Matrix4x4(std::initializer_list<T> values, bool columnMajor = false) {
148  size_t l = 0, c = 0;
149  for (auto v : values) {
150  m_[l][c] = v;
151  if (columnMajor) {
152  l++;
153  if (l >= 4)
154  l = 0, c++;
155  } else {
156  c++;
157  if (c >= 4)
158  c = 0, l++;
159  }
160  }
161  }
164  HERMES_DEVICE_CALLABLE explicit Matrix4x4(const T mat[16], bool columnMajor = false) {
165  size_t k = 0;
166  if (columnMajor)
167  for (int c = 0; c < 4; c++)
168  for (auto &l : m_)
169  l[c] = mat[k++];
170  else
171  for (auto &l : m_)
172  for (int c = 0; c < 4; c++)
173  l[c] = mat[k++];
174  }
176  HERMES_DEVICE_CALLABLE explicit Matrix4x4(T mat[4][4]) {
177  for (int i = 0; i < 4; i++)
178  for (int j = 0; j < 4; j++)
179  m_[i][j] = mat[i][j];
180  }
197  HERMES_DEVICE_CALLABLE Matrix4x4(T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20,
198  T m21, T m22, T m23, T m30, T m31, T m32, T m33) {
199  m_[0][0] = m00;
200  m_[0][1] = m01;
201  m_[0][2] = m02;
202  m_[0][3] = m03;
203  m_[1][0] = m10;
204  m_[1][1] = m11;
205  m_[1][2] = m12;
206  m_[1][3] = m13;
207  m_[2][0] = m20;
208  m_[2][1] = m21;
209  m_[2][2] = m22;
210  m_[2][3] = m23;
211  m_[3][0] = m30;
212  m_[3][1] = m31;
213  m_[3][2] = m32;
214  m_[3][3] = m33;
215  }
216  // *******************************************************************************************************************
217  // OPERATORS
218  // *******************************************************************************************************************
219  // arithmetic
220  HERMES_DEVICE_CALLABLE Matrix4x4<T> operator*(const Matrix4x4<T> &B) const {
221  Matrix4x4<T> r;
222  for (int i = 0; i < 4; ++i)
223  for (int j = 0; j < 4; ++j)
224  r[i][j] = m_[i][0] * B[0][j] + m_[i][1] * B[1][j] +
225  m_[i][2] * B[2][j] + m_[i][3] * B[3][j];
226  return r;
227  }
228  HERMES_DEVICE_CALLABLE Vector4 <T> operator*(const Vector4 <T> &v) const {
229  Vector4<T> r;
230  for (int i = 0; i < 4; i++)
231  for (int j = 0; j < 4; j++)
232  r[i] += m_[i][j] * v[j];
233  return r;
234  }
235  HERMES_DEVICE_CALLABLE Matrix4x4<T> &operator*=(T s) {
236  for (int i = 0; i < 4; i++)
237  for (int j = 0; j < 4; j++)
238  m_[i][j] *= s;
239  return *this;
240  }
241  // boolean
242  HERMES_DEVICE_CALLABLE bool operator==(const Matrix4x4<T> &B) const {
243  for (int i = 0; i < 4; i++)
244  for (int j = 0; j < 4; j++)
245  if (!Check::is_equal(m_[i][j], B[i][j]))
246  return false;
247  return true;
248  }
249  HERMES_DEVICE_CALLABLE bool operator!=(const Matrix4x4<T> &B) const {
250  return !((*this) == B);
251  }
252  // *******************************************************************************************************************
253  // METHODS
254  // *******************************************************************************************************************
256  HERMES_DEVICE_CALLABLE void setIdentity() {
257 #if defined(__CUDA_ARCH__) && __CUDA_ARCH__ > 0
258  for (auto & i : m_)
259  for (int j = 0; j < 4; j++)
260  i[j] = 0.f;
261 #else
262  std::memset(m_, 0, sizeof(m_));
263 #endif
264  for (int i = 0; i < 4; i++)
265  m_[i][i] = 1.f;
266  }
269  HERMES_DEVICE_CALLABLE void row_major(T *a) const {
270  int k = 0;
271  for (auto &i : m_)
272  for (int j = 0; j < 4; j++)
273  a[k++] = i[j];
274  }
277  HERMES_DEVICE_CALLABLE void column_major(T *a) const {
278  int k = 0;
279  for (int i = 0; i < 4; i++)
280  for (auto &j : m_)
281  a[k++] = j[i];
282  }
285  [[nodiscard]] HERMES_DEVICE_CALLABLE bool isIdentity() const {
286  for (int i = 0; i < 4; i++)
287  for (int j = 0; j < 4; j++)
288  if ((i != j && !Check::is_equal(m_[i][j], 0.f)) ||
289  (i == j && !Check::is_equal(m_[i][j], 1.f)))
290  return false;
291  return true;
292  }
293  // *******************************************************************************************************************
294  // ACCESS
295  // *******************************************************************************************************************
296  HERMES_DEVICE_CALLABLE T *operator[](u32 row_index) { return m_[row_index]; }
297  HERMES_DEVICE_CALLABLE const T *operator[](u32 row_index) const { return m_[row_index]; }
298 
299 private:
300  T m_[4][4];
301 };
302 
303 // *********************************************************************************************************************
304 // Matrix3x3
305 // *********************************************************************************************************************
308 template<typename T> class Matrix3x3 : public MathElement<T, 9> {
309 public:
310  // *******************************************************************************************************************
311  // STATIC METHODS
312  // *******************************************************************************************************************
313  HERMES_DEVICE_CALLABLE static inline Matrix3x3 I() {
314  return {1, 0, 0,
315  0, 1, 0,
316  0, 0, 1};
317  }
318  HERMES_DEVICE_CALLABLE static inline Matrix3x3 diag(T m00, T m11, T m22) {
319  return Matrix3x3(m00, 0, 0, //
320  0, m11, 0,//
321  0, 0, m22);
322  }
323  // *******************************************************************************************************************
324  // CONSTRUCTORS
325  // *******************************************************************************************************************
326  HERMES_DEVICE_CALLABLE Matrix3x3() { std::memset(m_, 0, sizeof(m_)); }
327  HERMES_DEVICE_CALLABLE Matrix3x3(vec3 a, vec3 b, vec3 c)
328  : Matrix3x3(a.x, a.y, a.z, b.x, b.y, b.z, c.x, c.y, c.z) {}
329  HERMES_DEVICE_CALLABLE Matrix3x3(T m00, T m01, T m02, T m10, T m11, T m12, T m20, T m21, T m22) {
330  m_[0][0] = m00;
331  m_[0][1] = m01;
332  m_[0][2] = m02;
333  m_[1][0] = m10;
334  m_[1][1] = m11;
335  m_[1][2] = m12;
336  m_[2][0] = m20;
337  m_[2][1] = m21;
338  m_[2][2] = m22;
339  }
340  // *******************************************************************************************************************
341  // OPERATORS
342  // *******************************************************************************************************************
343  // arithmetic
344  HERMES_DEVICE_CALLABLE
345  Vector3 <T> operator*(const Vector3 <T> &v) const {
346  Vector3<T> r;
347  for (int i = 0; i < 3; i++)
348  for (int j = 0; j < 3; j++)
349  r[i] += m_[i][j] * v[j];
350  return r;
351  }
352  HERMES_DEVICE_CALLABLE
353  Matrix3x3<T> operator*(const Matrix3x3<T> &B) const {
354  Matrix3x3<T> r;
355  for (int i = 0; i < 3; ++i)
356  for (int j = 0; j < 3; ++j)
357  r[i][j] =
358  m_[i][0] * B[0][j] + m_[i][1] * B[1][j] + m_[i][2] * B[2][j];
359  return r;
360  }
361 
362 #define ARITHMETIC_OP(OP) \
363  HERMES_DEVICE_CALLABLE Matrix3x3<T> operator OP##=(const Matrix3x3<T> &B) const { \
364  for (int i = 0; i < 3; ++i) \
365  for (int j = 0; j < 3; ++j) \
366  m_[i][j] OP##= B[i][j]; \
367  } \
368  HERMES_DEVICE_CALLABLE Matrix3x3<T> operator OP(const Matrix3x3<T> &B) const { \
369  Matrix3x3<T> r; \
370  for (int i = 0; i < 3; ++i) \
371  for (int j = 0; j < 3; ++j) \
372  r[i][j] = m_[i][j] OP B[i][j]; \
373  return r; }
374  ARITHMETIC_OP(+)
375  ARITHMETIC_OP(-)
376 #undef ARITHMETIC_OP
377 
378  HERMES_DEVICE_CALLABLE
379  Matrix3x3<T> operator*(const T &f) const {
380  return Matrix3x3<T>(m_[0][0] * f, m_[0][1] * f, m_[0][2] * f, m_[1][0] * f,
381  m_[1][1] * f, m_[1][2] * f, m_[2][0] * f, m_[2][1] * f,
382  m_[2][2] * f);
383  }
384  HERMES_DEVICE_CALLABLE
385  Matrix3x3<T> operator/(const T &f) const {
386  return Matrix3x3<T>(m_[0][0] / f, m_[0][1] / f, m_[0][2] / f, m_[1][0] / f,
387  m_[1][1] / f, m_[1][2] / f, m_[2][0] / f, m_[2][1] / f,
388  m_[2][2] / f);
389  }
390  // *******************************************************************************************************************
391  // METHODS
392  // *******************************************************************************************************************
393  HERMES_DEVICE_CALLABLE void setIdentity() {
394  std::memset(m_, 0, sizeof(m_));
395  for (int i = 0; i < 3; i++)
396  m_[i][i] = 1.f;
397  }
398  HERMES_DEVICE_CALLABLE T determinant() const {
399  return m_[0][0] * m_[1][1] * m_[2][2] + m_[0][1] * m_[1][2] * m_[2][0] +
400  m_[0][2] * m_[1][0] * m_[2][1] - m_[2][0] * m_[1][1] * m_[0][2] -
401  m_[2][1] * m_[1][2] * m_[0][0] - m_[2][2] * m_[1][0] * m_[0][1];
402  }
403  // *******************************************************************************************************************
404  // ACCESS
405  // *******************************************************************************************************************
406  HERMES_DEVICE_CALLABLE T *operator[](u32 row_index) { return m_[row_index]; }
407  HERMES_DEVICE_CALLABLE const T *operator[](u32 row_index) const { return m_[row_index]; }
408 
409 private:
410  T m_[3][3];
411 };
412 
413 // *********************************************************************************************************************
414 // Matrix2x2
415 // *********************************************************************************************************************
418 template<typename T> class Matrix2x2 : public MathElement<T, 4> {
419 public:
420  // *******************************************************************************************************************
421  // CONSTRUCTORS
422  // *******************************************************************************************************************
423  HERMES_DEVICE_CALLABLE Matrix2x2() {
424  std::memset(m_, 0, sizeof(m_));
425  for (int i = 0; i < 2; i++)
426  m_[i][i] = 1.f;
427  }
428  HERMES_DEVICE_CALLABLE Matrix2x2(T m00, T m01, T m10, T m11) {
429  m_[0][0] = m00;
430  m_[0][1] = m01;
431  m_[1][0] = m10;
432  m_[1][1] = m11;
433  }
434  // *******************************************************************************************************************
435  // OPERATORS
436  // *******************************************************************************************************************
437  // arithmetic
438  HERMES_DEVICE_CALLABLE Vector2 <T> operator*(const Vector2 <T> &v) const {
439  Vector2<T> r;
440  for (int i = 0; i < 2; i++)
441  for (int j = 0; j < 2; j++)
442  r[i] += m_[i][j] * v[j];
443  return r;
444  }
445  HERMES_DEVICE_CALLABLE Matrix2x2<T> operator*(const Matrix2x2<T> &B) const {
446  Matrix2x2 <T> r;
447  for (int i = 0; i < 2; ++i)
448  for (int j = 0; j < 2; ++j)
449  r[i][j] = m_[i][0] * B[0][j] + m_[i][1] * B[1][j];
450  return r;
451  }
452  HERMES_DEVICE_CALLABLE Matrix2x2<T> operator*(T f) const {
453  return Matrix2x2<T>(m_[0][0] * f, m_[0][1] * f, m_[1][0] * f,
454  m_[1][1] * f);
455  }
456  // *******************************************************************************************************************
457  // METHODS
458  // *******************************************************************************************************************
459  HERMES_DEVICE_CALLABLE void setIdentity() {
460  std::memset(m_, 0, sizeof(m_));
461  for (int i = 0; i < 2; i++)
462  m_[i][i] = 1.f;
463  }
464  HERMES_DEVICE_CALLABLE T determinant() { return m_[0][0] * m_[1][1] - m_[0][1] * m_[1][0]; }
465  // *******************************************************************************************************************
466  // ACCESS
467  // *******************************************************************************************************************
468  HERMES_DEVICE_CALLABLE const T &operator()(u32 i, u32 j) const { return m_[i][j]; }
469  HERMES_DEVICE_CALLABLE T &operator()(u32 i, u32 j) { return m_[i][j]; }
470  HERMES_DEVICE_CALLABLE T *operator[](u32 row_index) { return m_[row_index]; }
471  HERMES_DEVICE_CALLABLE const T *operator[](u32 row_index) const { return m_[row_index]; }
472 
473 private:
474  T m_[2][2];
475 };
476 
477 // *********************************************************************************************************************
478 // EXTERNAL FUNCTIONS
479 // *********************************************************************************************************************
480 template<typename T>
481 HERMES_DEVICE_CALLABLE
482 Matrix4x4<T> rowReduce(const Matrix4x4<T> &p, const Matrix4x4<T> &q) {
483  Matrix4x4<T> l = p, r = q;
484  // TODO implement with gauss jordan elimination
485  return r;
486 }
487 template<typename T>
488 HERMES_DEVICE_CALLABLE Matrix4x4<T> transpose(const Matrix4x4<T> &m) {
489  return Matrix4x4<T>(m[0][0], m[1][0], m[2][0], m[3][0], m[0][1],
490  m[1][1], m[2][1], m[3][1], m[0][2], m[1][2],
491  m[2][2], m[3][2], m[0][3], m[1][3], m[2][3],
492  m[3][3]);
493 }
494 template<typename T>
495 HERMES_DEVICE_CALLABLE Matrix4x4<T> inverse(const Matrix4x4<T> &m) {
496  Matrix4x4<T> r;
497  T mm[16], inv[16];
498  m.row_major(mm);
499  if (gluInvertMatrix(mm, inv)) {
500  int k = 0;
501  for (int i = 0; i < 4; i++)
502  for (int j = 0; j < 4; j++)
503  r[i][j] = inv[k++];
504  return r;
505  }
506 
507  T det = m[0][0] * m[1][1] * m[2][2] * m[3][3] +
508  m[1][2] * m[2][3] * m[3][1] * m[1][3] +
509  m[2][1] * m[3][2] * m[1][1] * m[2][3] +
510  m[3][2] * m[1][2] * m[2][1] * m[3][3] +
511  m[1][3] * m[2][2] * m[3][1] * m[0][1] +
512  m[0][1] * m[2][3] * m[3][2] * m[0][2] +
513  m[2][1] * m[3][3] * m[0][3] * m[2][2] +
514  m[3][1] * m[0][1] * m[2][2] * m[3][3] +
515  m[0][2] * m[2][3] * m[3][1] * m[0][3] +
516  m[2][1] * m[3][2] * m[0][2] * m[0][1] +
517  m[1][2] * m[3][3] * m[0][2] * m[1][3] +
518  m[3][1] * m[0][3] * m[1][1] * m[3][2] -
519  m[0][1] * m[1][3] * m[3][2] * m[0][2] -
520  m[1][1] * m[3][3] * m[0][3] * m[1][2] -
521  m[3][1] * m[0][3] * m[0][1] * m[1][3] -
522  m[2][2] * m[0][2] * m[1][1] * m[2][3] -
523  m[0][3] * m[1][2] * m[2][1] * m[0][1] -
524  m[1][2] * m[2][3] * m[0][2] * m[1][3] -
525  m[2][1] * m[0][3] * m[1][1] * m[2][2] -
526  m[1][0] * m[1][0] * m[2][3] * m[3][2] -
527  m[1][2] * m[2][0] * m[3][3] * m[1][3] -
528  m[2][2] * m[3][0] * m[1][0] * m[2][2] -
529  m[3][3] * m[1][2] * m[2][3] * m[3][0] -
530  m[1][3] * m[2][0] * m[3][2] * m[1][1];
531  if (fabs(det) < 1e-8)
532  return r;
533 
534  r[0][0] =
535  (m[1][1] * m[2][2] * m[3][3] + m[1][2] * m[2][3] * m[3][1] +
536  m[1][3] * m[2][1] * m[3][2] - m[1][1] * m[2][3] * m[3][2] -
537  m[1][2] * m[2][1] * m[3][3] - m[1][3] * m[2][2] * m[3][1]) /
538  det;
539  r[0][1] =
540  (m[0][1] * m[2][3] * m[3][2] + m[0][2] * m[2][1] * m[3][3] +
541  m[0][3] * m[2][2] * m[3][1] - m[0][1] * m[2][2] * m[3][3] -
542  m[0][2] * m[2][3] * m[3][1] - m[0][3] * m[2][1] * m[3][2]) /
543  det;
544  r[0][2] =
545  (m[0][1] * m[1][2] * m[3][3] + m[0][2] * m[1][3] * m[3][1] +
546  m[0][3] * m[1][1] * m[3][2] - m[0][1] * m[1][3] * m[3][2] -
547  m[0][2] * m[1][1] * m[3][3] - m[0][3] * m[1][2] * m[3][1]) /
548  det;
549  r[0][3] =
550  (m[0][1] * m[1][3] * m[2][2] + m[0][2] * m[1][1] * m[2][3] +
551  m[0][3] * m[1][2] * m[2][1] - m[0][1] * m[1][2] * m[2][3] -
552  m[0][2] * m[1][3] * m[2][1] - m[0][3] * m[1][1] * m[2][2]) /
553  det;
554  r[1][0] =
555  (m[1][0] * m[2][3] * m[3][2] + m[1][2] * m[2][0] * m[3][3] +
556  m[1][3] * m[2][2] * m[3][0] - m[1][0] * m[2][2] * m[3][3] -
557  m[1][2] * m[2][3] * m[3][0] - m[1][3] * m[2][0] * m[3][2]) /
558  det;
559  r[1][1] =
560  (m[0][0] * m[2][2] * m[3][3] + m[0][2] * m[2][3] * m[3][0] +
561  m[0][3] * m[2][0] * m[3][2] - m[0][0] * m[2][3] * m[3][2] -
562  m[0][2] * m[2][0] * m[3][3] - m[0][3] * m[2][2] * m[3][0]) /
563  det;
564  r[1][2] =
565  (m[0][0] * m[1][3] * m[3][2] + m[0][2] * m[1][0] * m[3][3] +
566  m[0][3] * m[1][2] * m[3][0] - m[0][0] * m[1][2] * m[3][3] -
567  m[0][2] * m[1][3] * m[3][0] - m[0][3] * m[1][0] * m[3][2]) /
568  det;
569  r[1][3] =
570  (m[0][0] * m[1][2] * m[2][3] + m[0][2] * m[1][3] * m[2][0] +
571  m[0][3] * m[1][0] * m[2][2] - m[0][0] * m[1][3] * m[2][2] -
572  m[0][2] * m[1][0] * m[2][3] - m[0][3] * m[1][2] * m[2][0]) /
573  det;
574  r[2][0] =
575  (m[1][0] * m[2][1] * m[3][3] + m[1][1] * m[2][3] * m[3][0] +
576  m[1][3] * m[2][0] * m[3][1] - m[1][0] * m[2][3] * m[3][1] -
577  m[1][1] * m[2][0] * m[3][3] - m[1][3] * m[2][1] * m[3][0]) /
578  det;
579  r[2][1] =
580  (m[0][0] * m[2][3] * m[3][1] + m[0][1] * m[2][0] * m[3][3] +
581  m[0][3] * m[2][1] * m[3][0] - m[0][0] * m[2][1] * m[3][3] -
582  m[0][1] * m[2][3] * m[3][0] - m[0][3] * m[2][0] * m[3][1]) /
583  det;
584  r[2][2] =
585  (m[0][0] * m[1][1] * m[3][3] + m[0][1] * m[1][3] * m[3][0] +
586  m[0][3] * m[1][0] * m[3][1] - m[0][0] * m[1][3] * m[3][1] -
587  m[0][1] * m[1][0] * m[3][3] - m[0][3] * m[1][1] * m[3][0]) /
588  det;
589  r[2][3] =
590  (m[0][0] * m[1][3] * m[2][1] + m[0][1] * m[1][0] * m[2][3] +
591  m[0][3] * m[1][1] * m[2][0] - m[0][0] * m[1][1] * m[2][3] -
592  m[0][1] * m[1][3] * m[2][0] - m[0][3] * m[1][0] * m[2][1]) /
593  det;
594  r[3][0] =
595  (m[1][0] * m[2][2] * m[3][1] + m[1][1] * m[2][0] * m[3][2] +
596  m[1][2] * m[2][1] * m[3][0] - m[1][0] * m[2][1] * m[3][2] -
597  m[1][1] * m[2][2] * m[3][0] - m[1][2] * m[2][0] * m[3][1]) /
598  det;
599  r[3][1] =
600  (m[0][0] * m[2][1] * m[3][2] + m[0][1] * m[2][2] * m[3][0] +
601  m[0][2] * m[2][0] * m[3][1] - m[0][0] * m[2][2] * m[3][1] -
602  m[0][1] * m[2][0] * m[3][2] - m[0][2] * m[2][1] * m[3][0]) /
603  det;
604  r[3][2] =
605  (m[0][0] * m[1][2] * m[3][1] + m[0][1] * m[1][0] * m[3][2] +
606  m[0][2] * m[1][1] * m[3][0] - m[0][0] * m[1][1] * m[3][2] -
607  m[0][1] * m[1][2] * m[3][0] - m[0][2] * m[1][0] * m[3][1]) /
608  det;
609  r[3][3] =
610  (m[0][0] * m[1][1] * m[2][2] + m[0][1] * m[1][2] * m[2][0] +
611  m[0][2] * m[1][0] * m[2][1] - m[0][0] * m[1][2] * m[2][1] -
612  m[0][1] * m[1][0] * m[2][2] - m[0][2] * m[1][1] * m[2][0]) /
613  det;
614 
615  return r;
616 }
617 template<typename T>
618 HERMES_DEVICE_CALLABLE void decompose(const Matrix4x4<T> &m, Matrix4x4<T> &r, Matrix4x4<T> &s) {
619  // extract rotation r from transformation matrix
620  T norm;
621  int count = 0;
622  r = m;
623  do {
624  // compute next matrix in series
625  Matrix4x4<T> Rnext;
626  Matrix4x4<T> Rit = inverse(transpose(r));
627  for (int i = 0; i < 4; i++)
628  for (int j = 0; j < 4; j++)
629  Rnext[i][j] = .5f * (r[i][j] + Rit[i][j]);
630  // compute norm difference between R and Rnext
631  norm = 0.f;
632  for (int i = 0; i < 3; i++) {
633  T n = fabsf(r[i][0] - Rnext[i][0]) +
634  fabsf(r[i][1] - Rnext[i][1]) + fabsf(r[i][2] - Rnext[i][2]);
635  norm = std::max(norm, n);
636  }
637  } while (++count < 100 && norm > .0001f);
638  // compute scale S using rotation and original matrix
639  s = inverse(r) * m;
640 }
641 // arithmetic
642 template<typename T>
643 HERMES_DEVICE_CALLABLE Matrix2x2<T> operator*(T f, const Matrix2x2<T> &m) {
644  return m * f;
645 }
646 // algebra
647 template<typename T>
648 HERMES_DEVICE_CALLABLE Matrix2x2<T> inverse(const Matrix2x2<T> &m) {
649  Matrix2x2 <T> r;
650  T det = m[0][0] * m[1][1] - m[0][1] * m[1][0];
651  if (det == 0.f)
652  return r;
653  T k = 1.f / det;
654  r[0][0] = m[1][1] * k;
655  r[0][1] = -m[0][1] * k;
656  r[1][0] = -m[1][0] * k;
657  r[1][1] = m[0][0] * k;
658  return r;
659 }
660 
661 template<typename T>
662 HERMES_DEVICE_CALLABLE Matrix2x2<T> transpose(const Matrix2x2<T> &m) {
663  return Matrix2x2<T>(m[0][0], m[1][0], m[0][1], m[1][1]);
664 }
665 // algebra
666 template<typename T>
667 HERMES_DEVICE_CALLABLE Matrix3x3<T> inverse(const Matrix3x3<T> &m) {
668  Matrix3x3<T> r;
669  // r.setIdentity();
670  T det =
671  m[0][0] * m[1][1] * m[2][2] + m[1][0] * m[2][1] * m[0][2] +
672  m[2][0] * m[0][1] * m[1][2] - m[0][0] * m[2][1] * m[1][2] -
673  m[2][0] * m[1][1] * m[0][2] - m[1][0] * m[0][1] * m[2][2];
674  if (std::fabs(det) < 1e-8)
675  return r;
676  r[0][0] = (m[1][1] * m[2][2] - m[1][2] * m[2][1]) / det;
677  r[0][1] = (m[0][2] * m[2][1] - m[0][1] * m[2][2]) / det;
678  r[0][2] = (m[0][1] * m[1][2] - m[0][2] * m[1][1]) / det;
679  r[1][0] = (m[1][2] * m[2][0] - m[1][0] * m[2][2]) / det;
680  r[1][1] = (m[0][0] * m[2][2] - m[0][2] * m[2][0]) / det;
681  r[1][2] = (m[0][2] * m[1][0] - m[0][0] * m[1][2]) / det;
682  r[2][0] = (m[1][0] * m[2][1] - m[1][1] * m[2][0]) / det;
683  r[2][1] = (m[0][1] * m[2][0] - m[0][0] * m[2][1]) / det;
684  r[2][2] = (m[0][0] * m[1][1] - m[0][1] * m[1][0]) / det;
685  return r;
686 }
687 template<typename T>
688 HERMES_DEVICE_CALLABLE Matrix3x3<T> transpose(const Matrix3x3<T> &m) {
689  return Matrix3x3<T>(m[0][0], m[1][0], m[2][0], m[0][1], m[1][1],
690  m[2][1], m[0][2], m[1][2], m[2][2]);
691 }
692 template<typename T>
693 HERMES_DEVICE_CALLABLE Matrix3x3<T> star(const Vector3 <T> a) {
694  return Matrix3x3<T>(0, -a[2], a[1], a[2], 0, -a[0], -a[1], a[0], 0);
695 }
696 // arithmetic
697 template<typename T>
698 HERMES_DEVICE_CALLABLE Matrix3x3<T> operator*(T f, const Matrix3x3<T> &m) {
699  return m * f;
700 }
701 template<typename T>
702 HERMES_DEVICE_CALLABLE Matrix4x4<T> operator*(T f, const Matrix4x4<T> &m) {
703  return m * f;
704 }
705 // *********************************************************************************************************************
706 // IO
707 // *********************************************************************************************************************
708 template<typename T>
709 std::ostream &operator<<(std::ostream &os, const Matrix4x4<T> &m) {
710  for (int i = 0; i < 4; i++) {
711  for (int j = 0; j < 4; j++)
712  os << m[i][j] << " ";
713  os << std::endl;
714  }
715  return os;
716 }
717 template<typename T>
718 std::ostream &operator<<(std::ostream &os, const Matrix3x3<T> &m) {
719  for (int i = 0; i < 3; i++) {
720  for (int j = 0; j < 3; j++)
721  os << m[i][j] << " ";
722  os << std::endl;
723  }
724  return os;
725 }
726 template<typename T>
727 std::ostream &operator<<(std::ostream &os, const Matrix2x2<T> &m) {
728  for (int i = 0; i < 2; i++) {
729  for (int j = 0; j < 2; j++)
730  os << m[i][j] << " ";
731  os << std::endl;
732  }
733  return os;
734 }
735 
736 // *********************************************************************************************************************
737 // TYPEDEFS
738 // *********************************************************************************************************************
739 typedef Matrix4x4<real_t> mat4;
740 typedef Matrix3x3<real_t> mat3;
741 typedef Matrix2x2<real_t> mat2;
742 
743 } // namespace hermes
744 
745 #endif
746 
hermes::MathElement
Interface used by all basic geometric entities.
Definition: math_element.h:44
hermes::Matrix2x2
2x2 Matrix representation
Definition: matrix.h:418
hermes::Matrix3x3
3x3 Matrix representation
Definition: matrix.h:308
hermes::Matrix4x4
4x4 Matrix representation
Definition: matrix.h:121
hermes::Matrix4x4::isIdentity
HERMES_DEVICE_CALLABLE bool isIdentity() const
Definition: matrix.h:285
hermes::Matrix4x4::Matrix4x4
HERMES_DEVICE_CALLABLE Matrix4x4(std::initializer_list< T > values, bool columnMajor=false)
Definition: matrix.h:147
hermes::Matrix4x4::column_major
HERMES_DEVICE_CALLABLE void column_major(T *a) const
Definition: matrix.h:277
hermes::Matrix4x4::Matrix4x4
HERMES_DEVICE_CALLABLE Matrix4x4(T mat[4][4])
Definition: matrix.h:176
hermes::Matrix4x4::Matrix4x4
HERMES_DEVICE_CALLABLE Matrix4x4(const T mat[16], bool columnMajor=false)
Definition: matrix.h:164
hermes::Matrix4x4::Matrix4x4
HERMES_DEVICE_CALLABLE Matrix4x4(bool is_identity=false)
Definition: matrix.h:139
hermes::Matrix4x4::Matrix4x4
HERMES_DEVICE_CALLABLE Matrix4x4(T m00, T m01, T m02, T m03, T m10, T m11, T m12, T m13, T m20, T m21, T m22, T m23, T m30, T m31, T m32, T m33)
Definition: matrix.h:197
hermes::Matrix4x4::row_major
HERMES_DEVICE_CALLABLE void row_major(T *a) const
Definition: matrix.h:269
hermes::Matrix4x4::I
static HERMES_DEVICE_CALLABLE Matrix4x4 I()
Creates an identity matrix.
Definition: matrix.h:128
hermes::Vector2
Geometric 2-dimensional vector (x, y)
Definition: vector.h:54
hermes::Vector3::x
T x
0-th component
Definition: vector.h:415
hermes::Vector3::z
T z
2-th component
Definition: vector.h:417
hermes::Vector3::y
T y
1-th component
Definition: vector.h:416
hermes::cuda_utils::operator<<
std::ostream & operator<<(std::ostream &o, const LaunchInfo &info)
LaunchInfo support for std::ostream << operator.
Definition: cuda_utils.h:234
HERMES_DEVICE_CALLABLE
#define HERMES_DEVICE_CALLABLE
Specifies that the function can be called from both host and device sides.
Definition: defs.h:45
u32
uint32_t u32
32 bit size unsigned integer type
Definition: defs.h:88
ARITHMETIC_OP
#define ARITHMETIC_OP(OP)
asd
Definition: index.h:408
math_element.h
Base class for all geometric objects.
hermes::gluInvertMatrix
HERMES_DEVICE_CALLABLE bool gluInvertMatrix(const T m[16], T invOut[16])
Inverts a 4x4 matrix.
Definition: matrix.h:51
hermes::Check::is_equal
static constexpr HERMES_DEVICE_CALLABLE bool is_equal(T a, T b)
Checks if two numbers are at most 1e-8 apart.
Definition: numeric.h:847
vector.h
Geometric vector classes.