#ifndef __MATRIX4_H__
#define __MATRIX4_H__

#include <iostream>
#include <cassert>
#include <cmath>

#include "vector3.h"
#include "vector4.h"
#include "bbox3.h"

#include "public.h"


namespace xrt {

class UnMatrix4;

// The Matrix4 class is for a 4x4 matrix.
//
// We GUARANTEE that the memory layout and size of a Matrix4 is
// identical to a float[4][4], and thus it's presumed relatively safe
// to cast between a (float *) and a (Matrix4 *).  We also presume
// that the components of matrices (in either Matrix4 or float[]) has
// the translation portion in elements 12,13,14 -- this is as expected
// by both OpenGL and RenderMan.
//
// There are two Matrix4 constants, given by Matrix4::Zero() and
// Matrix4::Ident().  It may be handy to use these by name to
// provide fastest copying of the values.
//
// Component access is by [int][int].
//
// Multiplication with floats or matrices works as expected.
// Routines for transpose and inverse, as well as point/vector
// (implied 1 and 0 as 4th component, respectively) transformation
// are supplied.
//

class PUBLIC Matrix4 {

public:

    // Special Matrix4's -- Matrix4::Zero(), Matrix4::Identity()
    static const Matrix4& Zero() { 
        static const Matrix4 _zero (0.0);
	return _zero;
    }
    static const Matrix4& Identity() { 
        static const Matrix4 _ident (1.0);
	return _ident;
    }

    //////////////
    // Constructors

    // Constructor with no args - leave uninitialized
    Matrix4 (void) { }

    // Construct from a float f makes f*Ident
    Matrix4 (float x) {
	m[ 0] = x;     m[ 1] = 0.0f;  m[ 2] = 0.0f;  m[ 3] = 0.0f;
	m[ 4] = 0.0f;  m[ 5] = x;     m[ 6] = 0.0f;  m[ 7] = 0.0f;
	m[ 8] = 0.0f;  m[ 9] = 0.0f;  m[10] = x;     m[11] = 0.0f;
	m[12] = 0.0f;  m[13] = 0.0f;  m[14] = 0.0f;  m[15] = x;
    }

    // Construct directly from 16 floats
    Matrix4 (float m0, float m1, float m2, float m3,
	     float m4, float m5, float m6, float m7,
	     float m8, float m9, float m10, float m11,
	     float m12, float m13, float m14, float m15)
    {
	m[ 0] = m0;  m[ 1] = m1;  m[ 2] = m2;  m[ 3] = m3;
	m[ 4] = m4;  m[ 5] = m5;  m[ 6] = m6;  m[ 7] = m7;
	m[ 8] = m8;  m[ 9] = m9;  m[10] = m10; m[11] = m11;
	m[12] = m12; m[13] = m13; m[14] = m14; m[15] = m15;
    }

    // Construct from float* presumes 16 contigious floats
    Matrix4 (const float *xarray) {
	memcpy (m, xarray, 16*sizeof(float));
    }

    // Simple assignment
    const Matrix4& operator= (float x) {
	m[ 0] = x;     m[ 1] = 0.0f;  m[ 2] = 0.0f;  m[ 3] = 0.0f;
	m[ 4] = 0.0f;  m[ 5] = x;     m[ 6] = 0.0f;  m[ 7] = 0.0f;
	m[ 8] = 0.0f;  m[ 9] = 0.0f;  m[10] = x;     m[11] = 0.0f;
	m[12] = 0.0f;  m[13] = 0.0f;  m[14] = 0.0f;  m[15] = x;
	return *this;
    }
    // Use the default assignment operator (bitwise copy)

    // Direct access to the data -- use with caution
    const float *data(void) const { return m; }


    // Component access -- [] addresses a row as a float*, 
    // so [][] addresses a float
    const float* operator[] (int i) const {
	assert(i>=0 && i < 4);  // range check -- only in DEBUG mode
	return m+i*4;
    }
    float* operator[] (int i) {
	assert(i>=0 && i < 4);  // range check -- only in DEBUG mode
	return m+i*4;
    }

    // Comparisons between matrices
    friend bool operator== (const Matrix4 &a, const Matrix4 &b) {
	return memcmp (&a, &b, sizeof(Matrix4)) == 0;
    }
    friend bool operator!= (const Matrix4 &a, const Matrix4 &b) {
	return memcmp (&a, &b, sizeof(Matrix4)) != 0;
    }

    // Stream output -- very handy for debugging
    friend std::ostream& operator<< (std::ostream& out, const Matrix4& a) {
	out << '[';
	for (int i = 0;  i < 15;  ++i)
	    out << a.m[i] << ' ';
	out << a.m[15] << ']';
	return out;
    }

    // Addition
    friend Matrix4 operator+ (const Matrix4& a, const Matrix4& b) {
        Matrix4 r;
        for (int i = 0;  i < 16;  ++i)
            r.m[i] = a.m[i] + b.m[i];
        return r;
    }

    const Matrix4& operator+= (const Matrix4& b) {
	for (int i = 0;  i < 16;  ++i)
            m[i] += b.m[i];
	return *this;
    }

    // Subtraction
    friend Matrix4 operator- (const Matrix4& a, const Matrix4& b) {
        Matrix4 r;
        for (int i = 0;  i < 4;  ++i)
            r.m[i] = a.m[i] - b.m[i];
        return r;
    }
    const Matrix4& operator-= (const Matrix4& b) {
	for (int i = 0;  i < 16;  ++i)
            m[i] -= b.m[i];
	return *this;
    }

	// Matrix multiplication
    friend Matrix4 operator* (const Matrix4& a, const Matrix4& b) {
	Matrix4 result;
	for (int j = 0;  j < 16;  ++j) {
	    int col = (j & 0x03);
	    int row = (j & 0x0c);
#ifdef VECMAT_LOOPS
            result.m[j] = 0.0f;
            for (int i = 0;  i < 4;  ++i)
                result.m[j] += a.m[row+i] * b.m[col+4*i];
#else
	    result.m[j] = a.m[row] * b.m[col]
		        + a.m[row+1] * b.m[col+4]
		        + a.m[row+2] * b.m[col+8]
		        + a.m[row+3] * b.m[col+12];
#endif
	}
	return result;
    }

    // Matrix multiplication by a scalar
    friend Matrix4 operator* (const Matrix4& M, float f) {
	Matrix4 result;
	for (int i = 0;  i < 16;  ++i)
	    result.m[i] = M.m[i] * f;
	return result;
    }
    friend Matrix4 operator* (float f, const Matrix4& M) { return M*f; }
    const Matrix4& operator*= (float f) {
	for (int i = 0;  i < 16;  ++i)
	    m[i] *= f;
	return *this;
    }

    // Matrix division by a scalar -- up to the caller to beware of 
    // division-by-zero
    friend Matrix4 operator/ (const Matrix4& M, float f) {
	Matrix4 result;
	for (int i = 0;  i < 16;  ++i)
	    result.m[i] = M.m[i] / f;
	return result;
    }
    const Matrix4& operator/= (float f) {
	for (int i = 0;  i < 16;  ++i)
	    m[i] /= f;
	return *this;
    }

    // float/matrix is the same as float*inverse(matrix)
    friend Matrix4 operator/ (float f, const Matrix4& M) {
	return f * M.inverse();
    }

    // matrix1/matrix2 is the same as matrix1*inverse(matrix2)
    friend Matrix4 operator/ (const Matrix4 &M1, const Matrix4& M2) {
	return M1 * M2.inverse();
    }

	EXPORT friend Matrix4 slerp (float t, const Matrix4& a, const Matrix4& b); 

		// Transpose the rows & columns
    Matrix4 transpose (void) const { 
	return Matrix4 (m[0], m[4], m[ 8], m[12],
			m[1], m[5], m[ 9], m[13],
			m[2], m[6], m[10], m[14],
			m[3], m[7], m[11], m[15]);
    }

    // Matrix inverse
    Matrix4 inverse (float *determinant=NULL) const;
    // Matrix determinant
    float determinant (void) const;

    // Vector3 the Z axis of the current matrix toward the lookat point
    void lookat (const Vector3& lookat, const Vector3& up);
    
    ///////////////////
    // Generate commonly-used transformation matrices
    static Matrix4 Translation (const Vector3& trans);
    static Matrix4 Rotation (float angle, const Vector3& axis);
    static Matrix4 Scale (float sx, float sy, float sz);
    static Matrix4 Lookat (const Vector3& camera,
        const Vector3& lookat, const Vector3& up);
    static Matrix4 Perspective (float fov, float aspect,
        float hither, float yon);
    static Matrix4 Ortho (float left, float right, float bottom,
        float top, float hither, float yon);
    static Matrix4 Frustum (float left, float right, float bottom, 
        float top, float hither, float yon);
    
    // transformp transforms a hpoint 
    // by a matrix.
    Vector4
    transformp (const Vector4 &P) const {
	Vector4 result;
#ifdef VECMAT_LOOPS
        assert (&P != &result);
        for (int i = 0;  i < 4;  ++i)
            result[i] = P[0] * m[4*i+0] + P[1] * m[4*i+1] + P[2] * m[4*i+2] + P[3] * m[4*i+3];
#else
	register float x, y, z, w;
	x = P[0] * m[0] + P[1] * m[1] + P[2] * m[2] + P[3] * m[3];
	y = P[0] * m[4] + P[1] * m[5] + P[2] * m[6] + P[3] * m[7];
	z = P[0] * m[8] + P[1] * m[9] + P[2] * m[10] + P[3] * m[11];
	w = P[0] * m[12] + P[1] * m[13] + P[2] * m[14] + P[3] * m[15];
	result[0] = x;
	result[1] = y;
	result[2] = z;
	result[3] = w;
#endif
	return result;
    }

     // transformp transforms a point (with assumed 4th component == 1)
    // by a matrix.
    Vector3
    transformp (const Vector3 &P) const {
	Vector3 result;
#ifdef VECMAT_LOOPS
        float r[4];
#if 1
        for (int i = 0;  i < 4;  ++i)
            result[i] = P[0] * m[4*i+0] + P[1] * m[4*i+1] + P[2] * m[4*i+2] + m[4*i+3];
#else  /* alternate - may be better on certain architectures */
        for (int i = 0;  i < 4;  ++i)
            r[i] = m[4*i+3];
        for (int j = 0; j < 3;  ++j)
            for (int i = 0;  i < 4;  ++i)
                r[i] += P[j] * m[4*i+j];
#endif

	if (r[3] != 1.0f) {
	    float w = 1.0f/r[3];
            for (int i = 0;  i < 3;  ++i)
                result[i] = r[i] * w;
        } else {
            for (int i = 0;  i < 3;  ++i)
                result[i] = r[i];
        }
#else
	register float x, y, z, w;
	x = P[0] * m[0] + P[1] * m[1] + P[2] * m[2] + m[3];
	y = P[0] * m[4] + P[1] * m[5] + P[2] * m[6] + m[7];
	z = P[0] * m[8] + P[1] * m[9] + P[2] * m[10] + m[11];
	w = P[0] * m[12] + P[1] * m[13] + P[2] * m[14] + m[15];
	if (w != 1.0f) {
	    w = 1.0f/w;
	    x *= w;
	    y *= w;
	    z *= w;
	}
	result[0] = x;
	result[1] = y;
	result[2] = z;
#endif
	return result;
	}

	// transformv transforms a vector (with assumed 4th component == 0)
    // by a matrix.
    Vector3
    transformv (const Vector3 &P) const {
	Vector3 result;
#ifdef VECMAT_LOOPS
	float r[4];
#if 1
        for (int i = 0;  i < 3;  ++i)
            result[i] = P[0] * m[4*i+0] + P[1] * m[4*i+1] + P[2] * m[4*i+2];
#else  /* alternate - may be better on some architectures */
        for (int i = 0;  i < 4;  ++i)
            r[i] = 0.0f;
        for (int j = 0; j < 3;  ++j)
            for (int i = 0;  i < 3;  ++i)
                r[i] += P[j] * m[4*i+j];
#endif
        for (int i = 0;  i < 3;  ++i)
            result[i] = r[i];
#else
	register float x, y, z;
	x = P[0] * m[0] + P[1] * m[1] + P[2] * m[2];
	y = P[0] * m[4] + P[1] * m[5] + P[2] * m[6];
	z = P[0] * m[8] + P[1] * m[9] + P[2] * m[10];
	result[0] = x;
	result[1] = y;
	result[2] = z;
#endif
	return result;
	}

	// transformvT transforms a vector (with assumed 4th component == 0)
    // by the transpose of a matrix (handy for transforming normals).
    Vector3
    transformvT (const Vector3 &P) const {
	Vector3 result;
#ifdef VECMAT_LOOPS
        float r[4];
        for (int i = 0;  i < 3;  ++i)
            r[i] = P[0] * m[i+0] + P[1] * m[i+4] + P[2] * m[i+8];
        for (int i = 0;  i < 3;  ++i)
            result[i] = r[i];
#else
	register float x, y, z;
	x = P[0] * m[0] + P[1] * m[4] + P[2] * m[8];
	y = P[0] * m[1] + P[1] * m[5] + P[2] * m[9];
	z = P[0] * m[2] + P[1] * m[6] + P[2] * m[10];
	result[0] = x;
	result[1] = y;
	result[2] = z;
#endif
	return result;
	}

	 // Transform normal -- must use the inverse transpose
    Vector3
    transformn (const Vector3 &P) const {
        Matrix4 Minv = inverse();
        return Minv.transformvT (P);
    }

    // Transform an array of points into a new array.
    void transformp (const Vector3 *P, Vector3 *R,
                            int n) const;
    void transformv (const Vector3 *P, Vector3 *R,
                            int n) const;
    void transformvT (const Vector3 *P, Vector3 *R,
                            int n) const;
    void transformn (const Vector3 *P, Vector3 *R,
                            int n) const;
    void transformp (const Vector4 *P, Vector4 *R,
                            int n) const;

	// Transform an array of points or vectors in place
	void transformp (Vector3 *P, int n) const { transformp (P, P, n); }
    void transformv (Vector3 *P, int n) const { transformv (P, P, n); }
    void transformvT (Vector3 *P, int n) const { transformvT (P, P, n); }
	void transformn (Vector3 *P, int n) const { transformn (P, P, n); }
    void transformp (Vector4 *P, int n) const { transformp (P, P, n); }

    // Return just the z component of point P transformed by M.  If you
    // only need z, this is cheaper than the whole matrix multiply.
    float
    transformp_z (const Vector3 &P) const {
	register float z, w;
	z = P[0] * m[8] + P[1] * m[9] + P[2] * m[10] + m[11];
	w = P[0] * m[12] + P[1] * m[13] + P[2] * m[14] + m[15];
	if (w != 1.0f)
	    z /= w;
        return z;
    }

    // Return just the z component of vector P transformed by M.  If you
    // only need z, this is cheaper than the whole matrix multiply.
    float
    transformv_z (const Vector3 &P) const {
	return P[0] * m[8] + P[1] * m[9] + P[2] * m[10];
    }

    // Return just the xy components of point P transformed by M.  If you
    // only need x and y, this is cheaper than the whole matrix multiply.
    // CAVEAT: result[2] is not modified.
    friend void
    transformp_xy (const Matrix4& M, const Vector3 &P, Vector3 &result) {
	register float x, y, w;
	x = P[0] * M.m[0] + P[1] * M.m[1] + P[2] * M.m[2] + M.m[3];
	y = P[0] * M.m[4] + P[1] * M.m[5] + P[2] * M.m[6] + M.m[7];
	w = P[0] * M.m[12] + P[1] * M.m[13] + P[2] * M.m[14] + M.m[15];
	if (w != 1.0f) {
	    w = 1.0f/w;
	    x *= w;
	    y *= w;
	}
	result[0] = x;
	result[1] = y;
    }

    friend float transform_length (const Matrix4& M, const float length) {
        // transform the length as three separate vectors, (length, 0, 0),
        // (0, length, 0) and (0, 0, length) and compute the average length
        // of the three transformed vectors.  Note that v*M for each of
        // the above vectors corresponds to a row in the matrix.
        float d = ( M.m[0] * M.m[0] + M.m[1] * M.m[1] + M.m[2] * M.m[2] +
                    M.m[4] * M.m[4] + M.m[5] * M.m[5] + M.m[6] * M.m[6] +
                    M.m[8] * M.m[8] + M.m[9] * M.m[9] + M.m[10] * M.m[10] ) /3;
        return length * sqrt (d);
    }

    // btransform transforms a bbox by a matrix.
    BBox3 btransform (const BBox3 &b) const;
 
	// transform and compute only the z range of the bbox
    void btransform_z (const BBox3 &b, float &zmin, float &zmax)  const {
        zmin = zmax = transformp_z (b.corner (0));
        for (int i = 1 /*ONE*/;  i < 8;  ++i) {
            float z = transformp_z (b.corner (i));
            // parens around function names so min/max macros aren't a problem
            zmin = (std::min) (zmin, z);
            zmax = (std::max) (zmax, z);
        }
    }

	bool isLeftHanded (void) const;

	bool hasIsotropicScale (float& scaleFactor) const;

	UnMatrix4 decompose() const;

private:
    float m[16];

};


};  /* end namespace xrt */


#endif  /* !defined(__MATRIX4_H__) */