// Magic Software, Inc.
// http://www.magic-software.com
// http://www.wild-magic.com
// Copyright (c) 2003.  All Rights Reserved
//
// The Wild Magic Library (WML) source code is supplied under the terms of
// the license agreement http://www.magic-software.com/License/WildMagic.pdf
// and may not be copied or disclosed except in accordance with the terms of
// that agreement.

#ifndef WMLMATRIX3_H
#define WMLMATRIX3_H

// The (x,y,z) coordinate system is assumed to be right-handed.  Coordinate
// axis rotation matrices are of the form
//   RX =    1       0       0
//           0     cos(t) -sin(t)
//           0     sin(t)  cos(t)
// where t > 0 indicates a counterclockwise rotation in the yz-plane
//   RY =  cos(t)    0     sin(t)
//           0       1       0
//        -sin(t)    0     cos(t)
// where t > 0 indicates a counterclockwise rotation in the zx-plane
//   RZ =  cos(t) -sin(t)    0
//         sin(t)  cos(t)    0
//           0       0       1
// where t > 0 indicates a counterclockwise rotation in the xy-plane.

#include "WmlMatrix.h"
#include "WmlVector3.h"
#include "WmlMatrix4.h"

namespace Wml
{
  template<class Real>
  class WML_ITEM Matrix3: public Matrix<3, Real>
  {
  public:
    // construction
                                Matrix3();
                                Matrix3(const Matrix3& rkM);
                                Matrix3(const Matrix<3, Real>& rkM);

    // input Mrc is in row r, column c.
                                Matrix3(Real fM00,
                                        Real fM01,
                                        Real fM02,
                                        Real fM10,
                                        Real fM11,
                                        Real fM12,
                                        Real fM20,
                                        Real fM21,
                                        Real fM22);

                                operator Matrix4<Real>() const;

    // Create a matrix from an array of numbers.  The input array is
    // interpreted based on the Boolean input as
    //   true:  entry[0..8]={m00,m01,m02,m10,m11,m12,m20,m21,m22} [row major]
    //   false: entry[0..8]={m00,m10,m20,m01,m11,m21,m02,m12,m22} [col major]
                                Matrix3(const Real afEntry[9], bool bRowMajor);

    // Create matrices based on vector input.  The Boolean is interpreted as
    //   true: vectors are columns of the matrix
    //   false: vectors are rows of the matrix
                                Matrix3(const Vector3<Real>& rkU,
                                        const Vector3<Real>& rkV,
                                        const Vector3<Real>& rkW,
                                        bool bColumns);
                                Matrix3(const Vector3<Real>* akV, bool bColumns);

    // create a tensor product U*V^T
                                Matrix3(const Vector3<Real>& rkU, const Vector3<Real>& rkV);
    void                        MakeTensorProduct(const Vector3<Real>& rkU, const Vector3<Real>& rkV);

    // create a diagonal matrix
                                Matrix3(Real fM00, Real fM11, Real fM22);
    void                        MakeDiagonal(Real fM00, Real fM11, Real fM22);

    // Create rotation matrices (positive angle - counterclockwise).  The
    // angle must be in radians, not degrees.
                                Matrix3(const Vector3<Real>& rkAxis, Real fAngle);
    void                        FromAxisAngle(const Vector3<Real>& rkAxis, Real fAngle);

    // assignment
    Matrix3&                    operator=(const Matrix3& rkM);
    Matrix3&                    operator=(const Matrix<3, Real>& rkM);

    // matrix operations
    Matrix3                     Inverse() const;
    Matrix3                     Adjoint() const;
    Real                        Determinant() const;

    // The matrix must be a rotation for these functions to be valid.  The
    // last function uses Gram-Schmidt orthonormalization applied to the
    // columns of the rotation matrix.  The angle must be in radians, not
    // degrees.
    void                        ToAxisAngle(Vector3<Real>& rkAxis, Real& rfAngle) const;
    void                        Orthonormalize();

    // The matrix must be orthonormal.  The decomposition is yaw*pitch*roll
    // where yaw is rotation about the Up vector, pitch is rotation about the
    // Right axis, and roll is rotation about the Direction axis.
    bool                        ToEulerAnglesXYZ(Real& rfYAngle, Real& rfPAngle, Real& rfRAngle) const;
    bool                        ToEulerAnglesXZY(Real& rfYAngle, Real& rfPAngle, Real& rfRAngle) const;
    bool                        ToEulerAnglesYXZ(Real& rfYAngle, Real& rfPAngle, Real& rfRAngle) const;
    bool                        ToEulerAnglesYZX(Real& rfYAngle, Real& rfPAngle, Real& rfRAngle) const;
    bool                        ToEulerAnglesZXY(Real& rfYAngle, Real& rfPAngle, Real& rfRAngle) const;
    bool                        ToEulerAnglesZYX(Real& rfYAngle, Real& rfPAngle, Real& rfRAngle) const;
    void                        FromEulerAnglesXYZ(Real fYAngle, Real fPAngle, Real fRAngle);
    void                        FromEulerAnglesXZY(Real fYAngle, Real fPAngle, Real fRAngle);
    void                        FromEulerAnglesYXZ(Real fYAngle, Real fPAngle, Real fRAngle);
    void                        FromEulerAnglesYZX(Real fYAngle, Real fPAngle, Real fRAngle);
    void                        FromEulerAnglesZXY(Real fYAngle, Real fPAngle, Real fRAngle);
    void                        FromEulerAnglesZYX(Real fYAngle, Real fPAngle, Real fRAngle);

    // SLERP (spherical linear interpolation) without quaternions.  Computes
    // R(t) = R0*(Transpose(R0)*R1)^t.  If Q is a rotation matrix with
    // unit-length axis U and angle A, then Q^t is a rotation matrix with
    // unit-length axis U and rotation angle t*A.
    static Matrix3              Slerp(Real fT, const Matrix3& rkR0, const Matrix3& rkR1);

    // The matrix must be symmetric.  Factor M = R * D * R^T where
    // R = [u0|u1|u2] is a rotation matrix with columns u0, u1, and u2 and
    // D = diag(d0,d1,d2) is a diagonal matrix whose diagonal entries are d0,
    // d1, and d2.  The eigenvector u[i] corresponds to eigenvector d[i].
    // The eigenvalues are ordered as d0 <= d1 <= d2.
    void                        EigenDecomposition(Matrix3& rkRot, Matrix3& rkDiag) const;

    // Singular value decomposition, M = L*S*R, where L and R are orthogonal
    // and S is a diagonal matrix whose diagonal entries are nonnegative.
    void                        SingularValueDecomposition(Matrix3& rkL, Matrix3& rkS, Matrix3& rkR) const;
    void                        SingularValueComposition(const Matrix3& rkL, const Matrix3& rkS, const Matrix3& rkR);

    // factor M = Q*D*U with orthogonal Q, diagonal D, upper triangular U
    void                        QDUDecomposition(Matrix3& rkQ, Matrix3& rkD, Matrix3& rkU) const;

    void                        MakeRotationX(Real angle);
    void                        MakeRotationY(Real angle);
    void                        MakeRotationZ(Real angle);

    // special matrices
    static const Matrix3        ZERO;
    static const Matrix3        IDENTITY;

  protected:
    // support for eigendecomposition
    void                        Tridiagonalize(Real afDiag[3], Real afSubDiag[3]);
    bool                        QLAlgorithm(Real afDiag[3], Real afSubDiag[3]);

    // support for singular value decomposition
    static void                 Bidiagonalize(Matrix3& rkA, Matrix3& rkL, Matrix3& rkR);
    static void                 GolubKahanStep(Matrix3& rkA, Matrix3& rkL, Matrix3& rkR);
  };

  typedef Matrix3<float> Matrix3f;
  typedef Matrix3<double> Matrix3d;
}

#endif
