framework/common/tcuMatrix.hpp

   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
  65         Vector<T, Rows>&                                operator[]                      (int ndx)                                       { return getColumn(ndx);        }
  66         const Vector<T, Rows>&                  operator[]                      (int ndx) const                         { return getColumn(ndx);        }
  67
  68         inline const T&                                 operator()                      (int row, int col) const        { return m_data[col][row];      }
  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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 356 Matrix<T, Rows, Cols>::~Matrix (void)
 357 {
 358 }
 359
 360 // Assignment operator.
 361 template <typename T, int Rows, int Cols>
 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>
 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>
 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>
 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>
 392 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>
 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>
 407 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>
 413 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>
 424 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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>
 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