RSAGL/Math/Matrix.lhs

\section{RSAGL.Matrix} \begin{code}


module RSAGL.Math.Matrix
    (Matrix,
     matrix,
     transformHomogenous,
     rowMajorForm,
     colMajorForm,
     rowAt,
     matrixAt,
     identity_matrix,
     translationMatrix,
     rotationMatrix,
     scaleMatrix,
     xyzMatrix,
     matrixAdd,
     matrixMultiply,
     matrixTranspose,
     matrixInverse,
     determinant,
     matrixInversePrim,
     matrixTransposePrim,
     matrixInverseTransposePrim,
     determinantPrim)
    where

import Data.List as List
import RSAGL.Math.Angle
import RSAGL.Math.Vector
import Data.Vec as Vec
import RSAGL.Types

\end{code} A 4-by-4 matrix with cached inverse, transpose, inverse transpose, and determinant. \begin{code}

data Matrix = Matrix { matrix_data :: !(Mat44 RSdouble),
                       matrix_inverse :: Matrix,
                       matrix_transpose :: Matrix,
                       matrix_determinant :: RSdouble }

instance Eq Matrix where
    x == y = matrix_data x == matrix_data y

instance Show Matrix where
    show m = show $ rowMajorForm m

rowMajorForm :: Matrix -> [[RSdouble]]
rowMajorForm = matToLists . matrix_data

colMajorForm :: Matrix -> [[RSdouble]]
colMajorForm = List.transpose . rowMajorForm

\end{code} rowAt answers the nth row of a matrix. \begin{code}

rowAt :: Matrix -> Int -> [RSdouble]
rowAt m n = (rowMajorForm m) !! n

\end{code} matrixAt answers the (i'th,j'th) element of a matrix. \begin{code}

matrixAt :: Matrix -> (Int,Int) -> RSdouble
matrixAt m (i,j) = rowAt m i !! j

\end{code} \subsection{Constructing matrices} matrix constructs a matrix from row major list form. (Such a list form can be formatted correctly in monospaced font and haskell syntax, so that it looks like a matrix as it would be normally written.) \begin{code}

matrix :: [[RSdouble]] -> Matrix
matrix = uncheckedMatrix . matFromLists

-- | Generate a column matrix of length 4, perform an affine transformation on it, and produce the resulting value.
{-# INLINE transformHomogenous #-}
transformHomogenous :: RSdouble -> RSdouble -> RSdouble -> RSdouble -> (RSdouble -> RSdouble -> RSdouble -> a) -> Matrix -> a
transformHomogenous x y z w f m = f x' y' z'
    where (x':.y':.z':._) = multmv (matrix_data m) (x:.y:.z:.w:.())

uncheckedMatrix :: Mat44 RSdouble -> Matrix
uncheckedMatrix dats = m
    where m_inverse = (matrixInversePrim m) { matrix_inverse = m, matrix_transpose = m_inverse_transpose, matrix_determinant = recip m_det }
          m_transpose = (matrixTransposePrim m) { matrix_inverse = m_inverse_transpose, matrix_transpose = m, matrix_determinant = m_det }
          m_inverse_transpose = (matrixTransposePrim m_inverse) { matrix_inverse = m_transpose, matrix_transpose = m_inverse, matrix_determinant = recip m_det }
          m_det = determinantPrim m
          m = Matrix { matrix_data=dats,
                       matrix_inverse = m_inverse,
                       matrix_transpose = m_transpose,
                       matrix_determinant = m_det }

\end{code} identityMatrix constructs the n by n identity matrix \begin{code}

identity_matrix :: Matrix
identity_matrix = Matrix {
    matrix_data = identity,
    matrix_inverse = identity_matrix,
    matrix_transpose = identity_matrix,
    matrix_determinant = 1 }

\end{code} \begin{code}

translationMatrix :: Vector3D -> Matrix
translationMatrix (Vector3D x y z) = uncheckedMatrix $ (1:.0:.0:.x:.()) :.
                                                       (0:.1:.0:.y:.()) :.
                                                       (0:.0:.1:.z:.()) :.
                                                       (0:.0:.0:.1:.()) :. ()

\end{code} \begin{code}

rotationMatrix :: Vector3D -> Angle -> Matrix
rotationMatrix vector angle = let s = sine angle
				  c = cosine angle
				  c' = 1 - c
				  (Vector3D x y z) = vectorNormalize vector
				  in uncheckedMatrix $ ((c+c'*x*x)   :. (c'*y*x - s*z) :. (c'*z*x + s*y) :. 0 :. ()) :.
					               ((c'*x*y+s*z) :. (c+c'*y*y)     :. (c'*z*y - s*x) :. 0 :. ()) :.
					               ((c'*x*z-s*y) :. (c'*y*z+s*x)   :. (c+c'*z*z)     :. 0 :. ()) :.
					               (0            :. 0              :. 0              :. 1 :. ()) :. ()

\end{code} \begin{code}

scaleMatrix :: Vector3D -> Matrix
scaleMatrix (Vector3D x y z) = uncheckedMatrix $ (x :. 0 :. 0 :. 0) :.
				                 (0 :. y :. 0 :. 0) :.
				                 (0 :. 0 :. z :. 0) :.
				                 (0 :. 0 :. 0 :. 1) :. ()

\end{code} \texttt{xyzMatrix} constructs the matrix in which the x y and z axis are transformed to point in the direction of the specified vectors. \begin{code}

xyzMatrix :: Vector3D -> Vector3D -> Vector3D -> Matrix
xyzMatrix (Vector3D x1 y1 z1) (Vector3D x2 y2 z2) (Vector3D x3 y3 z3) =
    uncheckedMatrix $ (x1 :. x2 :. x3 :. 0) :.
                      (y1 :. y2 :. y3 :. 0) :.
                      (z1 :. z2 :. z3 :. 0) :.
                      (0  :. 0  :. 0  :. 1.0) :. ()

\end{code} \subsection{Matrix arithmetic} \begin{code}

matrixAdd :: Matrix -> Matrix -> Matrix
matrixAdd m n = uncheckedMatrix $ Vec.zipWith (+) (matrix_data m) (matrix_data n)

\end{code} \label{howtoUseMatrixMultiplyImpl} \begin{code}

matrixMultiply :: Matrix -> Matrix -> Matrix
matrixMultiply m n = uncheckedMatrix $ multmm (matrix_data m) (matrix_data n)

matrixTranspose :: Matrix -> Matrix
matrixTranspose = matrix_transpose

\end{code} \subsection{Inverse and determinant of a matrix} \begin{code}

matrixInverse :: Matrix -> Matrix
matrixInverse = matrix_inverse

\end{code} \begin{code}

determinant :: Matrix -> RSdouble
determinant = matrix_determinant

\end{code} \begin{code}

matrixInverseTransposePrim :: Matrix -> Matrix
matrixInverseTransposePrim = uncheckedMatrix . Vec.transpose . fst . invertAndDet . matrix_data

matrixInversePrim :: Matrix -> Matrix
matrixInversePrim = uncheckedMatrix . fst . invertAndDet . matrix_data

determinantPrim :: Matrix -> RSdouble
determinantPrim = det . matrix_data

matrixTransposePrim :: Matrix -> Matrix
matrixTransposePrim = uncheckedMatrix . Vec.transpose . matrix_data

\end{code}