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 
23 namespace 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).
29 template<int M, int N, typename T> class matMxN
30 {
31 public:
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 private:
107  vecN<M, T> c[N];
108 };
109 
110 
111 // Convenience typedefs (same names as used by GLSL).
115 
116 
117 // -- Matrix functions --
118 
119 // matrix equality / inequality
120 template<int M, int N, typename T>
121 [[nodiscard]] constexpr bool operator==(const matMxN<M, N, T>& A, const matMxN<M, N, T>& B)
122 {
123  for (auto i : xrange(N)) if (A[i] != B[i]) return false;
124  return true;
125 }
126 template<int M, int N, typename T>
127 [[nodiscard]] constexpr bool operator!=(const matMxN<M, N, T>& A, const matMxN<M, N, T>& B)
128 {
129  return !(A == B);
130 }
131 
132 // matrix + matrix
133 template<int M, int N, typename T>
134 [[nodiscard]] constexpr matMxN<M, N, T> operator+(const matMxN<M, N, T>& A, const matMxN<M, N, T>& B)
135 {
137  for (auto i : xrange(N)) R[i] = A[i] + B[i];
138  return R;
139 }
140 
141 // matrix - matrix
142 template<int M, int N, typename T>
143 [[nodiscard]] constexpr matMxN<M, N, T> operator-(const matMxN<M, N, T>& A, const matMxN<M, N, T>& B)
144 {
146  for (auto i : xrange(N)) R[i] = A[i] - B[i];
147  return R;
148 }
149 
150 // matrix negation
151 template<int M, int N, typename T>
152 [[nodiscard]] constexpr matMxN<M, N, T> operator-(const matMxN<M, N, T>& A)
153 {
154  return matMxN<M, N, T>(vecN<(M < N ? M : N), T>()) - A;
155 }
156 
157 // scalar * matrix
158 template<int M, int N, typename T>
159 [[nodiscard]] constexpr matMxN<M, N, T> operator*(T x, const matMxN<M, N, T>& A)
160 {
162  for (auto i : xrange(N)) R[i] = x * A[i];
163  return R;
164 }
165 
166 // matrix * scalar
167 template<int M, int N, typename T>
168 [[nodiscard]] constexpr matMxN<M, N, T> operator*(const matMxN<M, N, T>& A, T x)
169 {
171  for (auto i : xrange(N)) R[i] = A[i] * x;
172  return R;
173 }
174 
175 // matrix * column-vector
176 template<int M, int N, typename T>
177 [[nodiscard]] constexpr vecN<M, T> operator*(const matMxN<M, N, T>& A, const vecN<N, T>& x)
178 {
179  vecN<M, T> r;
180  for (auto i : xrange(N)) r += A[i] * x[i];
181  return r;
182 }
183 
184 // matrix * matrix
185 template<int M, int N, int O, typename T>
186 [[nodiscard]] constexpr matMxN<M, O, T> operator*(const matMxN<M, N, T>& A, const matMxN<N, O, T>& B)
187 {
189  for (auto i : xrange(O)) R[i] = A * B[i];
190  return R;
191 }
192 
193 // matrix transpose
194 template<int M, int N, typename T>
195 [[nodiscard]] constexpr matMxN<N, M, T> transpose(const matMxN<M, N, T>& A)
196 {
198  for (auto i : xrange(N)) {
199  for (auto j : xrange(M)) {
200  R[j][i] = A[i][j];
201  }
202  }
203  return R;
204 }
205 
206 // determinant of a 2x2 matrix
207 template<typename T>
208 [[nodiscard]] constexpr T determinant(const matMxN<2, 2, T>& A)
209 {
210  return A[0][0] * A[1][1] - A[0][1] * A[1][0];
211 }
212 
213 // determinant of a 3x3 matrix
214 template<typename T>
215 [[nodiscard]] constexpr T determinant(const matMxN<3, 3, T>& A)
216 {
217  return A[0][0] * (A[1][1] * A[2][2] - A[1][2] * A[2][1])
218  - A[1][0] * (A[0][1] * A[2][2] - A[0][2] * A[2][1])
219  + A[2][0] * (A[0][1] * A[1][2] - A[0][2] * A[1][1]);
220 }
221 
222 // determinant of a 4x4 matrix
223 template<typename T>
224 [[nodiscard]] constexpr T determinant(const matMxN<4, 4, T>& A)
225 {
226  // Implementation based on glm: http://glm.g-truc.net
227  T f0 = A[2][2] * A[3][3] - A[3][2] * A[2][3];
228  T f1 = A[2][1] * A[3][3] - A[3][1] * A[2][3];
229  T f2 = A[2][1] * A[3][2] - A[3][1] * A[2][2];
230  T f3 = A[2][0] * A[3][3] - A[3][0] * A[2][3];
231  T f4 = A[2][0] * A[3][2] - A[3][0] * A[2][2];
232  T f5 = A[2][0] * A[3][1] - A[3][0] * A[2][1];
233  vecN<4, T> c((A[1][1] * f0 - A[1][2] * f1 + A[1][3] * f2),
234  (A[1][2] * f3 - A[1][3] * f4 - A[1][0] * f0),
235  (A[1][0] * f1 - A[1][1] * f3 + A[1][3] * f5),
236  (A[1][1] * f4 - A[1][2] * f5 - A[1][0] * f2));
237  return dot(A[0], c);
238 }
239 
240 // inverse of a 2x2 matrix
241 template<typename T>
242 [[nodiscard]] constexpr matMxN<2, 2, T> inverse(const matMxN<2, 2, T>& A)
243 {
244  T d = T(1) / determinant(A);
245  return d * matMxN<2, 2, T>(vecN<2, T>( A[1][1], -A[0][1]),
246  vecN<2, T>(-A[1][0], A[0][0]));
247 }
248 
249 // inverse of a 3x3 matrix
250 template<typename T>
251 [[nodiscard]] constexpr matMxN<3, 3, T> inverse(const matMxN<3, 3, T>& A)
252 {
253  T d = T(1) / determinant(A);
254  return d * matMxN<3, 3, T>(
255  vecN<3, T>(A[1][1] * A[2][2] - A[1][2] * A[2][1],
256  A[0][2] * A[2][1] - A[0][1] * A[2][2],
257  A[0][1] * A[1][2] - A[0][2] * A[1][1]),
258  vecN<3, T>(A[1][2] * A[2][0] - A[1][0] * A[2][2],
259  A[0][0] * A[2][2] - A[0][2] * A[2][0],
260  A[0][2] * A[1][0] - A[0][0] * A[1][2]),
261  vecN<3, T>(A[1][0] * A[2][1] - A[1][1] * A[2][0],
262  A[0][1] * A[2][0] - A[0][0] * A[2][1],
263  A[0][0] * A[1][1] - A[0][1] * A[1][0]));
264 }
265 
266 // inverse of a 4x4 matrix
267 template<typename T>
268 [[nodiscard]] constexpr matMxN<4, 4, T> inverse(const matMxN<4, 4, T>& A)
269 {
270  // Implementation based on glm: http://glm.g-truc.net
271 
272  T c00 = A[2][2] * A[3][3] - A[3][2] * A[2][3];
273  T c02 = A[1][2] * A[3][3] - A[3][2] * A[1][3];
274  T c03 = A[1][2] * A[2][3] - A[2][2] * A[1][3];
275 
276  T c04 = A[2][1] * A[3][3] - A[3][1] * A[2][3];
277  T c06 = A[1][1] * A[3][3] - A[3][1] * A[1][3];
278  T c07 = A[1][1] * A[2][3] - A[2][1] * A[1][3];
279 
280  T c08 = A[2][1] * A[3][2] - A[3][1] * A[2][2];
281  T c10 = A[1][1] * A[3][2] - A[3][1] * A[1][2];
282  T c11 = A[1][1] * A[2][2] - A[2][1] * A[1][2];
283 
284  T c12 = A[2][0] * A[3][3] - A[3][0] * A[2][3];
285  T c14 = A[1][0] * A[3][3] - A[3][0] * A[1][3];
286  T c15 = A[1][0] * A[2][3] - A[2][0] * A[1][3];
287 
288  T c16 = A[2][0] * A[3][2] - A[3][0] * A[2][2];
289  T c18 = A[1][0] * A[3][2] - A[3][0] * A[1][2];
290  T c19 = A[1][0] * A[2][2] - A[2][0] * A[1][2];
291 
292  T c20 = A[2][0] * A[3][1] - A[3][0] * A[2][1];
293  T c22 = A[1][0] * A[3][1] - A[3][0] * A[1][1];
294  T c23 = A[1][0] * A[2][1] - A[2][0] * A[1][1];
295 
296  vecN<4, T> f0(c00, c00, c02, c03);
297  vecN<4, T> f1(c04, c04, c06, c07);
298  vecN<4, T> f2(c08, c08, c10, c11);
299  vecN<4, T> f3(c12, c12, c14, c15);
300  vecN<4, T> f4(c16, c16, c18, c19);
301  vecN<4, T> f5(c20, c20, c22, c23);
302 
303  vecN<4, T> v0(A[1][0], A[0][0], A[0][0], A[0][0]);
304  vecN<4, T> v1(A[1][1], A[0][1], A[0][1], A[0][1]);
305  vecN<4, T> v2(A[1][2], A[0][2], A[0][2], A[0][2]);
306  vecN<4, T> v3(A[1][3], A[0][3], A[0][3], A[0][3]);
307 
308  vecN<4, T> i0(v1 * f0 - v2 * f1 + v3 * f2);
309  vecN<4, T> i1(v0 * f0 - v2 * f3 + v3 * f4);
310  vecN<4, T> i2(v0 * f1 - v1 * f3 + v3 * f5);
311  vecN<4, T> i3(v0 * f2 - v1 * f4 + v2 * f5);
312 
313  vecN<4, T> sa(+1, -1, +1, -1);
314  vecN<4, T> sb(-1, +1, -1, +1);
315  matMxN<4, 4, T> inverse(i0 * sa, i1 * sb, i2 * sa, i3 * sb);
316 
317  vecN<4, T> row0(inverse[0][0], inverse[1][0], inverse[2][0], inverse[3][0]);
318  T OneOverDeterminant = static_cast<T>(1) / dot(A[0], row0);
319 
320  return inverse * OneOverDeterminant;
321 }
322 
323 // norm-2 squared
324 template<int M, int N, typename T>
325 [[nodiscard]] constexpr T norm2_2(const matMxN<M, N, T>& A)
326 {
327  vecN<M, T> t;
328  for (auto i : xrange(N)) t += A[i] * A[i];
329  return sum(t);
330 }
331 
332 // Textual representation. (Only) used to debug unittest.
333 template<int M, int N, typename T>
334 std::ostream& operator<<(std::ostream& os, const matMxN<M, N, T>& A)
335 {
336  for (auto j : xrange(M)) {
337  for (auto i : xrange(N)) {
338  os << A[i][j] << ' ';
339  }
340  os << '\n';
341  }
342  return os;
343 }
344 
345 } // namespace gl
346 
347 #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 const vecN< M, T > & operator[](unsigned i) const
Definition: gl_mat.hh:87
constexpr matMxN & operator-=(const matMxN &x)
Definition: gl_mat.hh:102
constexpr matMxN(const vecN<(M< N ? M :N), T > &d)
Definition: gl_mat.hh:45
constexpr matMxN & operator*=(const matMxN< N, N, T > &x)
Definition: gl_mat.hh:104
constexpr matMxN & operator+=(const matMxN &x)
Definition: gl_mat.hh:101
constexpr matMxN()
Definition: gl_mat.hh:35
constexpr matMxN(const matMxN< M2, N2, T > &x)
Definition: gl_mat.hh:55
constexpr vecN< M, T > & operator[](unsigned i)
Definition: gl_mat.hh:93
constexpr matMxN & operator*=(T x)
Definition: gl_mat.hh:103
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:325
constexpr matMxN< 2, 2, T > inverse(const matMxN< 2, 2, T > &A)
Definition: gl_mat.hh:242
constexpr T dot(const vecN< N, T > &x, const vecN< N, T > &y)
Definition: gl_vec.hh:324
std::ostream & operator<<(std::ostream &os, const matMxN< M, N, T > &A)
Definition: gl_mat.hh:334
constexpr matMxN< M, N, T > operator-(const matMxN< M, N, T > &A, const matMxN< M, N, T > &B)
Definition: gl_mat.hh:143
constexpr bool operator!=(const matMxN< M, N, T > &A, const matMxN< M, N, T > &B)
Definition: gl_mat.hh:127
constexpr matMxN< N, M, T > transpose(const matMxN< M, N, T > &A)
Definition: gl_mat.hh:195
constexpr T sum(const vecN< N, T > &x)
Definition: gl_vec.hh:310
constexpr T determinant(const matMxN< 2, 2, T > &A)
Definition: gl_mat.hh:208
constexpr bool operator==(const matMxN< M, N, T > &A, const matMxN< M, N, T > &B)
Definition: gl_mat.hh:121
constexpr matMxN< M, N, T > operator*(T x, const matMxN< M, N, T > &A)
Definition: gl_mat.hh:159
constexpr matMxN< M, N, T > operator+(const matMxN< M, N, T > &A, const matMxN< M, N, T > &B)
Definition: gl_mat.hh:134
constexpr unsigned N2
Definition: ResampleHQ.cc:230
constexpr unsigned N
Definition: ResampleHQ.cc:228
constexpr auto R
Definition: MSXMixer.cc:298
constexpr KeyMatrixPosition x
Keyboard bindings.
Definition: Keyboard.cc:118
constexpr auto xrange(T e)
Definition: xrange.hh:155