• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
1 /*
2  * Copyright (C) 2023 Huawei Device Co., Ltd.
3  * Licensed under the Apache License, Version 2.0 (the "License");
4  * you may not use this file except in compliance with the License.
5  * You may obtain a copy of the License at
6  *
7  *     http://www.apache.org/licenses/LICENSE-2.0
8  *
9  * Unless required by applicable law or agreed to in writing, software
10  * distributed under the License is distributed on an "AS IS" BASIS,
11  * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
12  * See the License for the specific language governing permissions and
13  * limitations under the License.
14  */
15 
16 #ifndef API_BASE_MATH_MATRIX_UTIL_H
17 #define API_BASE_MATH_MATRIX_UTIL_H
18 
19 #include <base/math/matrix.h>
20 #include <base/math/quaternion.h>
21 #include <base/math/vector_util.h>
22 #include <base/namespace.h>
23 #include <base/util/log.h>
24 
BASE_BEGIN_NAMESPACE()25 BASE_BEGIN_NAMESPACE()
26 namespace Math {
27 /** \addtogroup group_math_matrixutil
28  *  @{
29  */
30 /** Returns transpose of this matrix */
31 static inline constexpr Mat3X3 Transpose(const Mat3X3& m)
32 {
33     Mat3X3 result;
34 
35     result.x.x = m.x.x;
36     result.y.x = m.x.y;
37     result.z.x = m.x.z;
38 
39     result.x.y = m.y.x;
40     result.y.y = m.y.y;
41     result.z.y = m.y.z;
42 
43     result.x.z = m.z.x;
44     result.y.z = m.z.y;
45     result.z.z = m.z.z;
46 
47     return result;
48 }
49 
50 /** Returns transpose of this matrix */
51 static inline constexpr Mat4X4 Transpose(const Mat4X4& m)
52 {
53     Mat4X4 result;
54 
55     result.x.x = m.x.x;
56     result.y.x = m.x.y;
57     result.z.x = m.x.z;
58     result.w.x = m.x.w;
59 
60     result.x.y = m.y.x;
61     result.y.y = m.y.y;
62     result.z.y = m.y.z;
63     result.w.y = m.y.w;
64 
65     result.x.z = m.z.x;
66     result.y.z = m.z.y;
67     result.z.z = m.z.z;
68     result.w.z = m.z.w;
69 
70     result.x.w = m.w.x;
71     result.y.w = m.w.y;
72     result.z.w = m.w.z;
73     result.w.w = m.w.w;
74 
75     return result;
76 }
77 
78 /** Set scale of this matrix */
79 static inline constexpr Mat3X3 PostScale(const Mat3X3& mat, const Vec2& vec)
80 {
81     Mat3X3 result;
82     result.x = { mat.x.x * vec.x, mat.x.y * vec.y, mat.x.z };
83     result.y = { mat.y.x * vec.x, mat.y.y * vec.y, mat.y.z };
84     result.z = { mat.z.x * vec.x, mat.z.y * vec.y, mat.z.z };
85     return result;
86 }
87 
88 /** Set scale of this matrix */
89 static inline constexpr Mat4X4 PostScale(const Mat4X4& mat, const Vec3& vec)
90 {
91     Mat4X4 result;
92     result.x = { mat.x.x * vec.x, mat.x.y * vec.y, mat.x.z * vec.z, mat.x.w };
93     result.y = { mat.y.x * vec.x, mat.y.y * vec.y, mat.y.z * vec.z, mat.y.w };
94     result.z = { mat.z.x * vec.x, mat.z.y * vec.y, mat.z.z * vec.z, mat.z.w };
95     result.w = { mat.w.x * vec.x, mat.w.y * vec.y, mat.w.z * vec.z, mat.w.w };
96     return result;
97 }
98 
99 /** Set scale of this matrix */
100 static inline constexpr Mat3X3 Scale(const Mat3X3& mat, const Vec2& vec)
101 {
102     Mat3X3 result;
103     result.x = mat.x * vec.x;
104     result.y = mat.y * vec.y;
105     result.z = mat.z;
106     return result;
107 }
108 
109 /** Set scale of this matrix */
110 static inline constexpr Mat4X4 Scale(const Mat4X4& mat, const Vec3& vec)
111 {
112     Mat4X4 result;
113     result.x = mat.x * vec.x;
114     result.y = mat.y * vec.y;
115     result.z = mat.z * vec.z;
116     result.w = mat.w;
117     return result;
118 }
119 
120 /** Get column of the matrix */
121 static inline constexpr Vec3 GetColumn(const Mat3X3& mat, int index)
122 {
123     switch (index) {
124             // the first column of the matrix: Vec3(m00, m10, m20)
125         case 0:
126             return mat.x;
127             // the second column of the matrix:  Vec3(m01, m11, m21)
128         case 1:
129             return mat.y;
130             // the third column of the matrix:  Vec3(m02, m12, m22)
131         case 2:
132             return mat.z;
133 
134         default:
135             BASE_ASSERT(false);
136             return Vec3(0.0f, 0.0f, 0.0f);
137     }
138 }
139 
140 /** Get column of the matrix */
141 static inline constexpr Vec4 GetColumn(const Mat4X4& mat, int index)
142 {
143     switch (index) {
144             // the first column of the matrix: Vec4(m00, m10, m20, m30)
145         case 0:
146             return mat.x;
147             // the second column of the matrix:  Vec4(m01, m11, m21, m31)
148         case 1:
149             return mat.y;
150             // the third column of the matrix:  Vec4(m02, m12, m22, m32)
151         case 2:
152             return mat.z;
153             // the fourth column of the matrix:  Vec4(m03, m13, m23, m33)
154         case 3:
155             return mat.w;
156 
157         default:
158             BASE_ASSERT(false);
159             return Vec4(0.0f, 0.0f, 0.0f, 0.0f);
160     }
161 }
162 
163 /** Returns row of the matrix */
164 static inline constexpr Vec3 GetRow(const Mat3X3& mat, int index)
165 {
166     switch (index) {
167             // the first row of the matrix: Vec4(m00, m01, m02, m03)
168         case 0:
169             return Vec3(mat.x.x, mat.y.x, mat.z.x);
170             // the second row of the matrix: Vec4(m10, m11, m12, m13)
171         case 1:
172             return Vec3(mat.x.y, mat.y.y, mat.z.y);
173             // the third row of the matrix: Vec4(m20, m21, m22, m23)
174         case 2:
175             return Vec3(mat.x.z, mat.y.z, mat.z.z);
176         default:
177             BASE_ASSERT(false);
178             return Vec3(0.0f, 0.0f, 0.0f);
179     }
180 }
181 
182 /** Returns row of the matrix */
183 static inline constexpr Vec4 GetRow(const Mat4X4& mat, int index)
184 {
185     switch (index) {
186             // the first row of the matrix: Vec4(m00, m01, m02, m03)
187         case 0:
188             return Vec4(mat.x.x, mat.y.x, mat.z.x, mat.w.x);
189             // the second row of the matrix: Vec4(m10, m11, m12, m13)
190         case 1:
191             return Vec4(mat.x.y, mat.y.y, mat.z.y, mat.w.y);
192             // the third row of the matrix: Vec4(m20, m21, m22, m23)
193         case 2:
194             return Vec4(mat.x.z, mat.y.z, mat.z.z, mat.w.z);
195             // the fourth row of the matrix: Vec4(m30, m31, m32, m33)
196         case 3:
197             return Vec4(mat.x.w, mat.y.w, mat.z.w, mat.w.w);
198 
199         default:
200             BASE_ASSERT(false);
201             return Vec4(0.0f, 0.0f, 0.0f, 0.0f);
202     }
203 }
204 
205 /** Converts quaternion to matrix 3x3 */
206 static inline constexpr Mat3X3 Mat3Cast(const Quat& quaternion)
207 {
208     Mat3X3 result;
209     const float qXX(quaternion.x * quaternion.x);
210     const float qYY(quaternion.y * quaternion.y);
211     const float qZZ(quaternion.z * quaternion.z);
212     const float qXZ(quaternion.x * quaternion.z);
213     const float qXY(quaternion.x * quaternion.y);
214     const float qYZ(quaternion.y * quaternion.z);
215     const float qWX(quaternion.w * quaternion.x);
216     const float qWY(quaternion.w * quaternion.y);
217     const float qWZ(quaternion.w * quaternion.z);
218 
219     result.x.x = 1.0f - 2.0f * (qYY + qZZ);
220     result.x.y = 2.0f * (qXY + qWZ);
221     result.x.z = 2.0f * (qXZ - qWY);
222 
223     result.y.x = 2.0f * (qXY - qWZ);
224     result.y.y = 1.0f - 2.0f * (qXX + qZZ);
225     result.y.z = 2.0f * (qYZ + qWX);
226 
227     result.z.x = 2.0f * (qXZ + qWY);
228     result.z.y = 2.0f * (qYZ - qWX);
229     result.z.z = 1.0f - 2.0f * (qXX + qYY);
230     return result;
231 }
232 
233 /** Convert 3x3 matrix to 4x4 */
234 static inline constexpr Mat4X4 DimensionalShift(const Mat3X3& mat)
235 {
236     Mat4X4 result;
237     result.x = { mat.x.x, mat.x.y, 0, mat.x.z };
238     result.y = { mat.y.x, mat.y.y, 0, mat.y.z };
239     result.z = { 0, 0, 1, 0 };
240     result.w = { mat.z.x, mat.z.y, 0, mat.z.z };
241     return result;
242 }
243 
244 /** Converts quaternion to matrix 4X4 */
245 static inline constexpr Mat4X4 Mat4Cast(const Quat& quaternion)
246 {
247     return Mat4X4(Mat3Cast(quaternion));
248 }
249 
250 /** Translate by vector2 */
251 static inline constexpr Mat3X3 Translate(const Mat3X3& mat, const Vec2& vec)
252 {
253     Mat3X3 result(mat);
254     result.z = mat.x * vec.x + mat.y * vec.y + mat.z;
255     return result;
256 }
257 
258 /** Translate by vector2 */
259 static inline constexpr Mat4X4 Translate(const Mat4X4& mat, const Vec2& vec)
260 {
261     Mat4X4 result(mat);
262     result.w = mat.x * vec.x + mat.y * vec.y + mat.w;
263     return result;
264 }
265 /** Translate by vector3 */
266 static inline constexpr Mat4X4 Translate(const Mat4X4& mat, const Vec3& vec)
267 {
268     Mat4X4 result(mat);
269     result.w = mat.x * vec.x + mat.y * vec.y + mat.z * vec.z + mat.w;
270     return result;
271 }
272 
273 /** Creates translation, rotation and scaling matrix from translation vector, rotation quaternion and scale vector */
274 static inline constexpr Mat4X4 Trs(
275     const Vec3& translationVector, const Quat& rotationQuaternion, const Vec3& scaleVector)
276 {
277     return Scale(Translate(Mat4X4(1.f), translationVector) * Mat4Cast(rotationQuaternion), scaleVector);
278 }
279 
280 /** Transforms direction by this matrix */
281 static inline constexpr Vec3 MultiplyVector(const Mat4X4& mat, const Vec3& vec)
282 {
283     Vec3 res;
284     res.x = mat.x.x * vec.x + mat.y.x * vec.y + mat.z.x * vec.z;
285     res.y = mat.x.y * vec.x + mat.y.y * vec.y + mat.z.y * vec.z;
286     res.z = mat.x.z * vec.x + mat.y.z * vec.y + mat.z.z * vec.z;
287     return res;
288 }
289 
290 /** Transforms position by this matrix without perspective divide */
291 static inline constexpr Vec2 MultiplyPoint2X4(const Mat4X4& mat4X4, const Vec2& point)
292 {
293     Vec2 result;
294     result.x = mat4X4.x.x * point.x + mat4X4.y.x * point.y + mat4X4.w.x;
295     result.y = mat4X4.x.y * point.x + mat4X4.y.y * point.y + mat4X4.w.y;
296 
297     return result;
298 }
299 
300 /** Transforms position by this matrix without perspective divide */
301 static inline constexpr Vec2 MultiplyPoint2X3(const Mat3X3& mat3X3, const Vec2& point)
302 {
303     Vec2 result;
304     result.x = mat3X3.x.x * point.x + mat3X3.y.x * point.y + mat3X3.z.x;
305     result.y = mat3X3.x.y * point.x + mat3X3.y.y * point.y + mat3X3.z.y;
306 
307     return result;
308 }
309 
310 /** Transforms position by this matrix without perspective divide */
311 static inline constexpr Vec3 MultiplyPoint3X4(const Mat4X4& mat4X4, const Vec3& point)
312 {
313     Vec3 result;
314     result.x = mat4X4.x.x * point.x + mat4X4.y.x * point.y + mat4X4.z.x * point.z + mat4X4.w.x;
315     result.y = mat4X4.x.y * point.x + mat4X4.y.y * point.y + mat4X4.z.y * point.z + mat4X4.w.y;
316     result.z = mat4X4.x.z * point.x + mat4X4.y.z * point.y + mat4X4.z.z * point.z + mat4X4.w.z;
317 
318     return result;
319 }
320 
321 /** Inverse matrix with determinant output */
322 static inline constexpr Mat3X3 Inverse(Mat3X3 const& m, float& determinantOut)
323 {
324     determinantOut = m.x.x * (m.y.y * m.z.z - m.z.y * m.y.z) - m.x.y * (m.y.x * m.z.z - m.y.z * m.z.x) +
325                      m.x.z * (m.y.x * m.z.y - m.y.y * m.z.x);
326 
327     const float invdet = 1 / determinantOut;
328 
329     Mat3X3 inverse;
330     inverse.x.x = (m.y.y * m.z.z - m.z.y * m.y.z) * invdet;
331     inverse.x.y = (m.x.z * m.z.y - m.x.y * m.z.z) * invdet;
332     inverse.x.z = (m.x.y * m.y.z - m.x.z * m.y.y) * invdet;
333     inverse.y.x = (m.y.z * m.z.x - m.y.x * m.z.z) * invdet;
334     inverse.y.y = (m.x.x * m.z.z - m.x.z * m.z.x) * invdet;
335     inverse.y.z = (m.y.x * m.x.z - m.x.x * m.y.z) * invdet;
336     inverse.z.x = (m.y.x * m.z.y - m.z.x * m.y.y) * invdet;
337     inverse.z.y = (m.z.x * m.x.y - m.x.x * m.z.y) * invdet;
338     inverse.z.z = (m.x.x * m.y.y - m.y.x * m.x.y) * invdet;
339     return inverse;
340 }
341 
342 /** Inverse matrix */
343 static inline constexpr Mat3X3 Inverse(const Mat3X3& mat3X3)
344 {
345     float determinant = 0.0f;
346     return Inverse(mat3X3, determinant);
347 }
348 
349 /** Inverse matrix with determinant output */
350 static inline constexpr Mat4X4 Inverse(Mat4X4 const& m, float& determinantOut)
351 {
352     const float coef00 = m.z.z * m.w.w - m.w.z * m.z.w;
353     const float coef02 = m.y.z * m.w.w - m.w.z * m.y.w;
354     const float coef03 = m.y.z * m.z.w - m.z.z * m.y.w;
355 
356     const float coef04 = m.z.y * m.w.w - m.w.y * m.z.w;
357     const float coef06 = m.y.y * m.w.w - m.w.y * m.y.w;
358     const float coef07 = m.y.y * m.z.w - m.z.y * m.y.w;
359 
360     const float coef08 = m.z.y * m.w.z - m.w.y * m.z.z;
361     const float coef10 = m.y.y * m.w.z - m.w.y * m.y.z;
362     const float coef11 = m.y.y * m.z.z - m.z.y * m.y.z;
363 
364     const float coef12 = m.z.x * m.w.w - m.w.x * m.z.w;
365     const float coef14 = m.y.x * m.w.w - m.w.x * m.y.w;
366     const float coef15 = m.y.x * m.z.w - m.z.x * m.y.w;
367 
368     const float coef16 = m.z.x * m.w.z - m.w.x * m.z.z;
369     const float coef18 = m.y.x * m.w.z - m.w.x * m.y.z;
370     const float coef19 = m.y.x * m.z.z - m.z.x * m.y.z;
371 
372     const float coef20 = m.z.x * m.w.y - m.w.x * m.z.y;
373     const float coef22 = m.y.x * m.w.y - m.w.x * m.y.y;
374 
375     const float coef23 = m.y.x * m.z.y - m.z.x * m.y.y;
376 
377     const Vec4 fac0(coef00, coef00, coef02, coef03);
378     const Vec4 fac1(coef04, coef04, coef06, coef07);
379     const Vec4 fac2(coef08, coef08, coef10, coef11);
380     const Vec4 fac3(coef12, coef12, coef14, coef15);
381     const Vec4 fac4(coef16, coef16, coef18, coef19);
382     const Vec4 fac5(coef20, coef20, coef22, coef23);
383 
384     const Vec4 vec0(m.y.x, m.x.x, m.x.x, m.x.x);
385     const Vec4 vec1(m.y.y, m.x.y, m.x.y, m.x.y);
386     const Vec4 vec2(m.y.z, m.x.z, m.x.z, m.x.z);
387     const Vec4 vec3(m.y.w, m.x.w, m.x.w, m.x.w);
388 
389     const Vec4 inv0(vec1 * fac0 - vec2 * fac1 + vec3 * fac2);
390     const Vec4 inv1(vec0 * fac0 - vec2 * fac3 + vec3 * fac4);
391     const Vec4 inv2(vec0 * fac1 - vec1 * fac3 + vec3 * fac5);
392     const Vec4 inv3(vec0 * fac2 - vec1 * fac4 + vec2 * fac5);
393 
394     constexpr Vec4 signA(+1.f, -1.f, +1.f, -1.f);
395     constexpr Vec4 signB(-1.f, +1.f, -1.f, +1.f);
396     const Mat4X4 inverse(inv0 * signA, inv1 * signB, inv2 * signA, inv3 * signB);
397 
398     const Vec4 row0(inverse.x.x, inverse.y.x, inverse.z.x, inverse.w.x);
399 
400     const Vec4 dot0(m.x * row0);
401     determinantOut = (dot0.x + dot0.y) + (dot0.z + dot0.w);
402 
403     const float oneOverDeterminant = 1.0f / determinantOut;
404 
405     return inverse * oneOverDeterminant;
406 }
407 
408 /** Get determinant out from matrix 4X4 */
409 static inline constexpr float Determinant(const Mat4X4& mat4X4)
410 {
411     float determinant = 0.0f;
412     Inverse(mat4X4, determinant);
413     return determinant;
414 }
415 
416 /** Inverse matrix */
417 static inline constexpr Mat4X4 Inverse(const Mat4X4& mat4X4)
418 {
419     float determinant = 0.0f;
420     return Inverse(mat4X4, determinant);
421 }
422 
423 /** Inner product */
424 inline constexpr Vec4 operator*(const Mat4X4& m, const Vec4& v)
425 {
426     return { v.x * m.x.x + v.y * m.y.x + v.z * m.z.x + v.w * m.w.x,
427         v.x * m.x.y + v.y * m.y.y + v.z * m.z.y + v.w * m.w.y, v.x * m.x.z + v.y * m.y.z + v.z * m.z.z + v.w * m.w.z,
428         v.x * m.x.w + v.y * m.y.w + v.z * m.z.w + v.w * m.w.w };
429 }
430 
431 inline constexpr Vec3 operator*(const Mat3X3& m, const Vec3& v)
432 {
433     return { v.x * m.x.x + v.y * m.y.x + v.z * m.z.x, v.x * m.x.y + v.y * m.y.y + v.z * m.z.y,
434         v.x * m.x.z + v.y * m.y.z + v.z * m.z.z };
435 }
436 
437 /** Creates matrix for left handed symmetric perspective view frustum, near and far clip planes correspond to z
438  * normalized device coordinates of 0 and +1 */
439 static inline Mat4X4 PerspectiveLhZo(float fovy, float aspect, float zNear, float zFar)
440 {
441     BASE_ASSERT(abs(aspect) > EPSILON);
442 
443     float const tanHalfFovy = tan(fovy / 2.0f);
444 
445     Mat4X4 result(0.f);
446     if (const auto div = (aspect * tanHalfFovy); abs(div) > EPSILON) {
447         result.x.x = 1.0f / div;
448     } else {
449         result.x.x = HUGE_VALF;
450     }
451     if (abs(tanHalfFovy) > EPSILON) {
452         result.y.y = 1.0f / (tanHalfFovy);
453     } else {
454         result.y.y = HUGE_VALF;
455     }
456 
457     if (const auto div = (zFar - zNear); abs(div) > EPSILON) {
458         result.z.z = zFar / div;
459         result.w.z = -(zFar * zNear) / div;
460     } else {
461         result.z.z = HUGE_VALF;
462         result.w.z = HUGE_VALF;
463     }
464     result.z.w = 1.0f;
465     return result;
466 }
467 
468 /** Creates matrix for right handed symmetric perspective view frustum, near and far clip planes correspond to z
469  * normalized device coordinates of 0 and +1 */
470 static inline Mat4X4 PerspectiveRhZo(float fovy, float aspect, float zNear, float zFar)
471 {
472     Mat4X4 result = PerspectiveLhZo(fovy, aspect, zNear, zFar);
473     result.z.z = -result.z.z;
474     result.z.w = -result.z.w;
475 
476     return result;
477 }
478 
479 /** Creates matrix for left handed symmetric perspective view frustum, near and far clip planes correspond to z
480  * normalized device coordinates of -1 and +1 */
481 static inline Mat4X4 PerspectiveLhNo(float fovy, float aspect, float zNear, float zFar)
482 {
483     BASE_ASSERT(abs(aspect) > EPSILON);
484 
485     float const tanHalfFovy = tan(fovy / 2.0f);
486 
487     Mat4X4 result(0.f);
488     if (const auto div = (aspect * tanHalfFovy); abs(div) > EPSILON) {
489         result.x.x = 1.0f / div;
490     } else {
491         result.x.x = HUGE_VALF;
492     }
493     if (abs(tanHalfFovy) > EPSILON) {
494         result.y.y = 1.0f / (tanHalfFovy);
495     } else {
496         result.y.y = HUGE_VALF;
497     }
498 
499     if (const auto div = (zFar - zNear); abs(div) > EPSILON) {
500         result.z.z = (zFar + zNear) / div;
501         result.w.z = -(2.0f * zFar * zNear) / div;
502     } else {
503         result.z.z = HUGE_VALF;
504         result.w.z = HUGE_VALF;
505     }
506     result.z.w = 1.0f;
507     return result;
508 }
509 
510 /** Creates a matrix for right handed symmetric perspective view frustum, near and far clip planes correspond to z
511  * normalized device coordinates of -1 and +1 */
512 static inline Mat4X4 PerspectiveRhNo(float fovy, float aspect, float zNear, float zFar)
513 {
514     Mat4X4 result = PerspectiveLhNo(fovy, aspect, zNear, zFar);
515     result.z.z = -result.z.z;
516     result.z.w = -result.z.w;
517 
518     return result;
519 }
520 
521 /** Creates matrix for orthographic parallel viewing volume using left handed coordinates, near and far clip planes
522  * correspond to z normalized device coordinates of 0 and +1 */
523 static inline Mat4X4 OrthoLhZo(float left, float right, float bottom, float top, float zNear, float zFar)
524 {
525     Mat4X4 result(1.f);
526     result.x.x = 2.0f / (right - left);
527     result.y.y = 2.0f / (top - bottom);
528     result.z.z = 1.0f / (zFar - zNear);
529     result.w.x = -(right + left) / (right - left);
530     result.w.y = -(top + bottom) / (top - bottom);
531     result.w.z = -zNear / (zFar - zNear);
532     return result;
533 }
534 
535 /** Creates matrix for orthographic parallel viewing volume using left handed coordinates, near and far clip planes
536  * correspond to z normalized device coordinates of -1 and +1 */
537 static inline Mat4X4 OrthoLhNo(float left, float right, float bottom, float top, float zNear, float zFar)
538 {
539     Mat4X4 result(1.f);
540     result.x.x = 2.0f / (right - left);
541     result.y.y = 2.0f / (top - bottom);
542     result.z.z = 2.0f / (zFar - zNear);
543     result.w.x = -(right + left) / (right - left);
544     result.w.y = -(top + bottom) / (top - bottom);
545     result.w.z = -(zFar + zNear) / (zFar - zNear);
546     return result;
547 }
548 
549 /** Creates matrix for orthographic parallel viewing volume using right handed coordinates, near and far clip planes
550  * correspond to z normalized device coordinates of 0 and +1 */
551 static inline Mat4X4 OrthoRhZo(float left, float right, float bottom, float top, float zNear, float zFar)
552 {
553     Mat4X4 result = OrthoLhZo(left, right, bottom, top, zNear, zFar);
554     result.z.z = -result.z.z;
555     return result;
556 }
557 
558 /** Creates matrix for orthographic parallel viewing volume using right handed coordinates, near and far clip planes
559  * correspond to z normalized device coordinates of -1 and +1 */
560 static inline Mat4X4 OrthoRhNo(float left, float right, float bottom, float top, float zNear, float zFar)
561 {
562     Mat4X4 result = OrthoLhNo(left, right, bottom, top, zNear, zFar);
563     result.z.z = -result.z.z;
564     return result;
565 }
566 
567 /** Right handed look at function */
568 static inline Mat4X4 LookAtRh(Vec3 const& eye, Vec3 const& center, Vec3 const& up)
569 {
570     Vec3 const f(Normalize(center - eye));
571     Vec3 const s(Normalize(Cross(f, up)));
572     Vec3 const u(Cross(s, f));
573 
574     Mat4X4 result(1.f);
575     result.x.x = s.x;
576     result.y.x = s.y;
577     result.z.x = s.z;
578     result.x.y = u.x;
579     result.y.y = u.y;
580     result.z.y = u.z;
581     result.x.z = -f.x;
582     result.y.z = -f.y;
583     result.z.z = -f.z;
584     result.w.x = -Dot(s, eye);
585     result.w.y = -Dot(u, eye);
586     result.w.z = Dot(f, eye);
587     return result;
588 }
589 
590 /** Left handed look at function */
591 static inline Mat4X4 LookAtLh(Vec3 const& eye, Vec3 const& center, Vec3 const& up)
592 {
593     Vec3 const f(Normalize(center - eye));
594     Vec3 const s(Normalize(Cross(up, f)));
595     Vec3 const u(Cross(f, s));
596 
597     Mat4X4 result(1.f);
598     result.x.x = s.x;
599     result.y.x = s.y;
600     result.z.x = s.z;
601     result.x.y = u.x;
602     result.y.y = u.y;
603     result.z.y = u.z;
604     result.x.z = f.x;
605     result.y.z = f.y;
606     result.z.z = f.z;
607     result.w.x = -Dot(s, eye);
608     result.w.y = -Dot(u, eye);
609     result.w.z = -Dot(f, eye);
610     return result;
611 }
612 
613 /** Decompose matrix */
614 static inline bool Decompose(
615     Mat4X4 const& modelMatrix, Vec3& scale, Quat& orientation, Vec3& translation, Vec3& skew, Vec4& perspective)
616 {
617     Mat4X4 localMatrix(modelMatrix);
618 
619     if (localMatrix.w.w != 1.f) {
620         if (abs(localMatrix.w.w) < EPSILON) {
621             return false;
622         }
623 
624         for (int i = 0; i < 4; ++i) {
625             for (int j = 0; j < 4; ++j) {
626                 localMatrix[i][j] /= localMatrix.w.w;
627             }
628         }
629     }
630 
631     if (abs(localMatrix.x.w >= EPSILON) || abs(localMatrix.y.w >= EPSILON) || abs(localMatrix.z.w >= EPSILON)) {
632         Vec4 rightHandSide;
633         rightHandSide.x = localMatrix.x.w;
634         rightHandSide.y = localMatrix.y.w;
635         rightHandSide.z = localMatrix.z.w;
636         rightHandSide.w = localMatrix.w.w;
637 
638         Mat4X4 perspectiveMatrix(localMatrix);
639         for (int i = 0; i < 3; i++) {
640             perspectiveMatrix[i].w = 0.0f;
641         }
642         perspectiveMatrix.w.w = 1.0f;
643 
644         float determinant;
645         const Mat4X4 inversePerspectiveMatrix = Inverse(perspectiveMatrix, determinant);
646         if (abs(determinant) < EPSILON) {
647             return false;
648         }
649 
650         const Mat4X4 transposedInversePerspectiveMatrix = Transpose(inversePerspectiveMatrix);
651 
652         perspective = transposedInversePerspectiveMatrix * rightHandSide;
653 
654         localMatrix.x.w = localMatrix.y.w = localMatrix.z.w = 0.0f;
655         localMatrix.w.w = 1.0f;
656     } else {
657         perspective = Vec4(0.0f, 0.0f, 0.0f, 1.0f);
658     }
659 
660     translation = Vec3(localMatrix.w);
661     localMatrix.w = Vec4(0.0f, 0.0f, 0.0f, localMatrix.w.w);
662 
663     Vec3 row[3];
664     for (int i = 0; i < 3; ++i) {
665         for (int j = 0; j < 3; ++j) {
666             row[i][j] = localMatrix[i][j];
667         }
668     }
669 
670     scale.x = sqrt(Dot(row[0], row[0]));
671 
672     row[0] = Math::Scale(row[0], 1.0f);
673 
674     skew.z = Dot(row[0], row[1]);
675     row[1] = Combine(row[1], row[0], 1.0f, -skew.z);
676 
677     scale.y = sqrt(Dot(row[1], row[1]));
678     row[1] = Math::Scale(row[1], 1.0f);
679     skew.z /= scale.y;
680 
681     skew.y = Dot(row[0], row[2]);
682     row[2] = Combine(row[2], row[0], 1.0f, -skew.y);
683     skew.x = Dot(row[1], row[2]);
684     row[2] = Combine(row[2], row[1], 1.0f, -skew.x);
685 
686     scale.z = sqrt(Dot(row[2], row[2]));
687     row[2] = Math::Scale(row[2], 1.0f);
688     skew.y /= scale.z;
689     skew.x /= scale.z;
690 
691     const Vec3 pDum3 = Cross(row[1], row[2]);
692     if (Dot(row[0], pDum3) < 0) {
693         for (int i = 0; i < 3; i++) {
694             scale[i] *= -1.0f;
695             row[i] *= -1.0f;
696         }
697     }
698 
699     int i, j, k = 0;
700     float root;
701     const float trace = row[0].x + row[1].y + row[2].z;
702     if (trace > 0.0f) {
703         root = sqrt(trace + 1.0f);
704         orientation.w = 0.5f * root;
705         root = 0.5f / root; // root cannot be zero as it's square root of at least 1
706         orientation.x = root * (row[1].z - row[2].y);
707         orientation.y = root * (row[2].x - row[0].z);
708         orientation.z = root * (row[0].y - row[1].x);
709     } else { // End if > 0
710         constexpr const int next[3] = { 1, 2, 0 };
711         i = 0;
712         if (row[1].y > row[0].x) {
713             i = 1;
714         }
715         if (row[2].z > row[i][i]) {
716             i = 2;
717         }
718         j = next[i];
719         k = next[j];
720 
721         root = sqrt(row[i][i] - row[j][j] - row[k][k] + 1.0f);
722 
723         orientation[i] = 0.5f * root;
724         root = (abs(root) > EPSILON) ? 0.5f / root : HUGE_VALF;
725         orientation[j] = root * (row[i][j] + row[j][i]);
726         orientation[k] = root * (row[i][k] + row[k][i]);
727         orientation.w = root * (row[j][k] - row[k][j]);
728     } // End if <= 0
729 
730     return true;
731 }
732 
733 /** Check if matrix has any rotation/shear components */
734 static inline bool HasRotation(Mat3X3 const& m)
735 {
736     return (m.x.y != 0.f || m.y.x != 0.f);
737 }
738 
739 /** Check if matrix has any rotation/shear components */
740 static inline bool HasRotation(Mat4X4 const& m)
741 {
742     return (m.x.y != 0.f || m.x.z != 0.f || m.y.x != 0.f || m.y.z != 0.f || m.z.x != 0.f || m.z.y != 0.f);
743 }
744 
745 /** @} */
746 } // namespace Math
747 BASE_END_NAMESPACE()
748 
749 #endif // API_BASE_MATH_MATRIX_UTIL_H
750