RTVolumeRendering/Source/Matrix.h

Go to the documentation of this file.

#ifndef MATRIX_H
#define MATRIX_H

#include "Point3.h"
#include "Vector3.h"

template<class T>
class Matrix
{
public:
        
        union
        {
                struct
                {
                        T _11, _12, _13, _14;
                        T _21, _22, _23, _24;
                        T _31, _32, _33, _34;
                        T _41, _42, _43, _44;
                };
                T m[4][4];
        };
        
        Matrix() {}
        
        void MakeZero()
        {
                for(int i = 0; i < 4; i++)
                {
                        for(int j = 0; j < 4; j++)
                        {
                                m[i][j] = (T)0.0;
                        }
                }
        }
        
        void MakeIdentity()
        {
                MakeZero();
                _11 = _22 = _33 = _44 = (T)1.0;
        }

        void Transpose()
        {
                T tmp;
                tmp = _12; _12 = _21; _21 = tmp;
                tmp = _13; _13 = _31; _31 = tmp;
                tmp = _14; _14 = _41; _41 = tmp;
                tmp = _23; _23 = _32; _32 = tmp;
                tmp = _34; _34 = _43; _43 = tmp;
        }
        
        void MakeTranslation(T x, T y, T z)
        {
                MakeIdentity();
                _14 = x;
                _24 = y;
                _34 = z;
        }
        
        void MakeScale(T x, T y, T z)
        {
                MakeIdentity();
                _11 = x;
                _22 = y;
                _33 = z;
        }

        // NOTE: Expects angle in radians
        void MakeRotationX(T angle)
        {
                const T s = std::sin(angle);
                const T c = std::cos(angle);
                
                MakeIdentity();
                _22 = c; _23 = -s;
                _32 = s; _33 = c;
        }

        // NOTE: Expects angle in radians
        void MakeRotationY(T angle)
        {
                const T s = std::sin(angle);
                const T c = std::cos(angle);
                
                MakeIdentity();
                _11 = c; _13 = s;
                _31 = -s; _33 = c;
        }
        
        // NOTE: Expects angle in radians
        void MakeRotationZ(T angle)
        {
                const T s = std::sin(angle);
                const T c = std::cos(angle);
                
                MakeIdentity();
                _11 = c; _12 = -s;
                _21 = s; _22 = c;
        }


        // NOTE: Expects angle in radians
        void MakeRotation(const Vector3<T> axis, T angle)
        {
                // Normalizes x y and z
                const T len = axis.Length();
                const T x = axis.x / len;
                const T y = axis.y / len;
                const T z = axis.z / len;

                const T c = std::cos(angle);
                const T s = std::sin(angle);
                const T t = 1 - c;
                
                MakeIdentity();

                _11 = t*x*x + c;
                _12 = t*x*y - s*z;
                _13 = t*x*z + s*y;

                _21 = t*x*y + s*z;
                _22 = t*y*y + c;
                _23 = t*y*z - s*x;

                _31 = t*x*z - s*y;
                _32 = t*y*z + s*x;
                _33 = t*z*z + c;
        }
        
        Point3<T> Transform(const Point3<T> &p) const
        {
                Point3<T> out;
                out.x = p.x*_11 + p.y*_12 + p.z*_13 + _14;
                out.y = p.x*_21 + p.y*_22 + p.z*_23 + _24;
                out.z = p.x*_31 + p.y*_32 + p.z*_33 + _34;
                return out;
        }

        Vector3<T> Transform(const Vector3<T> &v) const
        {
                Vector3<T> out;
                out.x = v.x*_11 + v.y*_12 + v.z*_13;
                out.y = v.x*_21 + v.y*_22 + v.z*_23;
                out.z = v.x*_31 + v.y*_32 + v.z*_33;
                return out;
        }

        // Creates a view transform matrix with the camera centered at the eye directed towards the look position.
        void MakeView(const Point3<T> eye, const Point3<T> lookAt, const Vector3<T> up)
        {
                Vector3<T> zaxis = Normalize(lookAt - eye);
                Vector3<T> xaxis = Normalize(zaxis.Cross(up));
                Vector3<T> yaxis = Normalize(xaxis.Cross(zaxis));

                MakeView(eye, xaxis, yaxis, zaxis);
        }

        void MakeView(const Point3<T> &eye, const Vector3<T> &xaxis, const Vector3<T> &yaxis, const Vector3<T> &zaxis)
        {       
                /*
                        NOTE: Produces an equivalent matrix to gluLookAt()
                        
                        View Transform:

                        R = xaxis
                        U = yaxis
                        F = zaxis
                        E = eye

                        . = the dot product
                    
                            [  Rx   Ry   Rz  0 ]   [ 1  0  0  -Ex ]
                    V = [  Ux   Uy   Uz  0 ] * [ 0  1  0  -Ey ]
                            [ -Fx  -Fy  -Fz  0 ]   [ 0  0  1  -Ez ]
                            [   0    0    0  1 ]   [ 0  0  0    1 ]

                            [  Rx   Ry   Rz  -R.E ]
                    V = [  Ux   Uy   Uz  -U.E ]
                            [ -Fx  -Fy  -Fz   F.E ]
                            [   0    0    0     1 ]
                */

                MakeIdentity();
                _11 = xaxis.x;  _12 = xaxis.y;  _13 = xaxis.z;
                _21 = yaxis.x;  _22 = yaxis.y;  _23 = yaxis.z;
                _31 = -zaxis.x; _32 = -zaxis.y; _33 = -zaxis.z;

                _14 = -Dot(xaxis, eye);
                _24 = -Dot(yaxis, eye);
                _34 =  Dot(zaxis, eye);
        }

        // Creates a Perspective transform matrix
        void MakePerspective(T width, T height, T n, T f)
        {
                T aspect = width / height;

                T dw = aspect / ((T)2);
                T dh = ((T)0.5);
                MakePerspective(-dw, dw, -dh, dh, n, f);
        }

        void MakePerspective(T left, T right, T bottom, T top, T n, T f)
        {
                /*
                        NOTE: creates a glFrustrum() equivalent matrix.

                        The prespective transform

                        n = near
                        f = far
                        l = left
                        r = right
                        b = bottom
                        t = top

                        [ 2n/(r-l)    0       A    0 ]
                        [  0       2n/(t-b)   B    0 ]
                        [  0          0      -C   -D ]
                        [  0          0      -1    0 ]

                        A = (r + l) / (r - l)
                        B = (t + b) / (t - b)
                        C = (f + n) / (f - n)
                        D = (2*f*n) / (f - n)

                */

                T A = (right + left) / (right - left);
                T B = (top + bottom) / (top - bottom);
                T C = (f + n) / (f - n);
                T D = (((T)2)*f*n) / (f - n);

                MakeZero();
                _11 = ((T)2)*n / (right - left);
                _22 = ((T)2)*n / (top - bottom);
                _13 = A;
                _23 = B;
                _33 = -C;
                _34 = -D;
                _43 = ((T)-1);
        }

        T Det() const
        {
                T a1 = _11*_22*_33*_44 + _11*_23*_34*_42 + _11*_24*_32*_43;
                T a2 = _12*_21*_34*_43 + _12*_23*_31*_44 + _12*_24*_33*_41; 
                T a3 = _13*_21*_32*_44 + _13*_22*_34*_41 + _13*_24*_31*_42;
                T a4 = _14*_21*_32*_44 + _14*_22*_31*_43 + _14*_23*_32*_41;
                T a5 = _11*_22*_34*_43 + _11*_23*_32*_44 + _11*_24*_33*_42;
                T a6 = _12*_21*_33*_44 + _12*_23*_34*_41 + _12*_24*_31*_43;
                T a7 = _13*_21*_34*_42 + _13*_22*_31*_44 + _13*_24*_32*_41;
                T a8 = _14*_21*_32*_43 + _14*_22*_33*_41 + _14*_23*_31*_42;

                return a1 + a2 + a3 + a4 - a5 - a6 - a7 - a8;
        }

        // Fills m with the openGl equivalent of this matrix
        void GetOpenGLMatrix(T m[16]) const
        {
                // OpenGL stores matrices in column-major order
                m[0] = _11; m[4] = _12; m[8]  = _13; m[12] = _14;
                m[1] = _21; m[5] = _22; m[9]  = _23; m[13] = _24;
                m[2] = _31; m[6] = _32; m[10] = _33; m[14] = _34;
                m[3] = _41; m[7] = _42; m[11] = _43; m[15] = _44;
        }
        
        Matrix operator+(const Matrix &b) const
        {
                Matrix c;
                for(int i = 0; i < 4; i++)
                {
                        for(int j = 0; j < 4; j++)
                        {
                                c.m[i][j] = m[i][j] + b.m[i][j];
                        }
                }
                return c;
        }
        
        Matrix operator*(const Matrix &b) const
        {
                Matrix c;
                for(int i = 0; i < 4; i++)
                {
                        for(int j = 0; j < 4; j++)
                        {
                                T sum = (T)0.0;
                                for(int k = 0; k < 4; k++)
                                {
                                        sum += m[i][k] * b.m[k][j];
                                }
                                c.m[i][j] = sum;
                        }
                }
                return c;
        }
        
        Matrix operator*(T k) const
        {
                Matrix c;
                for(int i = 0; i < 4; i++)
                {
                        for(int j = 0; j < 4; j++)
                        {
                                c.m[i][j] = k * m[i][j];
                        }
                }
                return c;
        }
        
        friend Matrix operator*(T k, const Matrix &a)
        {
                return a*k;
        }
};

typedef Matrix<float> Matrixf;
typedef Matrix<double> Matrixd;

#endif

---