public class Matrix3D
{
    float[][] matrix;
    
    //        
    // Constructors        
    //
    
    public Matrix3D()                  // initialize with identity transform
    {
        matrix = new float[4][4];
        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                if (i == j)
                    matrix[j][i] = 1.0f;
                else
                    matrix[j][i] = 0.0f;
            }
        }
    }
    public Matrix3D(Matrix3D copy)     // initialize with copy of source
    {
        matrix = new float[4][4];
        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                matrix[j][i] = copy.matrix[j][i];
            }
        }
    }

                
    public Matrix3D(Raster r)          // initialize with a mapping from
    {                                  // canonical space to screen space
        this();
        
        float zmax = 200;
        
        translate(r.width/2, r.height/2, 0);       
        scale(r.width/2, r.height/2, zmax/2);
    }

    //        
    // General interface methods        
    //   
    
    public void set(int i, int j, float value)          // set element [j][i] to value
    {
        matrix[j][i] = value;
    }
    
    public float get(int i, int j)                      // return element [j][i]
    {
        return matrix[j][i];
    }
    
    //        
    // Transform points from the in array to the out array
    // using the current matrix. The subset of points transformed
    // begins at the start index and has the specified length        
    //
    //    for (i = 0; i < length; i++)
    //        out[start+i] = this * in[start+i]        
    //
    public void transform(Point3D in[], Point3D out[], int start, int length)
    {
        for (int i = 0; i < length; i++)
        {
            out[start + i] = ptimes(in[start + i]);
        }
    }

    public final void compose(Matrix3D src)                        // this = this * src
    {
        Matrix3D copy = new Matrix3D(this);
        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                set(i,j,copy.get(i,0)*src.get(0,j)+
                        copy.get(i,1)*src.get(1,j)+
                        copy.get(i,2)*src.get(2,j)+
                        copy.get(i,3)*src.get(3,j));
            }
        }
    }
    
    public void printout()
    {
        for (int i = 0; i < 4; i++)
        {
            System.out.println(this.matrix[0][i]+"   "+this.matrix[1][i]+"   "+this.matrix[2][i]+"   "+this.matrix[3][i]);
        }
        System.out.println(" ");
    }
    
    public void loadIdentity()                                      // this = identity
    {
        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                if (i == j)
                    matrix[j][i] = 1.0f;
                else
                    matrix[j][i] = 0.0f;
            }
        }
    }
    
    public void loadSymmetric(float ax, float ay, float az)         // this = symmetric
    {
        set(0, 0, ax*ax);  set(0, 1, ax*ay);  set(0, 2, ax*az);  set(0, 3, 0);    
        set(1, 0, ay*ax);  set(1, 1, ay*ay);  set(1, 2, ay*az);  set(1, 3, 0);    
        set(2, 0, az*ax);  set(2, 1, az*ay);  set(2, 2, az*az);  set(2, 3, 0);    
        set(3, 0, 0);      set(3, 1, 0);      set(3, 2, 0);      set(3, 3, 1);    
    }
    
    public void loadSkewSym(float ax, float ay, float az)           // this = skewSymmetric
    {
        set(0, 0, 0);      set(0, 1, -az);    set(0, 2, ay);     set(0, 3, 0);    
        set(1, 0, az);     set(1, 1, 0);      set(1, 2, -ax);    set(1, 3, 0);    
        set(2, 0, -ay);    set(2, 1, ax);     set(2, 2, 0);      set(2, 3, 0);    
        set(3, 0, 0);      set(3, 1, 0);      set(3, 2, 0);      set(3, 3, 1);    
    }
    
    public void stimes(float scalar)                                // this = this * scalar
    {
        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                matrix[j][i] *= scalar;
            }
        }
    }
    
    public Point3D ptimes(Point3D in)
    {
        float x, y, z, w;
        
        x = matrix[0][0] * in.x + matrix[1][0] * in.y + matrix[2][0] * in.z + matrix[3][0];
        y = matrix[0][1] * in.x + matrix[1][1] * in.y + matrix[2][1] * in.z + matrix[3][1];
        z = matrix[0][2] * in.x + matrix[1][2] * in.y + matrix[2][2] * in.z + matrix[3][2];
        w = matrix[0][3] * in.x + matrix[1][3] * in.y + matrix[2][3] * in.z + matrix[3][3];
        return (new Point3D(x/w, y/w, z/w));
    }  
    
    public Point3D normalize(Point3D in)
    {
        float mag = (float)Math.sqrt(in.x * in.x + in.y * in.y + in.z * in.z);
        return new Point3D(in.x / mag, in.y / mag, in.z / mag);
    }
    
    public Point3D cross(Point3D a, Point3D b)
    {
        return new Point3D(-a.z * b.y + a.y * b.z, a.z * b.x - a.x * b.z, -a.y * b.x + a.x * b.y);
    }
    
    public float dot(Point3D a, Point3D b)
    {
        return (a.x * b.x + a.y * b.y + a.z * b.z);
    }
       
    public void mplus(Matrix3D m)                                   // this = this + m
    {
        for (int i = 0; i < 4; i++)
        {
            for (int j = 0; j < 4; j++)
            {
                matrix[j][i] = matrix[j][i] + m.matrix[j][i];
            }
        }
    }
    
    public void translate(float tx, float ty, float tz)             // this = this * t
    {
        Matrix3D t = new Matrix3D();
        t.set(0, 3, tx);
        t.set(1, 3, ty);
        t.set(2, 3, tz);
        compose(t);
    }

    public void scale(float sx, float sy, float sz)                // this = this * scale
    {
        Matrix3D s = new Matrix3D();
        s.set(0, 0, sx);
        s.set(1, 1, sy);
        s.set(2, 2, sz);
        compose(s);
    }
    
    public void skew(float kxy, float kxz, float kyz)              // this = this * skew
    {
        Matrix3D k = new Matrix3D();
        k.set(0, 1, kxy);
        k.set(0, 2, kxz);
        k.set(1, 2, kyz);
        compose(k);
    }
    
    public void rotate(float ax, float ay, float az, float angle)       // this = this * rotate 
    {
        Point3D axis = normalize(new Point3D(ax, ay, az));
        
        Matrix3D sym = new Matrix3D();
        Matrix3D skewsym = new Matrix3D();
        Matrix3D rot = new Matrix3D();
        sym.loadSymmetric(axis.x, axis.y, axis.z);
        skewsym.loadSkewSym(axis.x, axis.y, axis.z);
        rot.stimes((float)Math.cos(angle));
        skewsym.stimes((float)Math.sin(angle));
        sym.stimes((float)(1 - Math.cos(angle)));
        rot.mplus(skewsym);
        rot.mplus(sym);
        rot.set(3,3,1);
        compose(rot);
    }
        
    public void lookAt(float eyex, float eyey, float eyez,
                       float atx,  float aty,  float atz,
                       float upx,  float upy,  float upz)               // this = this * lookat
    {
        Point3D l, r, u, eye;
        Matrix3D lA = new Matrix3D();

        l = normalize(new Point3D(atx - eyex, aty - eyey, atz - eyez));
        r = normalize(cross(l,new Point3D(upx, upy, upz)));
        u = normalize(cross(r, l));
        eye = new Point3D(eyex, eyey, eyez);
        
        lA.set(0, 0, r.x);      lA.set(0, 1, r.y);      lA.set(0, 2, r.z);      lA.set(0, 3, -dot(r, eye));
        lA.set(1, 0, u.x);      lA.set(1, 1, u.y);      lA.set(1, 2, u.z);      lA.set(1, 3, -dot(u, eye));
        lA.set(2, 0, -l.x);     lA.set(2, 1, -l.y);     lA.set(2, 2, -l.z);     lA.set(2, 3, dot(l, eye));
        compose(lA);
    }

    //        
    // Assume the following projection transformations
    // transform points into the canonical viewing space        
    //
    public void perspective(float left, float right,                // this = this * persp
                            float bottom, float top,
                            float near, float far)
    {
        Matrix3D persp = new Matrix3D();
        persp.set(0, 0, (2*near)/(right-left));     persp.set(0, 2, (-right-left)/(right-left));
        persp.set(1, 1, (2*near)/(bottom-top));     persp.set(1, 2, (-bottom-top)/(bottom-top));
        persp.set(2, 2, (far+near)/(far-near));     persp.set(2, 3, (-2*far*near)/(far-near));
        persp.set(3, 2, 1);
        persp.set(3, 3, 0);
        compose(persp);
    }

    
    public void orthographic(float left, float right,               // this = this * ortho
                             float bottom, float top,
                             float near, float far)
    {
        Matrix3D orth = new Matrix3D();
        orth.set(0, 0, 2/(right-left));     orth.set(0, 3, -(right+left)/(right-left));
        orth.set(1, 1, 2/(bottom-top));     orth.set(1, 3, -(bottom+top)/(bottom-top));
        orth.set(2, 2, 2/(far-near));       orth.set(2, 3, -(far+near)/(far-near));
        compose(orth);
    }
}