Hermes
Loading...
Searching...
No Matches
matrix.h
Go to the documentation of this file.
1
31
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
39namespace hermes {
40
41// *********************************************************************************************************************
42// AUXILIARY FUNCTIONS
43// *********************************************************************************************************************
50template<typename T>
51HERMES_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// *********************************************************************************************************************
121template<typename T> class Matrix4x4 : public MathElement<T, 16> {
122public:
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
299private:
300 T m_[4][4];
301};
302
303// *********************************************************************************************************************
304// Matrix3x3
305// *********************************************************************************************************************
308template<typename T> class Matrix3x3 : public MathElement<T, 9> {
309public:
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
409private:
410 T m_[3][3];
411};
412
413// *********************************************************************************************************************
414// Matrix2x2
415// *********************************************************************************************************************
418template<typename T> class Matrix2x2 : public MathElement<T, 4> {
419public:
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
473private:
474 T m_[2][2];
475};
476
477// *********************************************************************************************************************
478// EXTERNAL FUNCTIONS
479// *********************************************************************************************************************
480template<typename T>
481HERMES_DEVICE_CALLABLE
482Matrix4x4<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}
487template<typename T>
488HERMES_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}
494template<typename T>
495HERMES_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}
617template<typename T>
618HERMES_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
642template<typename T>
643HERMES_DEVICE_CALLABLE Matrix2x2<T> operator*(T f, const Matrix2x2<T> &m) {
644 return m * f;
645}
646// algebra
647template<typename T>
648HERMES_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
661template<typename T>
662HERMES_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
666template<typename T>
667HERMES_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}
687template<typename T>
688HERMES_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}
692template<typename T>
693HERMES_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
697template<typename T>
698HERMES_DEVICE_CALLABLE Matrix3x3<T> operator*(T f, const Matrix3x3<T> &m) {
699 return m * f;
700}
701template<typename T>
702HERMES_DEVICE_CALLABLE Matrix4x4<T> operator*(T f, const Matrix4x4<T> &m) {
703 return m * f;
704}
705// *********************************************************************************************************************
706// IO
707// *********************************************************************************************************************
708template<typename T>
709std::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}
717template<typename T>
718std::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}
726template<typename T>
727std::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// *********************************************************************************************************************
739typedef Matrix4x4<real_t> mat4;
740typedef Matrix3x3<real_t> mat3;
741typedef 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::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:59
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 HERMES_DEVICE_CALLABLE constexpr bool is_equal(T a, T b)
Checks if two numbers are at most 1e-8 apart.
Definition numeric.h:847
hermes::Index2
Holds 2-dimensional integer index coordinates.
Definition index.h:50
vector.h
Geometric vector classes.