import java.Math.*;

public class Matrix3D {
    
        public float[][] m4x4 = new float[4][4];
        
        //
        // Constructors
        //
        public Matrix3D() {
            for (int r=0;r<4;r++) {
                for (int c=0;c<4;c++) {
                    if (r==c) {
                        m4x4[r][c]=1;
                    }
                    else {m4x4[r][c]=0;}
                }
            }
        }
        
        public Matrix3D(Matrix3D copy) {
            this.m4x4 = copy.m4x4;
        }
        
        public Matrix3D(Raster r){// initialize with a mapping from canonical space to screen space
            
            for (int row=0;row<4;row++) {
                for (int col=0;col<4;col++) {
                    if (row==col) {
                        m4x4[row][col]=1;
                    }
                    else {m4x4[row][col]=0;}
                }
            }
            
            
            //adjust with scales
            this.m4x4[0][0]=(r.width-1)/2;
            this.m4x4[1][1]=(r.height-1)/2;
            this.m4x4[2][2]=0.5f*(r.Zfar-r.Znear);	
            
            //adjust with translation
            this.m4x4[0][3]=(r.width-1)/2;
            this.m4x4[1][3]=(r.height-1)/2;
            this.m4x4[2][3]=0.5f*(r.Zfar+r.Znear);
        }
        
        //
        // General interface methods
        //
        public void set(int i, int j, float value) { // set element [j][i] to value
            m4x4[j][i]=value;
        }
        
        public float get(int i, int j) { // return element [j][i]
            return (m4x4[j][i]);
        }

        public void transform(Point3D in[], Point3D out[], int start, int length)
        {
            for (int i=0;i<length;i++){
                    out[start+i].x = this.m4x4[0][0]*in[start+i].x+this.m4x4[0][1]*in[start+i].y+this.m4x4[0][2]*in[start+i].z+this.m4x4[0][3];
                    out[start+i].y = this.m4x4[1][0]*in[start+i].x+this.m4x4[1][1]*in[start+i].y+this.m4x4[1][2]*in[start+i].z+this.m4x4[1][3];
                    out[start+i].z = this.m4x4[2][0]*in[start+i].x+this.m4x4[2][1]*in[start+i].y+this.m4x4[2][2]*in[start+i].z+this.m4x4[2][3];
                    
                    //check for problems "at infinity"
                    float w;
                    w = (in[start+i].x*m4x4[3][0]+in[start+i].y*m4x4[3][1]+in[start+i].z*m4x4[3][2]+m4x4[3][3]);
                    if (!(w==1)) {
                        out[start+i].x/=w;
                        out[start+i].y/=w;
                        out[start+i].z/=w;
                    }
            }
            
        }
        
        public Point3D transform(Point3D in)
        {
                    Point3D out = new Point3D();
                    out.x = this.m4x4[0][0]*in.x+this.m4x4[0][1]*in.y+this.m4x4[0][2]*in.z+this.m4x4[0][3];
                    out.y = this.m4x4[1][0]*in.x+this.m4x4[1][1]*in.y+this.m4x4[1][2]*in.z+this.m4x4[1][3];
                    out.z = this.m4x4[2][0]*in.x+this.m4x4[2][1]*in.y+this.m4x4[2][2]*in.z+this.m4x4[2][3];
                    
                    //check for problems "at infinity"
                    float w;
                    w = (in.x*m4x4[3][0]+in.y*m4x4[3][1]+in.z*m4x4[3][2]+m4x4[3][3]);
                    if (!(w==1)) {
                        out.x/=w;
                        out.y/=w;
                        out.z/=w;
                    }
                    return (out);
            
        }
        
        public void transform(Vertex3D in[], Vertex3D out[], int length)
        {
            for (int i=0;i<length;i++){
                    out[i].x = this.m4x4[0][0]*in[i].x+this.m4x4[0][1]*in[i].y+this.m4x4[0][2]*in[i].z+this.m4x4[0][3];
                    out[i].y = this.m4x4[1][0]*in[i].x+this.m4x4[1][1]*in[i].y+this.m4x4[1][2]*in[i].z+this.m4x4[1][3];
                    out[i].z = this.m4x4[2][0]*in[i].x+this.m4x4[2][1]*in[i].y+this.m4x4[2][2]*in[i].z+this.m4x4[2][3];
                    
                    
                    //check for problems "at infinity"
                    float w;
                    w = (in[i].x*m4x4[3][0]+in[i].y*m4x4[3][1]+in[i].z*m4x4[3][2]+m4x4[3][3]);
                    out[i].w=w;
                    if (!(w==1)) {
                        out[i].x/=w;
                        out[i].y/=w;
                        out[i].z/=w;
                    }
            }
            
        }
        
        public void transformNormal(Vertex3D in[], Vertex3D out[], int length)
        {
            float m3x3[][] = new float[3][3];
            float det; //3x3 determinant
            //2x2 determinants det2x2[i][j] with row i column j eliminated
            float[][] det2x2=new float[3][3];
            
            //set m3x3 to A'
            for (int i=0;i<3;i++) {
                for (int j=0;j<3;j++){
                    m3x3[i][j]=m4x4[i][j];
                }
            }  
            
            
            //calculate inverse of A'
            
            //calculate determinant
            det = m3x3[0][0]*m3x3[1][1]*m3x3[2][2]+m3x3[0][1]*m3x3[1][2]*m3x3[2][0]+m3x3[0][2]*m3x3[1][0]*m3x3[2][1]
                    -m3x3[0][2]*m3x3[1][1]*m3x3[2][0]-m3x3[1][2]*m3x3[2][1]*m3x3[0][0]-m3x3[2][2]*m3x3[0][1]*m3x3[1][0];
            det2x2[0][0] = m3x3[1][1]*m3x3[2][2]-m3x3[1][2]*m3x3[2][1];
            det2x2[0][1] = m3x3[1][0]*m3x3[2][2]-m3x3[2][0]*m3x3[1][2];
            det2x2[0][2] = m3x3[1][0]*m3x3[2][1]-m3x3[1][2]*m3x3[2][0];
            det2x2[1][0] = m3x3[0][1]*m3x3[2][2]-m3x3[0][2]*m3x3[2][1];
            det2x2[1][1] = m3x3[0][0]*m3x3[2][2]-m3x3[2][0]*m3x3[0][2];
            det2x2[1][2] = m3x3[0][0]*m3x3[2][1]-m3x3[0][1]*m3x3[2][0];
            det2x2[2][0] = m3x3[0][1]*m3x3[1][2]-m3x3[1][1]*m3x3[0][2];
            det2x2[2][1] = m3x3[0][0]*m3x3[1][2]-m3x3[1][0]*m3x3[0][2];
            det2x2[2][2] = m3x3[0][0]*m3x3[1][1]-m3x3[1][0]*m3x3[0][1];
            int coeff=1;
            
             
            
            for (int i=0;i<3;i++) {
                for (int j=0;j<3;j++) {
                    if(((i+j)%2)==0){coeff=1;}
                    else {coeff=-1;}
                    m3x3[i][j]=(coeff*det2x2[j][i])/det;
                }
            }
            
            //transpose m3x3
            for (int i=0;i<3;i++) {
                for (int j=0;j<3;j++) {
                    float tmp=m3x3[i][j];
                    m3x3[i][j]=m3x3[j][i];
                    m3x3[j][i]=tmp;
                }
            }
            
                   
            
            //transform normals with m3x3 which = (((A')-1)T)
            for (int i=0;i<length;i++){
                    out[i].nx = m3x3[0][0]*in[i].nx+m3x3[0][1]*in[i].ny+m3x3[0][2]*in[i].nz;
                    out[i].ny = m3x3[1][0]*in[i].nx+m3x3[1][1]*in[i].ny+m3x3[1][2]*in[i].nz;
                    out[i].nz = m3x3[2][0]*in[i].nx+m3x3[2][1]*in[i].ny+m3x3[2][2]*in[i].nz;       
            }
            
        }

        public final void compose(Matrix3D src) {// this = this * src
            int r,c;
            float[][] tmp = new float[4][4];
            for (r=0;r<4;r++) {
                for (c=0;c<4;c++) {
                    tmp[r][c]=this.m4x4[r][0]*src.m4x4[0][c]+
                            this.m4x4[r][1]*src.m4x4[1][c]+
                            this.m4x4[r][2]*src.m4x4[2][c]+
                            this.m4x4[r][3]*src.m4x4[3][c];
                }
            }
            for (r=0;r<4;r++) {
                for (c=0;c<4;c++) {
                    this.m4x4[r][c]=tmp[r][c];
                }
            }
        }//compose
        
        public void loadIdentity() { // this = identity
            Matrix3D tmp = new Matrix3D();
            this.m4x4 = tmp.m4x4;
        }
        
        public void translate(float tx, float ty, float tz){// this = this * t
            Matrix3D tmp = new Matrix3D();
            tmp.m4x4[0][3]=tx;
            tmp.m4x4[1][3]=ty;
            tmp.m4x4[2][3]=tz;
            this.compose(tmp);
        }
        
        public void scale(float sx, float sy, float sz){ // this = this * scale
            Matrix3D tmp = new Matrix3D();
            tmp.m4x4[0][0]=sx;
            tmp.m4x4[1][1]=sy;
            tmp.m4x4[2][2]=sz;
            this.compose(tmp);
        }
        
        public void skew(float kxy, float kxz, float kyz){ // this = this * skew
            Matrix3D tmp = new Matrix3D();
            tmp.m4x4[0][1]=kxy;
            tmp.m4x4[0][2]=kxz;
            tmp.m4x4[1][2]=kyz;
            this.compose(tmp);
        }
        
        public void rotate(float ax, float ay, float az, float angle){ // this = this * rotate
            Matrix3D identity = new Matrix3D();
            Matrix3D symmetric = new Matrix3D();
            Matrix3D skew = new Matrix3D();
            Matrix3D tmp = new Matrix3D();
            
            float t=(float)Math.sqrt(ax*ax+ay*ay+az*az);
            ax/=t;
            ay/=t;
            az/=t;
            
            symmetric.m4x4[0][0] = ax*ax;
            symmetric.m4x4[0][1] = ax*ay;
            symmetric.m4x4[0][2] = ax*az;
            symmetric.m4x4[1][0] = ay*ax;
            symmetric.m4x4[1][1] = ay*ay;
            symmetric.m4x4[1][2] = ay*az;
            symmetric.m4x4[2][0] = az*ax;
            symmetric.m4x4[2][1] = az*ay;
            symmetric.m4x4[2][2] = az*az;
            
            
            skew.m4x4[0][0]=0;
            skew.m4x4[0][1]=-az;
            skew.m4x4[0][2]=ay;
            skew.m4x4[1][0]=az;
            skew.m4x4[1][1]=0;
            skew.m4x4[1][2]=-ax;
            skew.m4x4[2][0]=-ay;
            skew.m4x4[2][1]=ax;
            skew.m4x4[2][2]=0;
            
            
            for (int r=0;r<3;r++) {
                for (int c=0;c<3;c++) {
                    skew.m4x4[r][c] *= Math.sin(angle); 
                    symmetric.m4x4[r][c] *= (1-Math.cos(angle));
                    identity.m4x4[r][c] *= Math.cos(angle);
                    tmp.m4x4[r][c] = skew.m4x4[r][c]+symmetric.m4x4[r][c]+identity.m4x4[r][c];
                }
            }
            this.compose(tmp);
        }//rotate
        
        public void lookAt(float eyex, float eyey, float eyez,
                            float atx,  float aty,  float atz,
                            float upx,  float upy,  float upz){ // this = this * lookat
            float lx,ly,lz,normlx,normly,normlz;
            float rx,ry,rz,normrx,normry,normrz;
            float ux,uy,uz,normux,normuy,normuz;
            Matrix3D tmp = new Matrix3D();
            
            lx = atx-eyex;
            ly = aty-eyey;
            lz = atz-eyez;
            
            normlx = lx/(float)Math.sqrt(lx*lx+ly*ly+lz*lz);
            normly = ly/(float)Math.sqrt(lx*lx+ly*ly+lz*lz);
            normlz = lz/(float)Math.sqrt(lx*lx+ly*ly+lz*lz);
            
            //r = l x up
            rx = -lz*upy+ly*upz;
            ry = lz*upx-lx*upz;
            rz = -ly*upx+lx*upy;
            
            normrx = rx/(float)Math.sqrt(rx*rx+ry*ry+rz*rz);
            normry = ry/(float)Math.sqrt(rx*rx+ry*ry+rz*rz);
            normrz = rz/(float)Math.sqrt(rx*rx+ry*ry+rz*rz);
            
            //u = r x l
            ux = -rz*ly+ry*lz;
            uy = rz*lx-rx*lz;
            uz = -ry*lx+rx*ly;
            
            normux = ux/(float)Math.sqrt(ux*ux+uy*uy+uz*uz);
            normuy = uy/(float)Math.sqrt(ux*ux+uy*uy+uz*uz);
            normuz = uz/(float)Math.sqrt(ux*ux+uy*uy+uz*uz);
            
            tmp.m4x4[0][0]=normrx;
            tmp.m4x4[0][1]=normry;
            tmp.m4x4[0][2]=normrz;
            
            tmp.m4x4[1][0]=normux;
            tmp.m4x4[1][1]=normuy;
            tmp.m4x4[1][2]=normuz;
            
            tmp.m4x4[2][0]=-normlx;
            tmp.m4x4[2][1]=-normly;
            tmp.m4x4[2][2]=-normlz;
            
            tmp.m4x4[3][0]=0;
            tmp.m4x4[3][1]=0;
            tmp.m4x4[3][2]=0;
            tmp.m4x4[3][3]=1;
            
            tmp.m4x4[0][3]=-(normrx*eyex+normry*eyey+normrz*eyez);
            tmp.m4x4[1][3]=-(normux*eyex+normuy*eyey+normuz*eyez);
            tmp.m4x4[2][3]=(normlx*eyex+normly*eyey+normlz*eyez);
            
            this.compose(tmp);
        }//lookAt
        
        
        
        
        //
        // Assume the following projection transformations
        // transform points into the canonical viewing space
        //
        public void perspective(float left, float right,
                            float bottom, float top, 
                            float near, float far){
            
            Matrix3D tmp = new Matrix3D();
            tmp.m4x4[0][0]=2*near/(right-left);
            tmp.m4x4[1][1]=2*near/(bottom-top);
            tmp.m4x4[2][2]=(far+near)/(far-near);
            tmp.m4x4[3][3]=0;
            tmp.m4x4[3][2]=1;
            tmp.m4x4[0][2]=-(right+left)/(right-left);
            tmp.m4x4[1][2]=-(bottom+top)/(bottom-top);
            tmp.m4x4[2][3]=-2*far*near/(far-near);
            
            this.compose(tmp);
        }//perspective
        
        public void orthographic(float left, float right,
                                float bottom, float top,
                                float near, float far)
        {
                                    
            Matrix3D tmp = new Matrix3D();
            tmp.m4x4[0][0]=2/(right-left);
            tmp.m4x4[1][1]=2/(bottom-top);
            tmp.m4x4[2][2]=2/(far-near);
            tmp.m4x4[0][3]=-(right+left)/(right-left);
            tmp.m4x4[1][3]=-(bottom+top)/(bottom-top);
            tmp.m4x4[2][3]=-(far+near)/(far-near);
            this.compose(tmp);
        }//orthographic
    }//class Matrix3D