root/trunk/helix/linalgebra.d

/*
Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions
are met:

    Redistributions of source code must retain the above copyright
    notice, this list of conditions and the following disclaimer.

    Redistributions in binary form must reproduce the above
    copyright notice, this list of conditions and the following
    disclaimer in the documentation and/or other materials provided
    with the distribution.

    Neither name of Victor Nakoryakov nor the names of
    its contributors may be used to endorse or promote products
    derived from this software without specific prior written
    permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
REGENTS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION)
HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT,
STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
OF THE POSSIBILITY OF SUCH DAMAGE.

Copyright (C) 2006. Victor Nakoryakov.
*/
/**
Module consists of basic mathematical objects oriented to working with 3D
graphics.

Those are 2,3,4-D vectors, quaternion, 3x3 and 4x4 matrices. In case of
specialization for 3D graphics there are always some features and deviation
from classical math. Here I summarize such features of helix'es linear algebra:
$(UL
    $(LI In helix paradigm of column-vector is taken. So multiplication of matrix
         by vector makes sense but multiplication of vector by matrix makes not.
         This approach conforms to rules accepted in classical math and coincides
         with OpenGL rules. However this is opposite to Direct3D paradigm where
         vector is a row. So, in helix, to combine sequence of transforms specified with
         matrices A, B, C in order A then B then C, you have to multiply them in
         back-to-front order order: M=C*B*A. )

    $(LI When an issue deal with euler angles following definitions are accepted.
         Yaw is rotation around Y, Pitch is rotaion around X, Roll is rotation around Z
         Rotations are always applied in order: Roll then Pitch then Yaw. )

    $(LI Helix matrices use column-major memory layout. I.e. matrix
        $(MAT33
            $(MAT33_ROW a, b, c)
            $(MAT33_ROW d, e, f)
            $(MAT33_ROW g, h, i)
        )
        in memory will looks like: a, d, g, b, e, h, c, f, i. This order is the same as
        in OpenGL API, but opposite to Direct3D API. However as mentioned above, Direct3D
        uses vector-row paradigm that is opposite to classical math, so D3D requires
        transposed matrix as compared to classical math to get desired transformation. As a
        result you haven't to transpose helix matrix while transmission to API even in Direct3D
        case. Normaly you haven't to remember about memory layout, just use it as in classical
        math, this feature is significant only in routines that operate with data pointers
        and plain array representation. There are reminders in such methods' documentation. )
)

Authors:
    Victor Nakoryakov, nail-mail[at]mail.ru
    
Macros:
    MAT33     = <table style="border-left: double 3px #666666; border-right: double 3px #666666;
                 margin-left: 3em;">$0</table>
    MAT33_ROW = <tr><td>$1</td><td>$2</td><td>$3</td></tr>
*/
module helix.linalgebra;

private import helix.basic,
               helix.config;

/** Defines ort names that are usualy used as indices. */
enum Ort
{
    X, ///
    Y, /// ditto
    Z, /// ditto
    W  /// ditto
}

/**
Wrapper template to provide possibility to use different float types
in implemented structs and routines.
*/
private template LinearAlgebra(float_t)
{
    private alias helix.basic.equal          equal;
    
    alias .LinearAlgebra!(float).Vector2     Vector2f;
    alias .LinearAlgebra!(float).Vector3     Vector3f;
    alias .LinearAlgebra!(float).Vector4     Vector4f;
    alias .LinearAlgebra!(float).Quaternion  Quaternionf;
    alias .LinearAlgebra!(float).Matrix33    Matrix33f;
    alias .LinearAlgebra!(float).Matrix44    Matrix44f;
    
    alias .LinearAlgebra!(double).Vector2    Vector2d;
    alias .LinearAlgebra!(double).Vector3    Vector3d;
    alias .LinearAlgebra!(double).Vector4    Vector4d;
    alias .LinearAlgebra!(double).Quaternion Quaterniond;
    alias .LinearAlgebra!(double).Matrix33   Matrix33d;
    alias .LinearAlgebra!(double).Matrix44   Matrix44d;
    
    alias .LinearAlgebra!(real).Vector2      Vector2r;
    alias .LinearAlgebra!(real).Vector3      Vector3r;
    alias .LinearAlgebra!(real).Vector4      Vector4r;
    alias .LinearAlgebra!(real).Quaternion   Quaternionr;
    alias .LinearAlgebra!(real).Matrix33     Matrix33r;
    alias .LinearAlgebra!(real).Matrix44     Matrix44r;

    /************************************************************************************
    Two dimensional vector.

    For formal definition of vector, meaning of possible operations and related
    information see $(LINK http://en.wikipedia.org/wiki/Vector_(spatial)).
    *************************************************************************************/
    struct Vector2
    {
        align(1)
        {
            float_t x; /// Components of vector.
            float_t y; /// ditto
        }
        
        static Vector2 nan = { float_t.nan, float_t.nan }; /// Vector with both components seted to NaN.
        static Vector2 unitX = { 1, 0 };                   /// Unit vector codirectional with corresponding axis.
        static Vector2 unitY = { 0, 1 };                   /// ditto
        
        /**
        Method to construct vector in C-like syntax.

        Examples:
        ------------
        Vector2 myVector = Vector2(1, 2);
        ------------
        */
        static Vector2 opCall(float_t x, float_t y)
        {
            Vector2 v;
            v.set(x, y);
            return v;
        }
        
        /** Sets values of components to passed values. */
        void set(float_t x, float_t y)
        {
            this.x = x;
            this.y = y;
        }
        
        /** Returns: Norm (also known as length, magnitude) of vector. */
        real norm()
        {
            return hypot(x, y);
        }
    
        /**
        Returns: Square of vector's norm.
        
        Since this method doesn't need calculation of square root it is better
        to use it instead of norm() when you can. For example, if you want just
        to know which of 2 vectors is longer it's better to compare their norm
        squares instead of their norm.
        */
        real normSquare()
        {
            return x*x + y*y;
        }
    
        /** Normalizes this vector. */
        void normalize()
        {
            *this /= norm;
        }
    
        /** Returns: Normalized copy of this vector. */
        Vector2 normalized()
        {
            real n = norm;
            return Vector2(x / n, y / n);
        }
    
        /**
        Returns: Whether this vector is unit.
        Params:
            relprec, absprec = Parameters passed to equal function while comparison of
                               norm square and 1. Have the same meaning as in equal function.
        */
        bool isUnit(int relprec = defrelprec, int absprec = defabsprec)
        {
            return equal( normSquare(), 1, relprec, absprec );
        }
    
        /** Returns: Axis for which projection of this vector on it will be longest. */
        Ort dominatingAxis()
        {
            return (x > y) ? Ort.X : Ort.Y;
        }
    
        /** Returns: Whether all components are normalized numbers. */
        bool isnormal()
        {
            return std.math.isnormal(x) && std.math.isnormal(y);
        }
    
        /** Returns: float_t pointer to x component of this vector. It's like a _ptr method for arrays. */
        float_t* ptr()
        {
            return cast(float_t*)this;
        }
    
        /** Returns: Component corresponded to passed index. */
        float_t opIndex(Ort ort)
        in { assert(ort <= Ort.Y); }
        body
        {
            return ptr[cast(int)ort];
        }
    
        /** Assigns new _value to component corresponded to passed index. */
        void opIndexAssign(float_t value, Ort ort)
        in { assert(ort <= Ort.Y); }
        body
        {
            ptr[cast(int)ort] = value;
        }
    
        /**
        Standard operators that have intuitive meaning, same as in classical math.
        
        Note that division operators do no cheks of value of k, so in case of division
        by 0 result vector will have infinity components. You can check this with isnormal()
        method.
        */
        bool opEquals(Vector2 v)
        {
            return x == v.x && y == v.y;
        }
    
        /** ditto */
        Vector2 opNeg()
        {
            return Vector2(-x, -y);
        }
    
        /** ditto */
        Vector2 opAdd(Vector2 v)
        {
            return Vector2(x + v.x, y + v.y);
        }
    
        /** ditto */
        void opAddAssign(Vector2 v)
        {
            x += v.x;
            y += v.y;
        }
    
        /** ditto */
        Vector2 opSub(Vector2 v)
        {
            return Vector2(x - v.x, y - v.y);
        }
    
        /** ditto */
        void opSubAssign(Vector2 v)
        {
            x -= v.x;
            y -= v.y;
        }
    
        /** ditto */
        Vector2 opMul(real k)
        {
            return Vector2(x * k, y * k);
        }
    
        /** ditto */
        void opMulAssign(real k)
        {
            x *= k;
            y *= k;
        }
    
        /** ditto */
        Vector2 opMul_r(real k)
        {
            return Vector2(x * k, y * k);
        }
    
        /** ditto */
        Vector2 opDiv(real k)
        {
            return Vector2(x / k, y / k);
        }
    
        /** ditto */
        void opDivAssign(real k)
        {
            x /= k;
            y /= k;
        }
    
        /** Returns: Copy of this vector with float type components */
        Vector2f toVector2f()
        {
            return Vector2f(cast(float)x, cast(float)y);
        }
        
        /** Returns: Copy of this vector with double type components */
        Vector2d toVector2d()
        {
            return Vector2d(cast(double)x, cast(double)y);
        }
        
        /** Returns: Copy of this vector with real type components */
        Vector2r toVector2r()
        {
            return Vector2r(cast(real)x, cast(real)y);
        }
    
        /**
        Routines known as swizzling.
        Returns:
            New vector constructed from this one and having component values
            that correspond to method name.
        */
        Vector3 xy0()    { return Vector3(x, y, 0); }
        Vector3 x0y()    { return Vector3(x, 0, y); } /// ditto
    }
    
    /** Returns: Dot product between passed vectors. */
    real dp(Vector2 a, Vector2 b)
    {
        return a.x * b.x + a.y * b.y;
    }
        
    alias EqualityByNorm!(Vector2).equal equal; /// Introduces approximate equality function for Vector2.
    alias Lerp!(Vector2).lerp lerp;             /// Introduces linear interpolaton function for Vector2.
    
    
    /************************************************************************************
    Three dimensional vector.

    For formal definition of vector, meaning of possible operations and related
    information see $(LINK http://en.wikipedia.org/wiki/Vector_(spatial)).
    *************************************************************************************/
    struct Vector3
    {
        align(1)
        {
            float_t x; /// Components of vector.
            float_t y; /// ditto
            float_t z; /// ditto
        }
        
        static Vector3 nan = { float_t.nan, float_t.nan, float_t.nan }; /// Vector with all components seted to NaN.
        static Vector3 unitX = { 1, 0, 0 };  /// Unit vector codirectional with corresponding axis.
        static Vector3 unitY = { 0, 1, 0 };  /// ditto
        static Vector3 unitZ = { 0, 0, 1 };  /// ditto
        
        /**
        Method to construct vector in C-like syntax.

        Examples:
        ------------
        Vector3 myVector = Vector3(1, 2, 3);
        ------------
        */
        static Vector3 opCall(float_t x, float_t y, float_t z)
        {
            Vector3 v;
            v.set(x, y, z);
            return v;
        }
        
        /** Sets values of components to passed values. */
        void set(float_t x, float_t y, float_t z)
        {
            this.x = x;
            this.y = y;
            this.z = z;
        }
    
        /** Returns: Norm (also known as length, magnitude) of vector. */
        real norm()
        {
            return sqrt(x*x + y*y + z*z);
        }
    
        /**
        Returns: Square of vector's norm.
        
        Since this method doesn't need calculation of square root it is better
        to use it instead of norm() when you can. For example, if you want just
        to know which of 2 vectors is longer it's better to compare their norm
        squares instead of their norm.
        */
        real normSquare()
        {
            return x*x + y*y + z*z;
        }
    
        /** Normalizes this vector. */
        void normalize()
        {
            *this /= norm;
        }
    
        /** Returns: Normalized copy of this vector. */
        Vector3 normalized()
        {
            real n = norm;
            return Vector3(x / n, y / n, z / n);
        }
    
        /**
        Returns: Whether this vector is unit.
        Params:
            relprec, absprec = Parameters passed to equal function while comparison of
                               norm square and 1. Have the same meaning as in equal function.
        */
        bool isUnit(int relprec = defrelprec, int absprec = defabsprec)
        {
            return equal( normSquare(), 1, relprec, absprec );
        }
    
        /** Returns: Axis for which projection of this vector on it will be longest. */
        Ort dominatingAxis()
        {
            if (x > y)
                return (x > z) ? Ort.X : Ort.Z;
            else
                return (y > z) ? Ort.Y : Ort.Z;
        }
    
        /** Returns: Whether all components are normalized numbers. */
        bool isnormal()
        {
            return std.math.isnormal(x) && std.math.isnormal(y) && std.math.isnormal(z);
        }
    
        /** Returns: float_t pointer to x component of this vector. It's like a _ptr method for arrays. */
        float_t* ptr()
        {
            return cast(float_t*)this;
        }
    
        /** Returns: Component corresponded to passed index. */
        float_t opIndex(Ort ort)
        in { assert(ort <= Ort.Z); }
        body
        {
            return ptr[cast(int)ort];
        }
    
        /** Assigns new _value to component corresponded to passed index. */
        void opIndexAssign(float_t value, Ort ort)
        in { assert(ort <= Ort.Z); }
        body
        {
            ptr[cast(int)ort] = value;
        }
    
        /**
        Standard operators that have intuitive meaning, same as in classical math.
        
        Note that division operators do no cheks of value of k, so in case of division
        by 0 result vector will have infinity components. You can check this with isnormal()
        method.
        */
        bool opEquals(Vector3 v)
        {
            return x == v.x && y == v.y && z == v.z;
        }
    
        /** ditto */
        Vector3 opNeg()
        {
            return Vector3(-x, -y, -z);
        }
    
        /** ditto */
        Vector3 opAdd(Vector3 v)
        {
            return Vector3(x + v.x, y + v.y, z + v.z);
        }
    
        /** ditto */
        void opAddAssign(Vector3 v)
        {
            x += v.x;
            y += v.y;
            z += v.z;
        }
    
        /** ditto */
        Vector3 opSub(Vector3 v)
        {
            return Vector3(x - v.x, y - v.y, z - v.z);
        }
    
        /** ditto */
        void opSubAssign(Vector3 v)
        {
            x -= v.x;
            y -= v.y;
            z -= v.z;
        }
    
        /** ditto */
        Vector3 opMul(real k)
        {
            return Vector3(x * k, y * k, z * k);
        }
    
        /** ditto */
        void opMulAssign(real k)
        {
            x *= k;
            y *= k;
            z *= k;
        }
    
        /** ditto */
        Vector3 opMul_r(real k)
        {
            return Vector3(x * k, y * k, z * k);
        }
    
        /** ditto */
        Vector3 opDiv(real k)
        {
            return Vector3(x / k, y / k, z / k);
        }
    
        /** ditto */
        void opDivAssign(real k)
        {
            x /= k;
            y /= k;
            z /= k;
        }
        
        /** Returns: Copy of this vector with float type components */
        Vector3f toVector3f()
        {
            return Vector3f(cast(float)x, cast(float)y, cast(float)z);
        }
        
        /** Returns: Copy of this vector with double type components */
        Vector3d toVector3d()
        {
            return Vector3d(cast(double)x, cast(double)y, cast(double)z);
        }

        /** Returns: Copy of this vector with real type components */
        Vector3r toVector3r()
        {
            return Vector3r(cast(real)x, cast(real)y, cast(real)z);
        }

    
        /**
        Routines known as swizzling.
        Returns:
            New vector constructed from this one and having component values
            that correspond to method name.
        */
        Vector4 xyz0()        { return Vector4(x,y,z,0); }
        Vector4 xyz1()        { return Vector4(x,y,z,1); } /// ditto
        Vector2 xy()          { return Vector2(x, y); }    /// ditto
        Vector2 xz()          { return Vector2(x, z); }    /// ditto
        
        /**
        Routines known as swizzling.
        Assigns new values to some components corresponding to method name.
        */
        void xy(Vector2 v)    { x = v.x; y = v.y; }
        void xz(Vector2 v)    { x = v.x; z = v.y; }        /// ditto
    }
    
    /** Returns: Dot product between passed vectors. */
    real dp(Vector3 a, Vector3 b)
    {
        return a.x * b.x + a.y * b.y + a.z * b.z;
    }
    
    /**
    Returns: Cross product between passed vectors. Result is vector c
    so that a, b, c forms right-hand tripple.
    */
    Vector3 cp(Vector3 a, Vector3 b)
    {
        return Vector3(
            a.y * b.z - b.y * a.z,
            a.z * b.x - b.z * a.x,
            a.x * b.y - b.x * a.y  );
    }
    
    /**
    Returns: Whether passed basis is orthogonal.
    Params:
        r, s, t =          Vectors that form a basis.
        relprec, absprec = Parameters passed to equal function while calculations.
                           Have the same meaning as in equal function.
    References:
        $(LINK http://en.wikipedia.org/wiki/Orthogonal_basis)
    */
    bool isBasisOrthogonal(Vector3 r, Vector3 s, Vector3 t, int relprec = defrelprec, int absprec = defabsprec)
    {
        return equal( cp(r, cp(s, t)).normSquare, 0, relprec, absprec );
    }
    
    /**
    Returns: Whether passed basis is orthonormal.
    Params:
        r, s, t =   Vectors that form a basis.
        relprec, absprec = Parameters passed to equal function while calculations.
                           Have the same meaning as in equal function.
    References:
        $(LINK http://en.wikipedia.org/wiki/Orthonormal_basis)
    */
    bool isBasisOrthonormal(Vector3 r, Vector3 s, Vector3 t, int relprec = defrelprec, int absprec = defabsprec)
    {
        return isBasisOrthogonal(r, s, t, relprec, absprec) &&
            r.isUnit(relprec, absprec) &&
            s.isUnit(relprec, absprec) &&
            t.isUnit(relprec, absprec);
    }
    
    alias EqualityByNorm!(Vector3).equal equal; /// Introduces approximate equality function for Vector3.
    alias Lerp!(Vector3).lerp lerp;             /// Introduces linear interpolation function for Vector3.
    
    
    /************************************************************************************
    4D vector.

    For formal definition of vector, meaning of possible operations and related
    information see $(LINK http://en.wikipedia.org/wiki/Vector_(spatial)),
    $(LINK http://en.wikipedia.org/wiki/Homogeneous_coordinates).
    *************************************************************************************/
    struct Vector4
    {
        align(1)
        {
            float_t x; /// Components of vector.
            float_t y; /// ditto
            float_t z; /// ditto
            float_t w; /// ditto
        }
        
        /// Vector with all components seted to NaN.
        static Vector4 nan = { float_t.nan, float_t.nan, float_t.nan, float_t.nan };
        static