Matrix And Quaternions Faq

~cpp 
  In this document (as in most math textbooks), all matrices are drawn
  in the standard mathematical manner.  Unfortunately graphics libraries
  like IrisGL, OpenGL and SGI's Performer all represent them with the
  rows and columns swapped.

  Hence, in this document you will see (for example) a 4x4 Translation
  matrix represented as follows:

          | 1  0  0  X |
          |            |
          | 0  1  0  Y |
      M = |            |
          | 0  0  1  Z |
          |            |
          | 0  0  0  1 |

  In Performer (for example) this would be populated as follows:

    M[0][1] = M[0][2] = M[0][3] =
    M[1][0] = M[1][2] = M[1][3] =
    M[2][0] = M[2][1] = M[2][3] = 0 ;
    M[0][0] = M[1][1] = M[2][2] = m[3][3] = 1 ;
    M[3][0] = X ;
    M[3][1] = Y ;
    M[3][2] = Z ;

  ie, the matrix is stored like this:

          | M[0][0]  M[1][0]  M[2][0]  M[3][0] |
          |                                    |
          | M[0][1]  M[1][1]  M[2][1]  M[3][1] |
      M = |                                    |
          | M[0][2]  M[1][2]  M[2][2]  M[3][2] |
          |                                    |
          | M[0][3]  M[1][3]  M[2][3]  M[3][3] |

  OpenGL uses a one-dimensional array to store matrices - but fortunately,
  the packing order results in the same layout of bytes in memory - so
  taking the address of a pfMatrix and casting it to a float* will allow
  you to pass it directly into routines like glLoadMatrixf.

  In the code snippets scattered throughout this document, a one-dimensional
  array is used to store a matrix. The ordering of the array elements is
  transposed with respect to OpenGL.

         This Document                  OpenGL

        | 0  1  2  3  |            | 0  4  8  12 |
        |             |            |             |
        | 4  5  6  7  |            | 1  5  9  13 |
    M = |             |        M = |             |
        | 8  9  10 11 |            | 2  6  10 14 |
        |             |            |             |
        | 12 13 14 15 |            | 3  7  11 15 |

~cpp 

  A matrix is a two dimensional array of numeric data, where each
  row or column consists of one or more numeric values.

  Arithmetic operations which can be performed with matrices include
  addition, subtraction, multiplication and division.

  The size of a matrix is defined in terms of the number of rows
  and columns.

  A matrix with M rows and N columns is defined as a MxN matrix.

  Individual elements of the matrix are referenced using two index
  values. Using mathematical notation these are usually assigned the
  variables 'i' and 'j'. The order is row first, column second

  For example, if a matrix M with order 4x4 exists, then the elements
  of the matrix are indexed by the following row:column pairs:

        | 00 10 20 30 |
    M = | 01 11 21 31 |
        | 02 12 22 32 |
        | 03 13 23 33 |

  The element at the top right of the matrix has i=0 and j=3
  This is referenced as follows:

    M    = M
     i,j    0,3

  In computer animation, the most commonly used matrices have either
  2, 3 or 4 rows and columns. These are referred to as 2x2, 3x3 and 4x4
  matrices respectively.

  2x2 matrices are used to perform rotations, shears and other types
  of image processing. General purpose NxN matrices can be used to
  perform image processing functions such as convolution.

  3x3 matrices are used to perform low-budget 3D animation. Operations
  such as rotation and multiplication can be performed using matrix
  operations, but perspective depth projection is performed using
  standard optimised into pure divide operations.

  4x4 matrices are used to perform high-end 3D animation. Operations
  such as multiplication and perspective depth projection can be
  performed using matrix mathematics.

~cpp 
  The "order" of a matrix is another name for the size of the matrix.
  A matrix with M rows and N columns is said to have order MxN.

~cpp 
  The simplest way of defining a matrix using the C/C++ programming
  languages is to make use of the "typedef" keyword. Both 3x3 and 4x4
  matrices may be defined in this way ie:

    typedef float MATRIX3[9];
    typedef float MATRIX4[16];

  Since each type of matrix has dimensions 3x3 and 4x4, this requires
  9 and 16 data elements respectively.

  At first glance, the use of a single linear array of data values may
  seem counter-intuitive. The use of two dimensional arrays may seem
  more convenient ie.

    typedef float MATRIX3[3][3];
    typedef float MATRIX4[4][4];

  However, the use of two reference systems for each matrix element
  very often leads to confusion. With mathemetics, the order is row
  first (i), column second (j) ie.

     Mij

  Using C/C++, this becomes

     matrix[j][i]

  Using two dimensional arrays also incurs a CPU performance penalty in
  that C compilers will often make use of multiplication operations to
  resolve array index operations.

  So, it is more efficient to stick with linear arrays. However, one issue
  still remains to be resolved. How is an two dimensional matrix mapped
  onto a linear array? Since there are only two methods (row first/column
  second or column first/row column).

  The performance differences between the two are subtle. If all for-next
  loops are unravelled, then there is very little difference in the
  performance for operations such as matrix-matrix multiplication.

  Using the C/C++ programming languages the linear ordering of each
  matrix is as follows:

  mat[0]  = M        mat[3]  = M
             00                 03

  mat[12] = M        mat[15] = M
             30                 33

        |  0  1  2  3 |
        |             |              | 0 1 2 |
        |  4  5  6  7 |              |       |
    M = |             |          M = | 3 4 5 |
        |  8  9 10 11 |              |       |
        |             |              | 6 7 8 |
        | 12 13 14 15 |

~cpp 
  One of the first questions asked about the use of matrices in computer
  animation is why they should be used at all in the first place.
  Intuitively, it would appear that the overhead of for-next loops and
  matrix multiplication would slow down an application.

  Arguments that resolve these objections can be pointed out. These include
  the use of CPU registers to handle loop counters on-board data caches
  to optimise memory accesses.

  Advantages can also be pointed out. By following a mathematical approach
  to defining 3D algorithms, it is possible to predict and plan the
  design of a 3D animation system. Such mathematical approaches allow
  for the implementation of character animation, spline curves and inverse
  kinematics.

  However, one objection that frequently comes up is that it would be
  quicker to just multiply each pair of coordinates by the rotation
  coefficients for that axis, rather than perform a full vector-matrix
  multiplication.

  ie. Rotation in X transforms Y and Z
      Rotation in Y transforms X and Z
      Rotation in Z transforms X and Y

  The argument to this goes as follows:

  Given a vertex V = (x,y,z), rotation angles (A,B and C) and translation
  (D,E,F). A  the algorithm
  is defined as follows:

    ---------------------------

    sx = sin(A)             // Setup - only done once
    cx = cos(A)
    sy = sin(B)
    cy = cos(B)
    sz = sin(C)
    cz = cos(C)

    x1 =  x * cz +  y * sz  // Rotation of each vertex
    y1 =  y * cz -  x * sz
    z1 =  z

    x2 = x1 * cy + z1 * sy
    y2 = z1
    z2 = z1 * cy - x1 * sy

    x3 = x2
    y3 = y2 * cx + z1 * sx
    z3 = z2 * cx - x1 * sx

    xr = x3 + D             // Translation of each vertex
    yr = y3 + E
    zr = z3 + F

    ---------------------------

  Altogether, this algorithm will use the following amounts of processing
  time:

    Set-up                                 Per-vertex
    -------------------------              ------------------------
    6 trigonometric functions
    6 assignment operations.               12 assignment
                                           12 multiplication
                                            9 addition
    -------------------------              ------------------------

  Assume that the same operations is being performed using matrix
  multiplication.

  With a 4x4 matrix, the procesing time is used as follows:

    Set-up                       Change    Per-vertex               Change
    --------------------------   ------    ------------------------ ------
    6  trigonometric functions    0                                  0
    18 assignment operation      -12        3  assignment           -9
    12 multiplication            +12        9  multiplication       -3
    6  subtraction               +6         6  addition             -3
    --------------------------   ------    ------------------------ ------

  Comparing the two tables, it can be seen that setting up a rotation
  matrix costs at least 12 multiplication calculations and an extra
  18 assignment calls.

  However, while this may seem extravagant, the savings come from
  processing each vertex. Using matrix multiplication, the savings made
  from processing just 4 vertices, will outweigh the additional set-up
  cost.

~cpp 
  With either 3x3 or 4x4 rotation, translation or shearing matrices, there
  is a simple relationship between each matrix and the resulting coordinate
  system.

  The first three columns of the matrix define the direction vector of the
  X, Y and Z axii respectively.

  If a 4x4 matrix is defined as:

        | A B C D |
    M = | E F G H |
        | I J K L |
        | M N O P |

  Then the direction vector for each axis is as follows:

     X-axis = [ A E I ]

     Y-axis = [ B F J ]

     Z-axis = [ C G K ]

~cpp 
  The identity matrix is matrix in which has an identical number of rows
  and columns. Also, all the elements in which i=j are set one. All others
  are set to zero. For example a 4x4 identity matrix is as follows:

        | 1 0 0 0 |
    M = | 0 1 0 0 |
        | 0 0 1 0 |
        | 0 0 0 1 |

~cpp 
  The major diagonal of a matrix is the set of elements where the
  row number is equal to the column number ie.

    M   where i=j
     ij

  In the case of the identity matrix, only the elements on the major
  diagonal are set to 1, while all others are set to 0.

~cpp 
  The transpose of matrix is the matrix generated when every element in
  the matrix is swapped with the opposite relative to the major diagonal

  This can be expressed as the mathematical operation:

    M'   = M
      ij    ji

  However, this can only be performed if a matrix has an equal number
  of rows and columns.

  If the matrix M is defined as:

        |  0.707 -0.866 |
    M = |               |
        |  0.866  0.707 |

  Then the transpose is equal to:

        |  0.707  0.866 |
    T = |               |
        | -0.866  0.707 |

  If the matrix is a rotation matrix, then the transpose is guaranteed
  to be the inverse of the matrix.

~cpp 
  The rule of thumb with adding two matrices together is:

    "add row and column to row and column"

  This can be expressed mathematically as:

    R   = M   + L
     ij    ij    ij

  However, both matrices must be identical in size.

  For example, if the 2x2 matrix M is added with the 2x2 matrix L then
  the result is as follow:

    R = M + L

        | A B C |   | J K L |
        |       |   |       |
      = | D E F | + | M N O |
        |       |   |       |
        | G H I |   | P Q R |

        | A+J B+K C+L |
        |             |
      = | D+M E+N F+O |
        |             |
        | G+P H+Q I+R |

~cpp 
  The rule of thumb with subtracting two matrices is:

    "subtract row and column from row and column"

  This can be expressed mathematically as:

    R   = M   - L
     ij    ij    ij

  However, both matrices must be identical in size.

  For example, if the 2x2 matrix L is subtracted from the 2x2 matrix M then
  the result is as follows:

    R = M - L

        | A B C |   | J K L |
        |       |   |       |
      = | D E F | - | M N O |
        |       |   |       |
        | G H I |   | P Q R |

        | A-J B-K C-L |
        |             |
      = | D-M E-N F-O |
        |             |
        | G-P H-Q I-R |

~cpp 
  The rule of thumb with multiplying two matrices together is:

    "multiply row into column and sum the result".

  This can be expressed mathematically as:

          n
          --
    R   = \   M   x L
     ij   /    ij    ji
          --
          i=1

  If the two matrices to be multiplied together have orders:

    M = AxB and L = CxD

  then the two values B and C must be identical.

  Also, the resulting matrix has an order of AxD

  Thus, it is possible to multiply a 4xN matrix with a 4x4 matrix
  but not the other way around.

  For example, if the 4x4 matrix M is defined as:

        | A B C D |
    M = | E F G H |
        | I J K L |
        | M N O P |

  and a 4x2 matrix L is defined as:

    L = | Q R |
        | S T |
        | U V |
        | W X |

  then the size of the resulting matrix is 2x4. The resulting matrix
  is defined as:

    R = M x L

        | A B C D |   | Q R |
      = | E F G H | x | S T |
        | I J K L |   | U V |
        | M N O P |   | W X |

        | AQ+BS+CU+DW  AR+BT+CV+DX |
      = | EQ+FS+GU+HW  ER+FT+GV+HX |
        | IQ+JS+KU+LW  IR+JT+KV+LX |
        | MQ+NS+OU+PW  MR+NT+OV+PX |

~cpp 
  A matrix may be squared or even raised to an integer power. However
  there are several restrictions. For all powers, the matrix must be
  orthogonal ie. have the same width and height

  For example,

     -1
    M   is the inverse of the matrix

     0
    M   generates the identity matrix

     1
    M   leaves the matrix undamaged.

     2
    M   squares the matrix and

     3
    M   generates the cube of the matrix

  Raising a matrix to a power greater than one involves multiplying a matrix
  by itself a specific number of times.

  For example,

     2
    M  = M . M

     3
    M  = M . M . M

  and so on.

  Raising the identity matrix to any power always generates the identity
  matrix ie.

     n
    I  = I

~cpp 
  The best way to perform this task is to treat the list of vectors as
  a single matrix, with each vector represented as a column vector.

  If N vectors are to be multiplied by a 4x4 matrix, then they can be
  treated as a single 4xN matrix:

  If the matrix is defined as:

        | A B C D |
    M = | E F G H |
        | I J K L |
        | M N O P |

  and the list of vectors is defined as:

        | x1 x2 x3 x4 x5|
    V = | y1 y2 y3 y4 y5|
        | z1 z2 z3 z4 z5|
        | 1  1  1  1    |

  Note that an additional row of constant terms is added to the vector
  list, all of which are set to 1.0. In real life, this row does not
  exist. It is simply used to make the orders of the matrix M and the
  vector list V match.

  Then the multiplication is performed as follows:

            M . V = V'

  | A B C D |   | x1 x2 x3 x4 x5 |   | A.x1+B.y1+C.z1+D A.x2+B.y2+C.z2+D ... |
  | E F G H | . | y1 y2 y3 y4 y5 | = | E.x1+F.y1+G.z1+H E.x2+F.y2+G.z2+H ... |
  | I J K L |   | z1 z2 y3 y4 z5 |   | I.x1+J.y1+K.z1+L I.x2+J.y2+K.z2+L ... |
  | M N O P |   | 1  1  1  1  1  |   | M.x1+N.y1+O.z1+P M.x2+N.y2+O.z2+P ... |

  For each vector in the list there will be a total of 12 multiplication
  16 addition and 1 division operation (for perspective).

  If the matrix is known not to be a rotation or translation matrix then the
  division operation can be skipped.

~cpp 
  The determinant of a matrix is a floating point value which is used to
  indicate whether the matrix has an inverse or not. If negative, then
  no inverse exists. If the determinant is positive, then an inverse
  exists.

  For an identity matrix, the determinant is always equal to one.
  Any matrix with a determinant of 1.0 is said to be isotropic.

  Thus all rotation matrices are said to be isotropic, since the
  determinant is always equal to 1.0.

  This can be proved as follows:

        | A B |   | cos X  -sin X |
    M = |     | = |               |
        | C D |   | sin X   cos X |

    D = AD - BC

    D = (cos X . cos X) -  (-sin X . sin X)

            2          2
    D = (cos X ) + (sin X)

            2       2
    But, cos X + sin X = 1

    Therefore,

    D = 1

~cpp 
  The determinant of a matrix is calculated using Kramer's rule, where
  the value can be calculated by breaking the matrix into smaller
  matrices.

  For a 2x2 matrix M, the determinant D is calculated as follows:

        | A B |
    M = |     |
        | C D |

    D = AD - BC

  For 3x3 and 4x4 matrices, this is more complicated, but can be solved
  by methods such as Kramer's Rule.

~cpp 
  An Isotropic matrix is one in which the sum of the squares of all
  three rows or columns add up to one.

  A matrix in which this is not the case, is said to be Anisotropic.

  When 3x3 or 4x4 matrices are used to rotate and scale an object, it
  is sometimes necessary to enlarge or shrink one axis more than the
  others.

  For example, with seismic surveys, it is convenient to enlarge the
  Z-axis by a factor or 50 or more, while letting the X and Y axii
  remain the same.

  Another example is the implementation of "squash" and "stretch"
  with character animation. When a character is hit by a heavy object
  eg. an anvil, the desired effect is to character stretched out
  sideways and squashed vertically:

  A suitable matrix would be as follows:

        |  2   0   0    0  |
    M = |  0   2   0    0  |
        |  0   0   0.5  0  |
        |  0   0   0    1  |

  However, there is problem looming ahead. While this matrix will cause
  no problems with the transformation of vertex data, it will cause
  problems with gouraud shading using outward normals.

  Because the transformation stage is implemented using matrix
  multiplication, both vertex data and outward normal data will be
  multiplied with this matrix.

  While this is not a problem with vertex data (it is the desired effect)
  it causes a major headache with the outward normal data.

  After raw multiplication, each outward normal will no longer be
  normalised and consequently will affect other calculations such as
  shading and back-face culling.

~cpp                                                                  -1
  Given a matrix M, then the inverse of that matrix, denoted as M  , is
  the matrix which satisfies the following expression:

         -1
    M . M   = I

  where I is the identity matrix.

  Thus, multiplying a matrix with its inverse will generate the identity
  matrix. However, several requirements must be satisfied before the
  inverse of a matrix can be calculated.

  These include that the width and height of the matrix are identical and
  that the determinant of the matrix is non-zero.

  Calculating the inverse of a matrix is a task often performed in order
  to implement inverse kinematics using spline curves.

~cpp 
  Depending upon the size of the matrix, the calculation of the inverse
  can be trivial or extremely complicated.

  For example, the inverse of a 1x1 matrix is simply the reciprical of
  the single element:

  ie. M = | x |

  Then the inverse is defined as:

     -1   | 1 |
    M   = | - |
          | x |

  Solving 2x2 matrices and larger can be achieved by using Kramer's Rule
  or by solving as a set of simultaneous equations.

  However, in certain cases, such as identity or rotation matrices, the
  inverse is already known or can be determined from taking the transpose
  of the matrix.

~cpp 
  Don't even bother. The inverse of an identity matrix is the identity
  matrix. ie.

         -1
    I . I   = I

  Any identity matrix will always have a determinant of +1.

~cpp 
  Since a rotation matrix always generates a determinant of +1,
  calculating the inverse is equivalent of calculating the transpose.

  Alternatively, if the rotation angle is known, then the rotation
  angle can be negated and used to calculate a new rotation matrix.

~cpp 
  Given a 3x3 matrix M:

        | A B C |
        |       |
    M = | D E F |
        |       |
        | G H I |

  Then the determinant is calculated as follows:

             n
            ---
            \                           i
    det M = /   M    * submat    M  * -1
            ---  0,i         0,i
            i=1

  where

    submat   M defines the matrix composed of all rows and columns of M
          ij

  excluding row i and column j. submat   may be called recursively.
                                      ij

  If the determinant is non-zero then the inverse of the matrix exists.
  In this case, the value of each matrix element is defined by:

     -1      1                            i+j
    M    = -----  *  det submat     M * -1
     j,i   det M               i,j

~cpp 
  For a 2x2 matrix, the calculation is slightly harder. If the matrix is
  defined as follows:

        | A B |
    M = |     |
        | C D |

  Then the determinant is defined as:

    det = AD - BC

  And the inverse is defined as:

     -1     1   |  D  -B  |
    M   =  ---  |         |
           det  | -C   A  |

  This can be proved using Kramer's rule. Given the matrix M:

        | A B |
    M = |     |
        | C D |

  Then the determinant is:

                           0                            1
    det =  M    * submat M    * -1 +  M    * submat M    * -1
            0,0           0,0          0,1           0,1

    <=>    M    * M    * 1         +  M    * M    * -1
            0,0    1,1                 0,1    1,0

    <=>    A  * D                  +  B    * C    * -1

    <=>    AD                      +  BC . -1

    <=>    AD - BC

    ==============

  And the inverse is derived from:

     -1                      0+0      -1
    M    = det submat    * -1    <=> M    = M    *  1 <=> D
     0,0             0,0              0,0    1,1

     -1                      1+0      -1
    M    = det submat    * -1    <=> M    = M    * -1 <=> C * -1
     0,1             1,0              0,1    1,0

     -1                      0+1      -1
    M    = det submat    * -1    <=> M    = M    * -1 <=> B * -1
     1,0             0,1              1,0    0,1

     -1                      1+1      -1
    M    = det submat    * -1    <=> M    = M    *  1 <=> A
     1,1             1,1              1,1    0,0

  Then the inverse matrix is equal to:

     -1    1  |  D  -C |
    M   = --- |        |
          det | -B   A |

  Providing that the determinant is not zero.

~cpp 
  For 3x3 matrices and larger, the inverse can be calculated by
  either applying Kramer's rule or by solving as a set of linear
  equations.

  If Kramer's rule is applied to a matrix M:

        | A B C |
    M = | D E F |
        | G H I |

  then the determinant is calculated as follows:

    det M = A * (EI - HF) - B * (DI - GF) + C * (DH - GE)

  Providing that the determinant is non-zero, then the inverse is
  calculated as:

     -1     1     |   EI-FH  -(BI-HC)   BF-EC  |
    M   = ----- . | -(DI-FG)   AI-GC  -(AF-DC) |
          det M   |   DH-GE  -(AH-GB)   AE-BD  |

  This can be implemented using a pair of 'C' functions:

    ---------------------------------

    VFLOAT m3_det( MATRIX3 mat )
      {
      VFLOAT det;

      det = mat[0] * ( mat[4]*mat[8] - mat[7]*mat[5] )
          - mat[1] * ( mat[3]*mat[8] - mat[6]*mat[5] )
          + mat[2] * ( mat[3]*mat[7] - mat[6]*mat[4] );

      return( det );
      }

   ----------------------------------

   void m3_inverse( MATRIX3 mr, MATRIX3 ma )
     {
     VFLOAT det = m3_det( ma );

     if ( fabs( det ) < 0.0005 )
       {
       m3_identity( ma );
       return;
       }

     mr[0] =    ma[4]*ma[8] - ma[5]*ma[7]   / det;
     mr[1] = -( ma[1]*ma[8] - ma[7]*ma[2] ) / det;
     mr[2] =    ma[1]*ma[5] - ma[4]*ma[2]   / det;

     mr[3] = -( ma[3]*ma[8] - ma[5]*ma[6] ) / det;
     mr[4] =    ma[0]*ma[8] - ma[6]*ma[2]   / det;
     mr[5] = -( ma[0]*ma[5] - ma[3]*ma[2] ) / det;

     mr[6] =    ma[3]*ma[7] - ma[6]*ma[4]   / det;
     mr[7] = -( ma[0]*ma[7] - ma[6]*ma[1] ) / det;
     mr[8] =    ma[0]*ma[4] - ma[1]*ma[3]   / det;
     }

    ---------------------------------

~cpp 
  As with 3x3 matrices, either Kramer's rule can be applied or the
  matrix can be solved as a set of linear equations.

  An efficient way is to make use of the existing 'C' functions defined
  to calculate the determinant and inverse of a 3x3 matrix.

  In order to implement Kramer's rule with 4x4 matrices, it is necessary
  to determine individual sub-matrices. This is achieved by the following
  routine:

    --------------------------
    void m4_submat( MATRIX4 mr, MATRIX3 mb, int i, int j )
      {
      int ti, tj, idst, jdst;

      for ( ti = 0; ti < 4; ti++ )
        {
        if ( ti < i )
          idst = ti;
        else
          if ( ti > i )
            idst = ti-1;

        for ( tj = 0; tj < 4; tj++ )
          {
          if ( tj < j )
            jdst = tj;
          else
            if ( tj > j )
              jdst = tj-1;

          if ( ti != i && tj != j )
            mb[idst*3 + jdst] = mr[ti*4 + tj ];
          }
        }
      }

    --------------------------

  The determinant of a 4x4 matrix can be calculated as follows:

    --------------------------

    VFLOAT m4_det( MATRIX4 mr )
      {
      VFLOAT  det, result = 0, i = 1;
      MATRIX3 msub3;
      int     n;

      for ( n = 0; n < 4; n++, i *= -1 )
        {
        m4_submat( mr, msub3, 0, n );

        det     = m3_det( msub3 );
        result += mr[n] * det * i;
        }

      return( result );
      }

    --------------------------

  And the inverse can be calculated as follows:

    --------------------------

    int m4_inverse( MATRIX4 mr, MATRIX4 ma )
      {
      VFLOAT  mdet = m4_det( ma );
      MATRIX3 mtemp;
      int     i, j, sign;

      if ( fabs( mdet ) < 0.0005 )
        return( 0 );

      for ( i = 0; i < 4; i++ )
        for ( j = 0; j < 4; j++ )
          {
          sign = 1 - ( (i +j) % 2 ) * 2;

          m4_submat( ma, mtemp, i, j );

          mr[i+j*4] = ( m3_det( mtemp ) * sign ) / mdet;
          }

      return( 1 );
      }

    --------------------------

  Having a function that can calculate the inverse of any 4x4 matrix is
  an incredibly useful tool. Application include being able to calculate
  the base matrix for splines, inverse rotations and rearranging matrix
  equations.

~cpp 
  If a matrix M exists, such that:

         | A B C |
    M  = | D E F |
         | G H I |

  then the inverse exists:

         | P Q R |
    M' = | S T U |
         | V W X |

  and the following expression is valid:

                     -1
        M     .     M     = I

    | A B C |   | P Q R |   | 1 0 0 |
    | D E F | . | S T U | = | 0 1 0 |
    | G H I |   | V W X |   | 0 0 1 |

  The inverse can then be calculated through the solution as a set of
  linear equations ie.:

    | AP + BS + CV |   | 1 |   Column 0 (X)
    | DP + ES + FV | = | 0 |
    | GP + HS + IV |   | 0 |

    | AQ + BT + CW |   | 0 |   Column 1 (Y)
    | DQ + ET + FW | = | 1 |
    | GQ + HT + IW |   | 0 |

    | AR + BU + CX |   | 0 |   Column 2 (Z)
    | DR + EU + FX | = | 0 |
    | GR + HU + IX |   | 1 |