maths/Matrix.h

   1 /***************************************************************************
   2  *   This file is part of Tecorrec.                                        *
   3  *   Copyright 2008 James Hogan <james@albanarts.com>                      *
   4  *                                                                         *
   5  *   Tecorrec is free software: you can redistribute it and/or modify      *
   6  *   it under the terms of the GNU General Public License as published by  *
   7  *   the Free Software Foundation, either version 2 of the License, or     *
   8  *   (at your option) any later version.                                   *
   9  *                                                                         *
  10  *   Tecorrec is distributed in the hope that it will be useful,           *
  11  *   but WITHOUT ANY WARRANTY; without even the implied warranty of        *
  12  *   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the         *
  13  *   GNU General Public License for more details.                          *
  14  *                                                                         *
  15  *   You should have received a copy of the GNU General Public License     *
  16  *   along with Tecorrec.  If not, write to the Free Software Foundation,  *
  17  *   Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.   *
  18  ***************************************************************************/
  19
  20 /**
  21  * @file Matrix.h
  22  * @brief Square matrix class.
  23  * @author James Hogan (james@albanarts.com)
  24  * @note Copyright (C) 2007
  25  *
  26  * @version 1 : From Computer Graphics and Visualisation Open Assessment.
  27  */
  28
  29 #ifndef _maths_Matrix_h
  30 #define _maths_Matrix_h
  31
  32 #include "Vector.h"
  33
  34 // For memcpy
  35 #include <cstring>
  36
  37 namespace maths
  38 {
  39   // This isn't always used with matricies. include Quaternion.h if you wish to use it.
  40   template <typename T>
  41       class Quaternion;
  42
  43   template <int N, typename T>
  44       class Matrix;
  45
  46   /// Get an arbitrary matrix where the z axis points in the direction of Z.
  47   Matrix<3,float> getTbnMatrix(Vector<3, float> vZ);
  48
  49   Matrix<4,float> RotateRadMatrix44(char axis, float rad);
  50   Matrix<4,float> RotateRadMatrix44(const maths::Vector<3,float> &axis, float rad);
  51   Matrix<4,float> TranslateMatrix44(const maths::Vector<3,float> & v);
  52   Matrix<4,float> ScaleMatrix44(const maths::Vector<4,float> & v);
  53   inline Matrix<4,float> ScaleMatrix44(const maths::Vector<3,float> & v);
  54
  55   /// Square tranformation matrix.
  56   /**
  57    * The matrix is stored in a column major fashion so as to be directly.
  58    *  compatible with OpenGL
  59    * @param N int Size of width and height of the matrix.
  60    * @param T typename Type of the components of the matrix.
  61    */
  62   template <int N, typename T>
  63       class Matrix
  64   {
  65     public:
  66       /// Columns
  67       Vector<N, T> col[N];
  68
  69       /// Default constructor (does not initialise values).
  70       inline Matrix()
  71       {
  72         // let the vector constructor deal with the initialisation
  73       }
  74
  75       /// Diagonal matrix constructor
  76       inline explicit Matrix(T diagonal)
  77       {
  78         for (int i = 0; i < N; ++i) {
  79           for (int j = 0; j < i; ++j) {
  80             col[i][j] = col[j][i] = maths::zero<T>();
  81           }
  82           col[i][i] = diagonal;
  83         }
  84       }
  85
  86       /// Diagonal matrix constructor
  87       inline explicit Matrix(Vector<N, T> diagonal)
  88       {
  89         for (int i = 0; i < N; ++i) {
  90           for (int j = 0; j < i; ++j) {
  91             col[i][j] = col[j][i] = maths::zero<T>();
  92           }
  93           col[i][i] = diagonal[i];
  94         }
  95       }
  96
  97       inline Matrix(Vector<N, T> x, Vector<N, T> y)
  98       {
  99         col[0] = x;
 100         col[1] = y;
 101         if (N > 2) {
 102           for (int i = 2; i < N; ++i) {
 103             col[i] = (Vector<N, T>)maths::zero<T>();
 104           }
 105         }
 106       }
 107       inline Matrix(Vector<N, T> x, Vector<N, T> y, Vector<N, T> z)
 108       {
 109         col[0] = x;
 110         col[1] = y;
 111         if (N > 2) {
 112           col[2] = z;
 113           if (N > 3) {
 114             for (int i = 3; i < N; ++i) {
 115               col[i] = (Vector<N, T>)maths::zero<T>();
 116             }
 117           }
 118         }
 119       }
 120       inline Matrix(Vector<N, T> x, Vector<N, T> y, Vector<N, T> z, Vector<N, T> w)
 121       {
 122         col[0] = x;
 123         col[1] = y;
 124         if (N > 2) {
 125           col[2] = z;
 126           if (N > 3) {
 127             col[3] = w;
 128             if (N > 4) {
 129               for (int i = 4; i < N; ++i) {
 130                 col[i] = (Vector<N, T>)maths::zero<T>();
 131               }
 132             }
 133           }
 134         }
 135       }
 136
 137       /// Copy the matrix directly from an array.
 138       /**
 139        * @param vec Pointer to array of values in column major format.
 140        */
 141       inline explicit Matrix(const T * vec)
 142       {
 143         std::memcpy(&col[0][0], vec, sizeof(T)*(N*N));
 144       }
 145
 146       /// Copy the matrix from an array of vectors.
 147       /**
 148        * @param columns Pointer to array of column vectors.
 149        */
 150       inline explicit Matrix(const Vector<N, T> * columns)
 151       {
 152         std::memcpy(col, columns.col, sizeof(T)*(N*N));
 153       }
 154       /// Convert quaternion to matrix
 155       explicit Matrix(const Quaternion<T> &);
 156
 157       /// From something else
 158       /**
 159        * Basically we want to prevent accidental double conversions between types.
 160        */
 161       template <typename U>
 162           explicit Matrix(const U &);
 163
 164       /// Convert, downsize or upsize another matrix.
 165       /**
 166        * If the matrix is being downsized, the extreme rows and columns are
 167        *  ignored.
 168        * If the matrix is being upsized, the extreme rows and columns are filled
 169        *  with the identity matrix.
 170        * @param copy The matrix object to construct from.
 171        */
 172       template <int M, typename U> inline explicit Matrix(const Matrix<M,U> & copy)
 173       {
 174         const int copyRange = (N < M ? N : M);
 175     // Copy as many columns as possible from copy
 176         int i = 0;
 177         for (; i < copyRange; ++i) {
 178       // Copy as many values as possible from the column
 179           int j = 0;
 180           for (; j < copyRange; ++j)
 181             col[i][j] = (T)copy.col[i][j];
 182           if (M < N) {
 183         // Generate any remaining values from the identity matrix.
 184             for (; j < N; ++j)
 185               col[i][j] = (i==j ? (T)1 : maths::zero<T>());
 186           }
 187         }
 188         if (M < N) {
 189       // Generate any remaining columns from the identity matrix.
 190           for (; i < N; ++i) {
 191             for (int j = 0; j < N; ++j)
 192               col[i][j] = (i==j ? (T)1 : maths::zero<T>());
 193           }
 194         }
 195       }
 196
 197       /// Index a column in the matrix.
 198       /**
 199        * @param columnIndex The index of the column starting at 0.
 200        * @return The column as a vector.
 201        */
 202       inline Vector<N, T> & operator [] (int columnIndex)
 203       {
 204         return col[columnIndex];
 205       }
 206       /// Index a column in the matrix.
 207       /**
 208        * @param columnIndex The index of the column starting at 0.
 209        * @return The column as a vector.
 210        */
 211       inline const Vector<N, T> & operator [] (int columnIndex) const
 212       {
 213         return col[columnIndex];
 214       }
 215
 216       /// Get a pointer to the values of the matrix.
 217       /**
 218        * @return Pointer to the values of the array in column major format.
 219        */
 220       inline operator T * ()
 221       {
 222         return (T*)col;
 223       }
 224       /// Get a pointer to the values of the matrix.
 225       /**
 226        * @return Pointer to the values of the array in column major format.
 227        */
 228       inline operator const T * () const
 229       {
 230         return (const T*)col;
 231       }
 232
 233       inline Matrix & operator += (const Matrix & rhs)
 234       {
 235         for (int i = 0; i < N; ++i)
 236           col[i] += rhs.col[i];
 237         return *this;
 238       }
 239       inline Matrix operator + (const Matrix & rhs) const
 240       {
 241         Matrix ret = *this;
 242         return ret += rhs;
 243       }
 244
 245       inline Matrix & operator -= (const Matrix & rhs)
 246       {
 247         for (int i = 0; i < N; ++i)
 248           col[i] -= rhs.col[i];
 249         return *this;
 250       }
 251       inline Matrix operator - (const Matrix & rhs) const
 252       {
 253         Matrix ret = *this;
 254         return ret -= rhs;
 255       }
 256
 257       inline Matrix & operator *= (const T rhs)
 258       {
 259         for (int i = 0; i < N; ++i)
 260           col[i] *= rhs;
 261         return *this;
 262       }
 263       inline Matrix operator * (const T rhs) const
 264       {
 265         Matrix ret = *this;
 266         return ret *= rhs;
 267       }
 268
 269       inline friend Matrix operator * (const T lhs, const Matrix & rhs)
 270       {
 271         return rhs*lhs;
 272       }
 273
 274       template <typename U>
 275       inline friend Vector<N,U> operator * (const Matrix & lhs, const Vector<N,U> & rhs)
 276       {
 277         Vector<N,U> ret;
 278         for (int row = 0; row < N; ++row) {
 279           ret[row] = maths::zero<U>();
 280           for (int i = 0; i < N; ++i) {
 281             ret[row] += lhs[i][row] * rhs[i];
 282           }
 283         }
 284         return ret;
 285       }
 286       template <typename U>
 287       inline friend Vector<N,U> operator * (const Vector<N,U> & lhs, const Matrix & rhs)
 288       {
 289         Vector<N,U> ret;
 290         for (int col = 0; col < N; ++col) {
 291           ret[col] = maths::zero<T>();
 292           for (int i = 0; i < N; ++i) {
 293             ret[col] += lhs[i] * rhs[col][i];
 294           }
 295         }
 296         return ret;
 297       }
 298
 299       inline friend Matrix operator * (const Matrix & lhs, const Matrix & rhs)
 300       {
 301         Matrix ret;
 302         for (int col = 0; col < N; ++col) {
 303           for (int row = 0; row < N; ++row) {
 304             T val = maths::zero<T>();
 305             for (int i = 0; i < N; ++i) {
 306               val += lhs[i][row] * rhs[col][i];
 307             }
 308             ret[col][row] = val;
 309           }
 310         }
 311         return ret;
 312       }
 313       inline Matrix & operator *= (const Matrix & rhs)
 314       {
 315         return (*this = *this * rhs);
 316       }
 317
 318       inline Matrix & operator /= (const T rhs)
 319       {
 320         for (int i = 0; i < N; ++i)
 321           col[i] /= rhs;
 322         return *this;
 323       }
 324       inline Matrix operator / (const T rhs) const
 325       {
 326         Matrix ret = *this;
 327         return ret /= rhs;
 328       }
 329
 330       /// Make the matrix an identity matrix.
 331       inline Matrix & identity()
 332       {
 333         for (int i = 0; i < N; ++i) {
 334           for (int j = 0; j < i; ++j) {
 335             col[i][j] = col[j][i] = maths::zero<T>();
 336           }
 337           col[i][i] = 1;
 338         }
 339         return *this;
 340       }
 341
 342       /// Transpose the matrix.
 343       inline Matrix & transpose()
 344       {
 345         // go through lower left cells and swap with upper right
 346         for (int i = 1; i < N; ++i) {
 347           for (int j = 0; j < i; ++j) {
 348             T temp = col[i][j];
 349             col[i][j] = col[j][i];
 350             col[j][i] = temp;
 351           }
 352         }
 353         return *this;
 354       }
 355
 356       /// Get the transposed matrix.
 357       inline Matrix transposed() const
 358       {
 359         Matrix ret;
 360         for (int i = 0; i < N; ++i) {
 361           for (int j = 0; j < N; ++j) {
 362             ret[j][i] = col[i][j];
 363           }
 364         }
 365         return ret;
 366       }
 367
 368       /// Invert the matrix.
 369       Matrix & invert();
 370
 371       /// Get the inverted matrix.
 372       inline Matrix inverted() const
 373       {
 374         Matrix ret = *this;
 375         return ret.invert();
 376       }
 377   }; // ::maths::Matrix
 378
 379   inline Matrix<4,float> ScaleMatrix44(const maths::Vector<3,float> & v)
 380   {
 381     return ScaleMatrix44((v, 1.0f));
 382   }
 383
 384   template <int N, typename T>
 385   inline maths::Matrix<N,T> outer(const maths::Vector<N,T> & vertical, const maths::Vector<N,T> & horizontal)
 386   {
 387     maths::Matrix<N,T> result;
 388     for (int i = 0; i < N; ++i) {
 389       result[i] = horizontal[i] * vertical;
 390     }
 391     return result;
 392   }
 393 } // ::maths
 394
 395 #endif