openMSX
gl_mat.hh
Go to the documentation of this file.
1#ifndef GL_MAT_HH
2#define GL_MAT_HH
3
4// This code implements a mathematical matrix, comparable in functionality
5// and syntax to the matrix types in GLSL.
6//
7// The code was tuned for 4x4 matrixes, but when it didn't complicate the
8// code, other sizes (even non-square) are supported as well.
9//
10// The functionality of this code is focused on openGL(ES). So only the
11// operation 'matrix x column-vector' is supported and not e.g. 'row-vector x
12// matrix'. The internal layout is also compatible with openGL(ES).
13//
14// In the past we had (some) manual SSE optimizations in this code. Though for
15// the functions that matter (matrix-vector and matrix-matrix multiplication),
16// the compiler's auto-vectorization has become as good as the manually
17// vectorized code.
18
19#include "gl_vec.hh"
20#include "xrange.hh"
21#include <cassert>
22
23namespace gl {
24
25// Matrix with M rows and N columns (M by N), elements have type T.
26// Internally elements are stored in column-major order. This is the same
27// convention as openGL (and Fortran), but different from the typical row-major
28// order in C++. This allows to directly pass these matrices to openGL(ES).
29template<int M, int N, typename T> class matMxN
30{
31public:
32 // Default copy-constructor and assignment operator.
33
34 // Construct identity matrix.
35 constexpr matMxN()
36 {
37 for (auto i : xrange(N)) {
38 for (auto j : xrange(M)) {
39 c[i][j] = (i == j) ? T(1) : T(0);
40 }
41 }
42 }
43
44 // Construct diagonal matrix.
45 constexpr explicit matMxN(const vecN<(M < N ? M : N), T>& d)
46 {
47 for (auto i : xrange(N)) {
48 for (auto j : xrange(M)) {
49 c[i][j] = (i == j) ? d[i] : T(0);
50 }
51 }
52 }
53
54 // Construct from larger matrix (higher order elements are dropped).
55 template<int M2, int N2> constexpr explicit matMxN(const matMxN<M2, N2, T>& x)
56 {
57 static_assert(((M2 > M) && (N2 >= N)) || ((M2 >= M) && (N2 > N)),
58 "matrix must have strictly larger dimensions");
59 for (auto i : xrange(N)) c[i] = vecN<M, T>(x[i]);
60 }
61
62 // Construct matrix from 2 given columns (only valid when N == 2).
63 constexpr matMxN(const vecN<M, T>& x, const vecN<M, T>& y)
64 {
65 static_assert(N == 2, "wrong #constructor arguments");
66 c[0] = x; c[1] = y;
67 }
68
69 // Construct matrix from 3 given columns (only valid when N == 3).
70 constexpr matMxN(const vecN<M, T>& x, const vecN<M, T>& y, const vecN<M, T>& z)
71 {
72 static_assert(N == 3, "wrong #constructor arguments");
73 c[0] = x; c[1] = y; c[2] = z;
74 }
75
76 // Construct matrix from 4 given columns (only valid when N == 4).
77 constexpr matMxN(const vecN<M, T>& x, const vecN<M, T>& y,
78 const vecN<M, T>& z, const vecN<M, T>& w)
79 {
80 static_assert(N == 4, "wrong #constructor arguments");
81 c[0] = x; c[1] = y; c[2] = z; c[3] = w;
82 }
83
84 // Access the i-the column of the matrix.
85 // Vectors are also indexable, so 'A[i][j]' returns the element
86 // at the i-th column, j-th row.
87 [[nodiscard]] constexpr const vecN<M, T>& operator[](unsigned i) const {
88 #ifdef DEBUG
89 assert(i < N);
90 #endif
91 return c[i];
92 }
93 [[nodiscard]] constexpr vecN<M, T>& operator[](unsigned i) {
94 #ifdef DEBUG
95 assert(i < N);
96 #endif
97 return c[i];
98 }
99
100 // Assignment version of the +,-,* operations defined below.
101 constexpr matMxN& operator+=(const matMxN& x) { *this = *this + x; return *this; }
102 constexpr matMxN& operator-=(const matMxN& x) { *this = *this - x; return *this; }
103 constexpr matMxN& operator*=(T x) { *this = *this * x; return *this; }
104 constexpr matMxN& operator*=(const matMxN<N, N, T>& x) { *this = *this * x; return *this; }
105
106 // gcc-10 mis-compiles this (fixed in gcc-11):
107 // [[nodiscard]] constexpr bool operator==(const matMxN&) const = default;
108 // For now still manually implement it.
109 [[nodiscard]] friend constexpr bool operator==(const matMxN& x, const matMxN& y)
110 {
111 for (auto i : xrange(N)) if (x[i] != y[i]) return false;
112 return true;
113 }
114
115private:
116 std::array<vecN<M, T>, N> c;
117};
118
119
120// Convenience typedefs (same names as used by GLSL).
124
125
126// -- Matrix functions --
127
128// matrix + matrix
129template<int M, int N, typename T>
130[[nodiscard]] constexpr matMxN<M, N, T> operator+(const matMxN<M, N, T>& A, const matMxN<M, N, T>& B)
131{
133 for (auto i : xrange(N)) R[i] = A[i] + B[i];
134 return R;
135}
136
137// matrix - matrix
138template<int M, int N, typename T>
139[[nodiscard]] constexpr matMxN<M, N, T> operator-(const matMxN<M, N, T>& A, const matMxN<M, N, T>& B)
140{
142 for (auto i : xrange(N)) R[i] = A[i] - B[i];
143 return R;
144}
145
146// matrix negation
147template<int M, int N, typename T>
148[[nodiscard]] constexpr matMxN<M, N, T> operator-(const matMxN<M, N, T>& A)
149{
150 return matMxN<M, N, T>(vecN<(M < N ? M : N), T>()) - A;
151}
152
153// scalar * matrix
154template<int M, int N, typename T>
155[[nodiscard]] constexpr matMxN<M, N, T> operator*(T x, const matMxN<M, N, T>& A)
156{
158 for (auto i : xrange(N)) R[i] = x * A[i];
159 return R;
160}
161
162// matrix * scalar
163template<int M, int N, typename T>
164[[nodiscard]] constexpr matMxN<M, N, T> operator*(const matMxN<M, N, T>& A, T x)
165{
167 for (auto i : xrange(N)) R[i] = A[i] * x;
168 return R;
169}
170
171// matrix * column-vector
172template<int M, int N, typename T>
173[[nodiscard]] constexpr vecN<M, T> operator*(const matMxN<M, N, T>& A, const vecN<N, T>& x)
174{
175 vecN<M, T> r;
176 for (auto i : xrange(N)) r += A[i] * x[i];
177 return r;
178}
179
180// matrix * matrix
181template<int M, int N, int O, typename T>
182[[nodiscard]] constexpr matMxN<M, O, T> operator*(const matMxN<M, N, T>& A, const matMxN<N, O, T>& B)
183{
185 for (auto i : xrange(O)) R[i] = A * B[i];
186 return R;
187}
188
189// matrix transpose
190template<int M, int N, typename T>
191[[nodiscard]] constexpr matMxN<N, M, T> transpose(const matMxN<M, N, T>& A)
192{
194 for (auto i : xrange(N)) {
195 for (auto j : xrange(M)) {
196 R[j][i] = A[i][j];
197 }
198 }
199 return R;
200}
201
202// determinant of a 2x2 matrix
203template<typename T>
204[[nodiscard]] constexpr T determinant(const matMxN<2, 2, T>& A)
205{
206 return A[0][0] * A[1][1] - A[0][1] * A[1][0];
207}
208
209// determinant of a 3x3 matrix
210template<typename T>
211[[nodiscard]] constexpr T determinant(const matMxN<3, 3, T>& A)
212{
213 return A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1])
214 - A[1][0] * (A[0][1] * A[2][2] - A[0][2] * A[2][1])
215 + A[2][0] * (A[0][1] * A[1][2] - A[0][2] * A[1][1]);
216}
217
218// determinant of a 4x4 matrix
219template<typename T>
220[[nodiscard]] constexpr T determinant(const matMxN<4, 4, T>& A)
221{
222 // Implementation based on glm: http://glm.g-truc.net
223 T f0 = A[2][2] * A[3][3] - A[3][2] * A[2][3];
224 T f1 = A[2][1] * A[3][3] - A[3][1] * A[2][3];
225 T f2 = A[2][1] * A[3][2] - A[3][1] * A[2][2];
226 T f3 = A[2][0] * A[3][3] - A[3][0] * A[2][3];
227 T f4 = A[2][0] * A[3][2] - A[3][0] * A[2][2];
228 T f5 = A[2][0] * A[3][1] - A[3][0] * A[2][1];
229 vecN<4, T> c((A[1][1] * f0 - A[1][2] * f1 + A[1][3] * f2),
230 (A[1][2] * f3 - A[1][3] * f4 - A[1][0] * f0),
231 (A[1][0] * f1 - A[1][1] * f3 + A[1][3] * f5),
232 (A[1][1] * f4 - A[1][2] * f5 - A[1][0] * f2));
233 return dot(A[0], c);
234}
235
236// inverse of a 2x2 matrix
237template<typename T>
238[[nodiscard]] constexpr matMxN<2, 2, T> inverse(const matMxN<2, 2, T>& A)
239{
240 T d = T(1) / determinant(A);
241 return d * matMxN<2, 2, T>(vecN<2, T>( A[1][1], -A[0][1]),
242 vecN<2, T>(-A[1][0], A[0][0]));
243}
244
245// inverse of a 3x3 matrix
246template<typename T>
247[[nodiscard]] constexpr matMxN<3, 3, T> inverse(const matMxN<3, 3, T>& A)
248{
249 T d = T(1) / determinant(A);
250 return d * matMxN<3, 3, T>(
251 vecN<3, T>(A[1][1] * A[2][2] - A[1][2] * A[2][1],
252 A[0][2] * A[2][1] - A[0][1] * A[2][2],
253 A[0][1] * A[1][2] - A[0][2] * A[1][1]),
254 vecN<3, T>(A[1][2] * A[2][0] - A[1][0] * A[2][2],
255 A[0][0] * A[2][2] - A[0][2] * A[2][0],
256 A[0][2] * A[1][0] - A[0][0] * A[1][2]),
257 vecN<3, T>(A[1][0] * A[2][1] - A[1][1] * A[2][0],
258 A[0][1] * A[2][0] - A[0][0] * A[2][1],
259 A[0][0] * A[1][1] - A[0][1] * A[1][0]));
260}
261
262// inverse of a 4x4 matrix
263template<typename T>
264[[nodiscard]] constexpr matMxN<4, 4, T> inverse(const matMxN<4, 4, T>& A)
265{
266 // Implementation based on glm: http://glm.g-truc.net
267
268 T c00 = A[2][2] * A[3][3] - A[3][2] * A[2][3];
269 T c02 = A[1][2] * A[3][3] - A[3][2] * A[1][3];
270 T c03 = A[1][2] * A[2][3] - A[2][2] * A[1][3];
271
272 T c04 = A[2][1] * A[3][3] - A[3][1] * A[2][3];
273 T c06 = A[1][1] * A[3][3] - A[3][1] * A[1][3];
274 T c07 = A[1][1] * A[2][3] - A[2][1] * A[1][3];
275
276 T c08 = A[2][1] * A[3][2] - A[3][1] * A[2][2];
277 T c10 = A[1][1] * A[3][2] - A[3][1] * A[1][2];
278 T c11 = A[1][1] * A[2][2] - A[2][1] * A[1][2];
279
280 T c12 = A[2][0] * A[3][3] - A[3][0] * A[2][3];
281 T c14 = A[1][0] * A[3][3] - A[3][0] * A[1][3];
282 T c15 = A[1][0] * A[2][3] - A[2][0] * A[1][3];
283
284 T c16 = A[2][0] * A[3][2] - A[3][0] * A[2][2];
285 T c18 = A[1][0] * A[3][2] - A[3][0] * A[1][2];
286 T c19 = A[1][0] * A[2][2] - A[2][0] * A[1][2];
287
288 T c20 = A[2][0] * A[3][1] - A[3][0] * A[2][1];
289 T c22 = A[1][0] * A[3][1] - A[3][0] * A[1][1];
290 T c23 = A[1][0] * A[2][1] - A[2][0] * A[1][1];
291
292 vecN<4, T> f0(c00, c00, c02, c03);
293 vecN<4, T> f1(c04, c04, c06, c07);
294 vecN<4, T> f2(c08, c08, c10, c11);
295 vecN<4, T> f3(c12, c12, c14, c15);
296 vecN<4, T> f4(c16, c16, c18, c19);
297 vecN<4, T> f5(c20, c20, c22, c23);
298
299 vecN<4, T> v0(A[1][0], A[0][0], A[0][0], A[0][0]);
300 vecN<4, T> v1(A[1][1], A[0][1], A[0][1], A[0][1]);
301 vecN<4, T> v2(A[1][2], A[0][2], A[0][2], A[0][2]);
302 vecN<4, T> v3(A[1][3], A[0][3], A[0][3], A[0][3]);
303
304 vecN<4, T> i0(v1 * f0 - v2 * f1 + v3 * f2);
305 vecN<4, T> i1(v0 * f0 - v2 * f3 + v3 * f4);
306 vecN<4, T> i2(v0 * f1 - v1 * f3 + v3 * f5);
307 vecN<4, T> i3(v0 * f2 - v1 * f4 + v2 * f5);
308
309 vecN<4, T> sa(+1, -1, +1, -1);
310 vecN<4, T> sb(-1, +1, -1, +1);
311 matMxN<4, 4, T> inverse(i0 * sa, i1 * sb, i2 * sa, i3 * sb);
312
313 vecN<4, T> row0(inverse[0][0], inverse[1][0], inverse[2][0], inverse[3][0]);
314 T OneOverDeterminant = static_cast<T>(1) / dot(A[0], row0);
315
316 return inverse * OneOverDeterminant;
317}
318
319// norm-2 squared
320template<int M, int N, typename T>
321[[nodiscard]] constexpr T norm2_2(const matMxN<M, N, T>& A)
322{
324 for (auto i : xrange(N)) t += A[i] * A[i];
325 return sum(t);
326}
327
328// Textual representation. (Only) used to debug unittest.
329template<int M, int N, typename T>
330std::ostream& operator<<(std::ostream& os, const matMxN<M, N, T>& A)
331{
332 for (auto j : xrange(M)) {
333 for (auto i : xrange(N)) {
334 os << A[i][j] << ' ';
335 }
336 os << '\n';
337 }
338 return os;
339}
340
341} // namespace gl
342
343#endif
TclObject t
constexpr matMxN(const vecN< M, T > &x, const vecN< M, T > &y, const vecN< M, T > &z, const vecN< M, T > &w)
Definition: gl_mat.hh:77
constexpr vecN< M, T > & operator[](unsigned i)
Definition: gl_mat.hh:93
friend constexpr bool operator==(const matMxN &x, const matMxN &y)
Definition: gl_mat.hh:109
constexpr matMxN(const vecN<(M< N ? M :N), T > &d)
Definition: gl_mat.hh:45
constexpr matMxN & operator*=(T x)
Definition: gl_mat.hh:103
constexpr const vecN< M, T > & operator[](unsigned i) const
Definition: gl_mat.hh:87
constexpr matMxN()
Definition: gl_mat.hh:35
constexpr matMxN(const matMxN< M2, N2, T > &x)
Definition: gl_mat.hh:55
constexpr matMxN & operator*=(const matMxN< N, N, T > &x)
Definition: gl_mat.hh:104
constexpr matMxN & operator-=(const matMxN &x)
Definition: gl_mat.hh:102
constexpr matMxN & operator+=(const matMxN &x)
Definition: gl_mat.hh:101
constexpr matMxN(const vecN< M, T > &x, const vecN< M, T > &y, const vecN< M, T > &z)
Definition: gl_mat.hh:70
constexpr matMxN(const vecN< M, T > &x, const vecN< M, T > &y)
Definition: gl_mat.hh:63
imat3 i3(ivec3(1, 2, 3), ivec3(4, 5, 6), ivec3(7, 8, 9))
Definition: gl_mat.hh:23
constexpr T norm2_2(const matMxN< M, N, T > &A)
Definition: gl_mat.hh:321
constexpr T dot(const vecN< N, T > &x, const vecN< N, T > &y)
Definition: gl_vec.hh:321
constexpr matMxN< M, N, T > operator*(T x, const matMxN< M, N, T > &A)
Definition: gl_mat.hh:155
constexpr matMxN< N, M, T > transpose(const matMxN< M, N, T > &A)
Definition: gl_mat.hh:191
constexpr T sum(const vecN< N, T > &x)
Definition: gl_vec.hh:307
std::ostream & operator<<(std::ostream &os, const matMxN< M, N, T > &A)
Definition: gl_mat.hh:330
constexpr matMxN< 2, 2, T > inverse(const matMxN< 2, 2, T > &A)
Definition: gl_mat.hh:238
constexpr matMxN< M, N, T > operator-(const matMxN< M, N, T > &A, const matMxN< M, N, T > &B)
Definition: gl_mat.hh:139
constexpr T determinant(const matMxN< 2, 2, T > &A)
Definition: gl_mat.hh:204
constexpr matMxN< M, N, T > operator+(const matMxN< M, N, T > &A, const matMxN< M, N, T > &B)
Definition: gl_mat.hh:130
std::array< const EDStorage, 4 > A
constexpr auto xrange(T e)
Definition: xrange.hh:133