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