public class Matrix3D {
        
	    ///////////////////////// Variables ////////////////////////////	

		public float matrix[][];
		


		///////////////////////// Constructors ////////////////////////////			
		
		/**
		 * Initialize with identity transform
		 * 
		 */
		public Matrix3D() {
 
			//Test if matrix is initialized with zeros!
			matrix = new float[4][4];
			loadIdentity();
				  
		}                  
		
		/**
		 * Initialize with copy of source
		 * 
		 */
		public Matrix3D(Matrix3D copy) {
							  
			matrix = new float[4][4];
			matrix = copy.matrix;
			
		}
		

		/**
		 * Initialize with copy of source matrix
		 * 
		 */
		public Matrix3D(float[][] copy) {
							  
			matrix = new float[4][4];
			matrix = copy;
			
		}

		
        /**
         * Initialize with a mapping from canonical space to screen space
         */
		public Matrix3D(Raster r) {
			
			this();
			
			matrix = new float[4][4];
			
			matrix[0][0] = (float)r.width/2;
			matrix[0][3] = (float)r.width/2;
			matrix[1][1] = (float)-r.height/2;
			matrix[1][3] = (float)r.height/2;
			matrix[3][3] = 1.0f;
			
            
			matrix[0][0]=r.width/2;
            matrix[0][3]=r.width/2;
            matrix[1][1]=-r.height/2;
            matrix[1][3]=r.height/2;
            matrix[2][2]=(float).5;  
            matrix[2][3]=(float).5;
            matrix[3][3]=1;
			
				
		}

		///////////////////////// Methods ////////////////////////////					
		
        /**
         * 
         * set element [j][i] to value
         */
		public void set(int i, int j, float value) {
			matrix[i][j] = value;
		}
		
        /**
         * 
         * return element [j][i]
         */
		public float get(int i, int j) {
			return matrix[i][j];
		}

        /**
         *  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
         */
		public void transform(Point3D in[], Point3D out[], int start, int length) {
			
			//Note: we must remember to use the 4th dimension in our transformation matrix,
			//keeping in mind that the 4th dimension coordinate of the incoming point is always 1.
			//Divide by W to guarantee correct perspective.
			float w;
			
			for (int i=0; i<length; i++) {
				
				w = matrix[3][0]*in[start+i].x +
				matrix[3][1]*in[start+i].y +
				matrix[3][2]*in[start+i].z +
				matrix[3][3];
				
				out[start+i].x = ((matrix[0][0] * in[start+i].x) +
								  (matrix[0][1] * in[start+i].y) +
								  (matrix[0][2] * in[start+i].z) +
								   matrix[0][3]) / w;
				
				out[start+i].y = ((matrix[1][0] * in[start+i].x) +
								  (matrix[1][1] * in[start+i].y) +
								  (matrix[1][2] * in[start+i].z) +
								   matrix[1][3]) / w;
				
				out[start+i].z = ((matrix[2][0] * in[start+i].x) +
								  (matrix[2][1] * in[start+i].y) +
								  (matrix[2][2] * in[start+i].z) +
								   matrix[2][3]) / w;
			}
				
			
		}
		
		
		/**
		 * Multiply the current matrix by another
		 */
		void multiply (float[][] matrix2) {
			
			float mult[][] = new float[4][4];
		
			for (int i=0; i<matrix2.length; i++)
				for (int j=0; j<(matrix2[0]).length; j++)
					mult[i][j] = (matrix2[i][0] * matrix[0][j]) +
								 (matrix2[i][1] * matrix[1][j]) +
								 (matrix2[i][2] * matrix[2][j]) +
								 (matrix2[i][3] * matrix[3][j]);
			
			matrix = mult;

		}


        /**
         * this = this * src
         */
		public final void compose(Matrix3D src) {
			
			src.multiply(matrix);
			matrix = src.matrix;
		} 
		
        /**
         * this = this * src
         */
		public final void compose(float src[][]) {
			
			compose(new Matrix3D(src));
		}		
        
		/**
		 * this = identity
		 * */
		public void loadIdentity() {
		
			for (int i=0; i<4; i++)
				for (int j=0; j<4; j++) 
					if (i==j) matrix[i][j] = 1.0f;
					else matrix[i][j] = 0.0f;
		}
		
		/**
		 * this = identity
		 * */
		public void loadIdentity(float[][] matrix) {
		
			for (int i=0; i<4; i++)
				for (int j=0; j<4; j++) 
					if (i==j) matrix[i][j] = 1.0f;
					else matrix[i][j] = 0.0f;
		}

			
        /**
         * this = this * t
         */
		
		public void translate(float tx, float ty, float tz) {
			
			matrix[0][3] += tx;
			matrix[1][3] += ty;
			matrix[2][3] += tz;
			
		}
		
        /**
         * this = this * scale
		 */
		public void scale(float sx, float sy, float sz) {

			//All we have to do here is create an indentity matrix whose diagonal
			//coefficients are multiplied by the scaling factor
			float scalingMatrix[][] = new float[4][4];
			scalingMatrix[0][0] = sx;
			scalingMatrix[1][1] = sy;
			scalingMatrix[2][2] = sz;
			scalingMatrix[3][3] = 1.0f;
			
			compose(scalingMatrix);
		}

		
		/**
		 * this = this * skew
		 */
		public void skew(float kxy, float kxz, float kyz) {
			
			float skewingMatrix[][] = new float[4][4];
			loadIdentity(skewingMatrix);
			skewingMatrix[1][0] = kxy;
			skewingMatrix[2][0] = kxz;
			skewingMatrix[2][1] = kyz;
			
			compose(skewingMatrix);			
			
		}
		
        /**
         * this = this * rotate 
         */
		public void rotate(float ax, float ay, float az, float angle) {
			
			float symmetrice[][] = new float[3][3];
			float cos = (float)Math.cos(angle);
			float sin = (float)Math.sin(angle);
			float multip = 1-cos;
  
  
			symmetrice[0][0] = ax * ax * multip;
			symmetrice[0][1] = ax * ay * multip;
			symmetrice[0][2] = ax * az * multip;
			symmetrice[1][0] = ax * ay * multip;
			symmetrice[1][1] = ay * ay * multip;
			symmetrice[1][2] = ay * az * multip;
			symmetrice[2][0] = ax * az * multip;
			symmetrice[2][1] = ay * az * multip;
			symmetrice[2][2] = az * az * multip;

			float skew[][] = new float[3][3];
			for (int i=0; i<=2; i++) {
				skew[i][i] = 0.0f;
			}
			  
			skew[0][1] = -az * sin;
			skew[0][2] = ay * sin;
			skew[1][0] = az * sin;
			skew[1][2] = -ax * sin;
			skew[2][0] = -ay * sin;
			skew[2][1] = ax * sin;
			  
			float id[][] = new float[3][3];
			for (int i=0; i<=2; i++) {
				for (int j=0; j<=2; j++) {
				if (i==j)
					id[j][i] = cos;
				else
					id[j][i]=0.0f;
				}
			}
			  
			  
			Matrix3D rotationMatrix = new Matrix3D();
			 
			for (int i=0; i<=2; i++) {
				for (int j=0; j<=2; j++) {
					rotationMatrix.set(i, j, (symmetrice[j][i] + skew[j][i] + id[j][i]));
				}
			}
  
			compose(rotationMatrix);
  
  
	
		}

		
		
		
        /**
         * this = this * lookat
         */
		public void lookAt(float eyex, float eyey, float eyez,
                           float atx,  float aty,  float atz,
						   float upx,  float upy,  float upz) {
			
			float lookAtMatrix[][] = new float[4][4];
			float R[][] = new float[4][4];
			float T[][] = new float[4][4];
    
			Point3D eye, up, n, u, v;
			

			eye = new Point3D (eyex, eyey, eyez);
			up = new Point3D (upx, upy, upz);
			n = new Point3D ((atx - eyex),
							 (aty - eyey),
							 (atz - eyez));
			u = new Point3D ();
			v = new Point3D ();

			
			u.crossProduct(n, up);

			v.crossProduct(u, n);

			n.normalize();
			v.normalize();
			u.normalize();
			
		
			//M = R*T
			R[0][0] = u.x;
			R[0][1] = u.y;
			R[0][2] = u.z;
			R[1][1] = v.x;
			R[1][1] = -v.y;
			R[1][2] = v.z;
			R[2][0] = n.x;
			R[2][1] = n.y;
			R[2][2] = n.z;
			R[3][3] = 1.0f;
			
			T[0][0] = 1.0f;
			T[0][3] = -eye.x;
			T[1][1] = 1.0f;
			T[1][3] = -eye.y;
			T[2][2] = 1.0f;
			T[2][3] = -eye.z;
			T[3][3] = 1.0f;
			
			
			lookAtMatrix = multiply(R,T);
			compose(lookAtMatrix);
			
			System.out.println("up, v, u, n, eye vectors:");
			printVector(up);
			printVector(v);
			printVector(u);
			printVector(n);
			printVector(eye);
			
			
		}		

        //
        // Assume the following projection transformations
        // transform points into the canonical viewing space
        //
        /**
         * this = this * persp
         */
		public void perspective(float left, float right,                 
                                float bottom, float top,
								float near, float far) {
			
			float perspMatrix[][] = new float[4][4];
			
			perspMatrix[0][0] = (2.0f * near) / (right - left);
			perspMatrix[0][2] = (right + left) / (right - left);

		    perspMatrix[1][1] = (2.0f * near) / (bottom - top);
			perspMatrix[1][2] = - (bottom + top) / (bottom - top);

			perspMatrix[2][2] = - (far + near) / (far - near);
			perspMatrix[2][3] = - (2.0f * far * near) / (far - near);
			
			perspMatrix[3][2] = -1.0f;
			
			compose(perspMatrix);
			
		}
		
		
        /**
         * this = this * ortho
         */
		public void orthographic(float left, float right,                
                                 float bottom, float top,
								 float near, float far) {
			
			float orthoMatrix[][] = new float[4][4];
			
			
			orthoMatrix[0][0] = 2.0f / (right - left);
			orthoMatrix[0][3] = - (right + left) / (right - left);
		
			orthoMatrix[1][1] = 2.0f / (bottom - top);
			orthoMatrix[1][3] = (bottom + top) / (top - bottom);
		
			orthoMatrix[2][2] = 2.0f / (far - near);
			orthoMatrix[2][3] = - (far + near) / (far - near);	
			
			orthoMatrix[3][3] = 1.0f;
			
			compose(orthoMatrix);
			

  
		}
		
		
		public void printMatrix(float[][] matrix) {
			
			for (int i=0; i<4; i++)
				System.out.println("[ " + matrix[i][0] + " " + matrix[i][1] + " " + matrix[i][2] + " " + matrix[i][3] + " ]");
		
			
		}
		
		public void printVector(Point3D p) {
			
				System.out.println("[ " + p.x + " " + p.y + " " + p.z + " ]");
		
			
		}
		
		/**
		 * Multiply a matrix by another
		 */
		public float[][] multiply (float[][] matrix1, float[][] matrix2) {
			
			float mult[][] = new float[4][4];
			
			for (int i=0; i<matrix1.length; i++)
				for (int j=0; j<(matrix1[0]).length; j++)
					mult[i][j] = (matrix1[i][0] * matrix2[0][j]) +
								 (matrix1[i][1] * matrix2[1][j]) +
								 (matrix1[i][2] * matrix2[2][j]) +
								 (matrix1[i][3] * matrix2[3][j]);
			
			return mult;
		}
		
		
	}