1 #ifndef _TCUMATRIX_HPP
2 #define _TCUMATRIX_HPP
3 /*-------------------------------------------------------------------------
4 * drawElements Quality Program Tester Core
5 * ----------------------------------------
6 *
7 * Copyright 2014 The Android Open Source Project
8 *
9 * Licensed under the Apache License, Version 2.0 (the "License");
10 * you may not use this file except in compliance with the License.
11 * You may obtain a copy of the License at
12 *
13 * http://www.apache.org/licenses/LICENSE-2.0
14 *
15 * Unless required by applicable law or agreed to in writing, software
16 * distributed under the License is distributed on an "AS IS" BASIS,
17 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
18 * See the License for the specific language governing permissions and
19 * limitations under the License.
20 *
21 *//*!
22 * \file
23 * \brief Templatized matrix class.
24 *//*--------------------------------------------------------------------*/
25
26 #include "tcuDefs.hpp"
27 #include "tcuVector.hpp"
28 #include "tcuArray.hpp"
29
30 namespace tcu
31 {
32
33 // Templated matrix class.
34 template <typename T, int Rows, int Cols>
35 class Matrix
36 {
37 public:
38 typedef Vector<T, Rows> Element;
39 typedef T Scalar;
40
41 enum
42 {
43 SIZE = Cols,
44 ROWS = Rows,
45 COLS = Cols,
46 };
47
48 Matrix (void);
49 explicit Matrix (const T& src);
50 explicit Matrix (const T src[Rows*Cols]);
51 Matrix (const Vector<T, Rows>& src);
52 Matrix (const Matrix<T, Rows, Cols>& src);
53 ~Matrix (void);
54
55 Matrix<T, Rows, Cols>& operator= (const Matrix<T, Rows, Cols>& src);
56 Matrix<T, Rows, Cols>& operator*= (const Matrix<T, Rows, Cols>& src);
57
58 void setRow (int rowNdx, const Vector<T, Cols>& vec);
59 void setColumn (int colNdx, const Vector<T, Rows>& vec);
60
61 Vector<T, Cols> getRow (int ndx) const;
62 Vector<T, Rows>& getColumn (int ndx);
63 const Vector<T, Rows>& getColumn (int ndx) const;
64
operator [](int ndx)65 Vector<T, Rows>& operator[] (int ndx) { return getColumn(ndx); }
operator [](int ndx) const66 const Vector<T, Rows>& operator[] (int ndx) const { return getColumn(ndx); }
67
operator ()(int row,int col) const68 inline const T& operator() (int row, int col) const { return m_data[col][row]; }
operator ()(int row,int col)69 inline T& operator() (int row, int col) { return m_data[col][row]; }
70
71 Array<T, Rows*Cols> getRowMajorData (void) const;
72 Array<T, Rows*Cols> getColumnMajorData (void) const;
73
74 private:
75 Vector<Vector<T, Rows>, Cols> m_data;
76 } DE_WARN_UNUSED_TYPE;
77
78 // Operators.
79
80 // Mat * Mat.
81 template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
82 Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b);
83
84 // Mat * Vec (column vector).
85 template <typename T, int Rows, int Cols>
86 Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec);
87
88 // Vec * Mat (row vector).
89 template <typename T, int Rows, int Cols>
90 Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx);
91
92 // Further operations
93
94 template <typename T, int Size>
95 struct SquareMatrixOps
96 {
97 static T doDeterminant (const Matrix<T, Size, Size>& mat);
98 static Matrix<T, Size, Size> doInverse (const Matrix<T, Size, Size>& mat);
99 };
100
101 template <typename T>
102 struct SquareMatrixOps<T, 2>
103 {
104 static T doDeterminant (const Matrix<T, 2, 2>& mat);
105 static Matrix<T, 2, 2> doInverse (const Matrix<T, 2, 2>& mat);
106 };
107
108 template <typename T>
109 struct SquareMatrixOps<T, 3>
110 {
111 static T doDeterminant (const Matrix<T, 3, 3>& mat);
112 static Matrix<T, 3, 3> doInverse (const Matrix<T, 3, 3>& mat);
113 };
114
115 template <typename T>
116 struct SquareMatrixOps<T, 4>
117 {
118 static T doDeterminant (const Matrix<T, 4, 4>& mat);
119 static Matrix<T, 4, 4> doInverse (const Matrix<T, 4, 4>& mat);
120 };
121
122 namespace matrix
123 {
124
125 template <typename T, int Size>
determinant(const Matrix<T,Size,Size> & mat)126 T determinant (const Matrix<T, Size, Size>& mat)
127 {
128 return SquareMatrixOps<T, Size>::doDeterminant(mat);
129 }
130
131 template <typename T, int Size>
inverse(const Matrix<T,Size,Size> & mat)132 Matrix<T, Size, Size> inverse (const Matrix<T, Size, Size>& mat)
133 {
134 return SquareMatrixOps<T, Size>::doInverse(mat);
135 }
136
137 } // matrix
138
139 // Template implementations.
140
141 template <typename T>
doDeterminant(const Matrix<T,2,2> & mat)142 T SquareMatrixOps<T, 2>::doDeterminant (const Matrix<T, 2, 2>& mat)
143 {
144 return mat(0,0) * mat(1,1) - mat(1,0) * mat(0,1);
145 }
146
147 template <typename T>
doDeterminant(const Matrix<T,3,3> & mat)148 T SquareMatrixOps<T, 3>::doDeterminant (const Matrix<T, 3, 3>& mat)
149 {
150 return + mat(0,0) * mat(1,1) * mat(2,2)
151 + mat(0,1) * mat(1,2) * mat(2,0)
152 + mat(0,2) * mat(1,0) * mat(2,1)
153 - mat(0,0) * mat(1,2) * mat(2,1)
154 - mat(0,1) * mat(1,0) * mat(2,2)
155 - mat(0,2) * mat(1,1) * mat(2,0);
156 }
157
158 template <typename T>
doDeterminant(const Matrix<T,4,4> & mat)159 T SquareMatrixOps<T, 4>::doDeterminant (const Matrix<T, 4, 4>& mat)
160 {
161 using matrix::determinant;
162
163 const T minorMatrices[4][3*3] =
164 {
165 {
166 mat(1,1), mat(2,1), mat(3,1),
167 mat(1,2), mat(2,2), mat(3,2),
168 mat(1,3), mat(2,3), mat(3,3),
169 },
170 {
171 mat(1,0), mat(2,0), mat(3,0),
172 mat(1,2), mat(2,2), mat(3,2),
173 mat(1,3), mat(2,3), mat(3,3),
174 },
175 {
176 mat(1,0), mat(2,0), mat(3,0),
177 mat(1,1), mat(2,1), mat(3,1),
178 mat(1,3), mat(2,3), mat(3,3),
179 },
180 {
181 mat(1,0), mat(2,0), mat(3,0),
182 mat(1,1), mat(2,1), mat(3,1),
183 mat(1,2), mat(2,2), mat(3,2),
184 }
185 };
186
187 return + mat(0,0) * determinant(Matrix<T, 3, 3>(minorMatrices[0]))
188 - mat(0,1) * determinant(Matrix<T, 3, 3>(minorMatrices[1]))
189 + mat(0,2) * determinant(Matrix<T, 3, 3>(minorMatrices[2]))
190 - mat(0,3) * determinant(Matrix<T, 3, 3>(minorMatrices[3]));
191 }
192
193 template <typename T>
doInverse(const Matrix<T,2,2> & mat)194 Matrix<T, 2, 2> SquareMatrixOps<T, 2>::doInverse (const Matrix<T, 2, 2>& mat)
195 {
196 using matrix::determinant;
197
198 const T det = determinant(mat);
199 Matrix<T, 2, 2> retVal;
200
201 retVal(0, 0) = mat(1, 1) / det;
202 retVal(0, 1) = -mat(0, 1) / det;
203 retVal(1, 0) = -mat(1, 0) / det;
204 retVal(1, 1) = mat(0, 0) / det;
205
206 return retVal;
207 }
208
209 template <typename T>
doInverse(const Matrix<T,3,3> & mat)210 Matrix<T, 3, 3> SquareMatrixOps<T, 3>::doInverse (const Matrix<T, 3, 3>& mat)
211 {
212 // Blockwise inversion
213 using matrix::inverse;
214
215 const T areaA[2*2] =
216 {
217 mat(0,0), mat(0,1),
218 mat(1,0), mat(1,1)
219 };
220 const T areaB[2] =
221 {
222 mat(0,2),
223 mat(1,2),
224 };
225 const T areaC[2] =
226 {
227 mat(2,0), mat(2,1),
228 };
229 const T areaD[1] =
230 {
231 mat(2,2)
232 };
233 const T nullField[4] = { T(0.0f) };
234
235 const Matrix<T, 2, 2> invA = inverse(Matrix<T, 2, 2>(areaA));
236 const Matrix<T, 2, 1> matB = Matrix<T, 2, 1>(areaB);
237 const Matrix<T, 1, 2> matC = Matrix<T, 1, 2>(areaC);
238 const Matrix<T, 1, 1> matD = Matrix<T, 1, 1>(areaD);
239
240 const T schurComplement = T(1.0f) / (matD - matC*invA*matB)(0,0);
241 const Matrix<T, 2, 2> zeroMat = Matrix<T, 2, 2>(nullField);
242
243 const Matrix<T, 2, 2> blockA = invA + invA*matB*schurComplement*matC*invA;
244 const Matrix<T, 2, 1> blockB = (zeroMat-invA)*matB*schurComplement;
245 const Matrix<T, 1, 2> blockC = matC*invA*(-schurComplement);
246 const T blockD = schurComplement;
247
248 const T result[3*3] =
249 {
250 blockA(0,0), blockA(0,1), blockB(0,0),
251 blockA(1,0), blockA(1,1), blockB(1,0),
252 blockC(0,0), blockC(0,1), blockD,
253 };
254
255 return Matrix<T, 3, 3>(result);
256 }
257
258 template <typename T>
doInverse(const Matrix<T,4,4> & mat)259 Matrix<T, 4, 4> SquareMatrixOps<T, 4>::doInverse (const Matrix<T, 4, 4>& mat)
260 {
261 // Blockwise inversion
262 using matrix::inverse;
263
264 const T areaA[2*2] =
265 {
266 mat(0,0), mat(0,1),
267 mat(1,0), mat(1,1)
268 };
269 const T areaB[2*2] =
270 {
271 mat(0,2), mat(0,3),
272 mat(1,2), mat(1,3)
273 };
274 const T areaC[2*2] =
275 {
276 mat(2,0), mat(2,1),
277 mat(3,0), mat(3,1)
278 };
279 const T areaD[2*2] =
280 {
281 mat(2,2), mat(2,3),
282 mat(3,2), mat(3,3)
283 };
284 const T nullField[4] = { T(0.0f) };
285
286 const Matrix<T, 2, 2> invA = inverse(Matrix<T, 2, 2>(areaA));
287 const Matrix<T, 2, 2> matB = Matrix<T, 2, 2>(areaB);
288 const Matrix<T, 2, 2> matC = Matrix<T, 2, 2>(areaC);
289 const Matrix<T, 2, 2> matD = Matrix<T, 2, 2>(areaD);
290
291 const Matrix<T, 2, 2> schurComplement = inverse(matD - matC*invA*matB);
292 const Matrix<T, 2, 2> zeroMat = Matrix<T, 2, 2>(nullField);
293
294 const Matrix<T, 2, 2> blockA = invA + invA*matB*schurComplement*matC*invA;
295 const Matrix<T, 2, 2> blockB = (zeroMat-invA)*matB*schurComplement;
296 const Matrix<T, 2, 2> blockC = (zeroMat-schurComplement)*matC*invA;
297 const Matrix<T, 2, 2> blockD = schurComplement;
298
299 const T result[4*4] =
300 {
301 blockA(0,0), blockA(0,1), blockB(0,0), blockB(0,1),
302 blockA(1,0), blockA(1,1), blockB(1,0), blockB(1,1),
303 blockC(0,0), blockC(0,1), blockD(0,0), blockD(0,1),
304 blockC(1,0), blockC(1,1), blockD(1,0), blockD(1,1),
305 };
306
307 return Matrix<T, 4, 4>(result);
308 }
309
310 // Initialize to identity.
311 template <typename T, int Rows, int Cols>
Matrix(void)312 Matrix<T, Rows, Cols>::Matrix (void)
313 {
314 for (int row = 0; row < Rows; row++)
315 for (int col = 0; col < Cols; col++)
316 (*this)(row, col) = (row == col) ? T(1) : T(0);
317 }
318
319 // Initialize to diagonal matrix.
320 template <typename T, int Rows, int Cols>
Matrix(const T & src)321 Matrix<T, Rows, Cols>::Matrix (const T& src)
322 {
323 for (int row = 0; row < Rows; row++)
324 for (int col = 0; col < Cols; col++)
325 (*this)(row, col) = (row == col) ? src : T(0);
326 }
327
328 // Initialize from data array.
329 template <typename T, int Rows, int Cols>
Matrix(const T src[Rows * Cols])330 Matrix<T, Rows, Cols>::Matrix (const T src[Rows*Cols])
331 {
332 for (int row = 0; row < Rows; row++)
333 for (int col = 0; col < Cols; col++)
334 (*this)(row, col) = src[row*Cols + col];
335 }
336
337 // Initialize to diagonal matrix.
338 template <typename T, int Rows, int Cols>
Matrix(const Vector<T,Rows> & src)339 Matrix<T, Rows, Cols>::Matrix (const Vector<T, Rows>& src)
340 {
341 DE_STATIC_ASSERT(Rows == Cols);
342 for (int row = 0; row < Rows; row++)
343 for (int col = 0; col < Cols; col++)
344 (*this)(row, col) = (row == col) ? src.m_data[row] : T(0);
345 }
346
347 // Copy constructor.
348 template <typename T, int Rows, int Cols>
Matrix(const Matrix<T,Rows,Cols> & src)349 Matrix<T, Rows, Cols>::Matrix (const Matrix<T, Rows, Cols>& src)
350 {
351 *this = src;
352 }
353
354 // Destructor.
355 template <typename T, int Rows, int Cols>
~Matrix(void)356 Matrix<T, Rows, Cols>::~Matrix (void)
357 {
358 }
359
360 // Assignment operator.
361 template <typename T, int Rows, int Cols>
operator =(const Matrix<T,Rows,Cols> & src)362 Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator= (const Matrix<T, Rows, Cols>& src)
363 {
364 for (int row = 0; row < Rows; row++)
365 for (int col = 0; col < Cols; col++)
366 (*this)(row, col) = src(row, col);
367 return *this;
368 }
369
370 // Multipy and assign op
371 template <typename T, int Rows, int Cols>
operator *=(const Matrix<T,Rows,Cols> & src)372 Matrix<T, Rows, Cols>& Matrix<T, Rows, Cols>::operator*= (const Matrix<T, Rows, Cols>& src)
373 {
374 *this = *this * src;
375 return *this;
376 }
377
378 template <typename T, int Rows, int Cols>
setRow(int rowNdx,const Vector<T,Cols> & vec)379 void Matrix<T, Rows, Cols>::setRow (int rowNdx, const Vector<T, Cols>& vec)
380 {
381 for (int col = 0; col < Cols; col++)
382 (*this)(rowNdx, col) = vec.m_data[col];
383 }
384
385 template <typename T, int Rows, int Cols>
setColumn(int colNdx,const Vector<T,Rows> & vec)386 void Matrix<T, Rows, Cols>::setColumn (int colNdx, const Vector<T, Rows>& vec)
387 {
388 m_data[colNdx] = vec;
389 }
390
391 template <typename T, int Rows, int Cols>
getRow(int rowNdx) const392 Vector<T, Cols> Matrix<T, Rows, Cols>::getRow (int rowNdx) const
393 {
394 Vector<T, Cols> res;
395 for (int col = 0; col < Cols; col++)
396 res[col] = (*this)(rowNdx, col);
397 return res;
398 }
399
400 template <typename T, int Rows, int Cols>
getColumn(int colNdx)401 Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx)
402 {
403 return m_data[colNdx];
404 }
405
406 template <typename T, int Rows, int Cols>
getColumn(int colNdx) const407 const Vector<T, Rows>& Matrix<T, Rows, Cols>::getColumn (int colNdx) const
408 {
409 return m_data[colNdx];
410 }
411
412 template <typename T, int Rows, int Cols>
getColumnMajorData(void) const413 Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getColumnMajorData (void) const
414 {
415 Array<T, Rows*Cols> a;
416 T* dst = a.getPtr();
417 for (int col = 0; col < Cols; col++)
418 for (int row = 0; row < Rows; row++)
419 *dst++ = (*this)(row, col);
420 return a;
421 }
422
423 template <typename T, int Rows, int Cols>
getRowMajorData(void) const424 Array<T, Rows*Cols> Matrix<T, Rows, Cols>::getRowMajorData (void) const
425 {
426 Array<T, Rows*Cols> a;
427 T* dst = a.getPtr();
428 for (int row = 0; row < Rows; row++)
429 for (int col = 0; col < Cols; col++)
430 *dst++ = (*this)(row, col);
431 return a;
432 }
433
434 // Multiplication of two matrices.
435 template <typename T, int Rows0, int Cols0, int Rows1, int Cols1>
operator *(const Matrix<T,Rows0,Cols0> & a,const Matrix<T,Rows1,Cols1> & b)436 Matrix<T, Rows0, Cols1> operator* (const Matrix<T, Rows0, Cols0>& a, const Matrix<T, Rows1, Cols1>& b)
437 {
438 DE_STATIC_ASSERT(Cols0 == Rows1);
439 Matrix<T, Rows0, Cols1> res;
440 for (int row = 0; row < Rows0; row++)
441 {
442 for (int col = 0; col < Cols1; col++)
443 {
444 T v = T(0);
445 for (int ndx = 0; ndx < Cols0; ndx++)
446 v += a(row,ndx) * b(ndx,col);
447 res(row,col) = v;
448 }
449 }
450 return res;
451 }
452
453 // Multiply of matrix with column vector.
454 template <typename T, int Rows, int Cols>
operator *(const Matrix<T,Rows,Cols> & mtx,const Vector<T,Cols> & vec)455 Vector<T, Rows> operator* (const Matrix<T, Rows, Cols>& mtx, const Vector<T, Cols>& vec)
456 {
457 Vector<T, Rows> res;
458 for (int row = 0; row < Rows; row++)
459 {
460 T v = T(0);
461 for (int col = 0; col < Cols; col++)
462 v += mtx(row,col) * vec.m_data[col];
463 res.m_data[row] = v;
464 }
465 return res;
466 }
467
468 // Multiply of matrix with row vector.
469 template <typename T, int Rows, int Cols>
operator *(const Vector<T,Rows> & vec,const Matrix<T,Rows,Cols> & mtx)470 Vector<T, Cols> operator* (const Vector<T, Rows>& vec, const Matrix<T, Rows, Cols>& mtx)
471 {
472 Vector<T, Cols> res;
473 for (int col = 0; col < Cols; col++)
474 {
475 T v = T(0);
476 for (int row = 0; row < Rows; row++)
477 v += mtx(row,col) * vec.m_data[row];
478 res.m_data[col] = v;
479 }
480 return res;
481 }
482
483 // Common typedefs.
484 typedef Matrix<float, 2, 2> Matrix2f;
485 typedef Matrix<float, 3, 3> Matrix3f;
486 typedef Matrix<float, 4, 4> Matrix4f;
487 typedef Matrix<double, 2, 2> Matrix2d;
488 typedef Matrix<double, 3, 3> Matrix3d;
489 typedef Matrix<double, 4, 4> Matrix4d;
490
491 // GLSL-style naming \note CxR.
492 typedef Matrix2f Mat2;
493 typedef Matrix<float, 3, 2> Mat2x3;
494 typedef Matrix<float, 4, 2> Mat2x4;
495 typedef Matrix<float, 2, 3> Mat3x2;
496 typedef Matrix3f Mat3;
497 typedef Matrix<float, 4, 3> Mat3x4;
498 typedef Matrix<float, 2, 4> Mat4x2;
499 typedef Matrix<float, 3, 4> Mat4x3;
500 typedef Matrix4f Mat4;
501
502 // Matrix-scalar operators.
503
504 template <typename T, int Rows, int Cols>
operator +(const Matrix<T,Rows,Cols> & mtx,T scalar)505 Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& mtx, T scalar)
506 {
507 Matrix<T, Rows, Cols> res;
508 for (int col = 0; col < Cols; col++)
509 for (int row = 0; row < Rows; row++)
510 res(row, col) = mtx(row, col) + scalar;
511 return res;
512 }
513
514 template <typename T, int Rows, int Cols>
operator -(const Matrix<T,Rows,Cols> & mtx,T scalar)515 Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& mtx, T scalar)
516 {
517 Matrix<T, Rows, Cols> res;
518 for (int col = 0; col < Cols; col++)
519 for (int row = 0; row < Rows; row++)
520 res(row, col) = mtx(row, col) - scalar;
521 return res;
522 }
523
524 template <typename T, int Rows, int Cols>
operator *(const Matrix<T,Rows,Cols> & mtx,T scalar)525 Matrix<T, Rows, Cols> operator* (const Matrix<T, Rows, Cols>& mtx, T scalar)
526 {
527 Matrix<T, Rows, Cols> res;
528 for (int col = 0; col < Cols; col++)
529 for (int row = 0; row < Rows; row++)
530 res(row, col) = mtx(row, col) * scalar;
531 return res;
532 }
533
534 template <typename T, int Rows, int Cols>
operator /(const Matrix<T,Rows,Cols> & mtx,T scalar)535 Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& mtx, T scalar)
536 {
537 Matrix<T, Rows, Cols> res;
538 for (int col = 0; col < Cols; col++)
539 for (int row = 0; row < Rows; row++)
540 res(row, col) = mtx(row, col) / scalar;
541 return res;
542 }
543
544 // Matrix-matrix component-wise operators.
545
546 template <typename T, int Rows, int Cols>
operator +(const Matrix<T,Rows,Cols> & a,const Matrix<T,Rows,Cols> & b)547 Matrix<T, Rows, Cols> operator+ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
548 {
549 Matrix<T, Rows, Cols> res;
550 for (int col = 0; col < Cols; col++)
551 for (int row = 0; row < Rows; row++)
552 res(row, col) = a(row, col) + b(row, col);
553 return res;
554 }
555
556 template <typename T, int Rows, int Cols>
operator -(const Matrix<T,Rows,Cols> & a,const Matrix<T,Rows,Cols> & b)557 Matrix<T, Rows, Cols> operator- (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
558 {
559 Matrix<T, Rows, Cols> res;
560 for (int col = 0; col < Cols; col++)
561 for (int row = 0; row < Rows; row++)
562 res(row, col) = a(row, col) - b(row, col);
563 return res;
564 }
565
566 template <typename T, int Rows, int Cols>
operator /(const Matrix<T,Rows,Cols> & a,const Matrix<T,Rows,Cols> & b)567 Matrix<T, Rows, Cols> operator/ (const Matrix<T, Rows, Cols>& a, const Matrix<T, Rows, Cols>& b)
568 {
569 Matrix<T, Rows, Cols> res;
570 for (int col = 0; col < Cols; col++)
571 for (int row = 0; row < Rows; row++)
572 res(row, col) = a(row, col) / b(row, col);
573 return res;
574 }
575
576 } // tcu
577
578 #endif // _TCUMATRIX_HPP
579