import Point3D;

    public class Matrix3D {
      //
      // Attributes
      //
      protected float matrix[][];

        //
        // Constructors
        //
        public Matrix3D()                  // initialize with identity transform
      {
	matrix = new float[4][4];

	loadIdentity();
      }

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

	for(int j=0;j<4;j++)
	  for(int i=0;i<4;i++)
	    matrix[j][i] = copy.matrix[j][i];
      }

        public Matrix3D(Raster r)          // initialize with a mapping from
      {                                    // canonical space to screen space
	matrix = new float[4][4];

	loadIdentity();

	//  scale to raster size
	float zmax = 1.0f;
	matrix[0][0] = r.getWidth()/2.0f;
	matrix[1][1] = r.getHeight()/2.0f;
	matrix[2][2] = zmax/2.0f;

	matrix[0][3] = r.getWidth()/2.0f;
	matrix[1][3] = r.getHeight()/2.0f;
	matrix[2][3] = zmax/2.0f;
      }

        //
        // 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 counter=start;counter<start+length;counter++)
	  {
	    out[counter].x = in[counter].x * matrix[0][0] + in[counter].y * matrix[0][1] + in[counter].z * matrix[0][2] + matrix[0][3];
	    out[counter].y = in[counter].x * matrix[1][0] + in[counter].y * matrix[1][1] + in[counter].z * matrix[1][2] + matrix[1][3];
	    out[counter].z = in[counter].x * matrix[2][0] + in[counter].y * matrix[2][1] + in[counter].z * matrix[2][2] + matrix[2][3];

            float fw = in[counter].x * matrix[3][0] + in[counter].y * matrix[3][1] + in[counter].z * matrix[3][2] + matrix[3][3];
            if(fw != 0.0f)  // protect against infinity 
	      {
                out[counter].x /= fw;
                out[counter].y /= fw;
                out[counter].z /= fw;
	      }
	  }
      }

        public final void compose(Matrix3D src)                        // this = this * src
      {
	Matrix3D temp = new Matrix3D(this);

	for(int i=0;i<4;i++)
	  for(int j=0;j<4;j++)
	    {
	      matrix[i][j] = 0.0f;
	      for(int k=0;k<4;k++)
		matrix[i][j] += temp.matrix[i][k] * src.matrix[k][j];
	    }
      }

        public void loadIdentity()                                     // this = identity
      {
	   for(int j=0;j<4;j++)
	     for(int i=0;i<4;i++)
	       matrix[j][i] = 0.0f;

	   for(int i=0;i<4;i++)
	     matrix[i][i] = 1.0f;
      }

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

        public void skew(float kxy, float kxz, float kyz)              // this = this * skew
      {

	//  PLEASE NOTE!   I thought there was some ambiguity between a 'shear' type
	//   matrix and a 'skew' matrix used in the rotate transformation- because of the
	//   arguments labeled kxy, etc, I chose to write this transform as the 'shear'
	//   matrix off of Lecture12/Slide24.html.

        Matrix3D tmpMat = new Matrix3D();
	tmpMat.matrix[0][1] = kxy;
	tmpMat.matrix[0][2] = kxz;
	tmpMat.matrix[1][2] = kyz;
        this.compose(tmpMat);
      }
	
        public void rotate(float ax, float ay, float az, float angle)  // this = this * rotate 
      {
	Matrix3D sym = new Matrix3D();  //  will be identity already

	//  First the symmetric
	sym.matrix[0][0] =  ax*ax * (1.0f - (float)Math.cos(angle));
	sym.matrix[0][1] =  ax*ay * (1.0f - (float)Math.cos(angle));
	sym.matrix[0][2] =  ax*az * (1.0f - (float)Math.cos(angle));
	sym.matrix[1][0] =  ay*ax * (1.0f - (float)Math.cos(angle));
	sym.matrix[1][1] =  ay*ay * (1.0f - (float)Math.cos(angle));
	sym.matrix[1][2] =  ay*az * (1.0f - (float)Math.cos(angle));
	sym.matrix[2][0] =  az*ax * (1.0f - (float)Math.cos(angle));
	sym.matrix[2][1] =  az*ay * (1.0f - (float)Math.cos(angle));
	sym.matrix[2][2] =  az*az * (1.0f - (float)Math.cos(angle));

	//  Then the skew
	sym.matrix[0][1] =  sym.matrix[0][ 1] + -az * (float)Math.sin(angle);
	sym.matrix[0][2] =  sym.matrix[0][ 2] + ay * (float)Math.sin(angle);
	sym.matrix[1][0] =  sym.matrix[1][ 0] + az * (float)Math.sin(angle);
	sym.matrix[1][2] =  sym.matrix[1][ 2] + -ax * (float)Math.sin(angle);
	sym.matrix[2][0] =  sym.matrix[2][ 0] + -ay * (float)Math.sin(angle);
	sym.matrix[2][1] =  sym.matrix[2][ 1] + ax * (float)Math.sin(angle);

	//  Finally the identity
	sym.matrix[0][0] =  sym.matrix[0][ 1] + (float)Math.cos(angle);
	sym.matrix[1][1] =  sym.matrix[0][2] + (float)Math.cos(angle);
	sym.matrix[2][2] =  sym.matrix[1][0] + (float)Math.cos(angle);

	//  Multiply
	this.compose(sym);
      }

        public void lookAt(float eyex, float eyey, float eyez,
                         float atx,  float aty,  float atz,
                         float upx,  float upy,  float upz)          // this = this * lookat
      {
	//  First make a normal vector along the look-at 
	Point3D ptTemp = new Point3D(0.0f, 0.0f, 0.0f);
	ptTemp.x = atx - eyex;
	ptTemp.y = aty - eyey;
	ptTemp.z = atz - eyez;
	float fNormal = (float)Math.sqrt((double)(ptTemp.x*ptTemp.x + ptTemp.y*ptTemp.y + ptTemp.z*ptTemp.z));
	ptTemp.x /= fNormal;
	ptTemp.y /= fNormal;
	ptTemp.z /= fNormal;

	//  Now make the 'right' vector
	Point3D ptTempRight = new Point3D(0.0f, 0.0f, 0.0f);
	ptTempRight.x = ptTemp.y * upz - ptTemp.z * upy;
	ptTempRight.y = ptTemp.z * upx - ptTemp.x * upz;
	ptTempRight.z = ptTemp.x * upy - ptTemp.y * upx;
	fNormal = (float)Math.sqrt((double)(ptTempRight.x*ptTempRight.x + ptTempRight.y*ptTempRight.y + ptTempRight.z*ptTempRight.z));
	ptTempRight.x /= fNormal;
	ptTempRight.y /= fNormal;
	ptTempRight.z /= fNormal;
	
	//  Now make the local 'up' vector
	Point3D ptTempUp = new Point3D(0.0f, 0.0f, 0.0f);
	ptTempUp.x = ptTempRight.y * ptTemp.z - ptTempRight.z * ptTemp.y;
	ptTempUp.y = ptTempRight.z * ptTemp.x - ptTempRight.x * ptTemp.z;
	ptTempUp.z = ptTempRight.x * ptTemp.y - ptTempRight.y * ptTemp.x;
	fNormal = (float)Math.sqrt((double)(ptTempUp.x*ptTempUp.x + ptTempUp.y*ptTempUp.y + ptTempUp.z*ptTempUp.z));
	ptTempUp.x /= fNormal;
	ptTempUp.y /= fNormal;
	ptTempUp.z /= fNormal;

	//  Put it into our matrix
	Matrix3D matTemp = new Matrix3D();
	matTemp.matrix[0][0] = ptTempRight.x;
	matTemp.matrix[0][1] = ptTempRight.y;
	matTemp.matrix[0][2] = ptTempRight.z;
	matTemp.matrix[0][3] = -(ptTempRight.x * eyex + ptTempRight.y * eyey + ptTempRight.z * eyez);
	matTemp.matrix[1][0] = ptTempUp.x;
	matTemp.matrix[1][1] = ptTempUp.y;
	matTemp.matrix[1][2] = ptTempUp.z;
	matTemp.matrix[1][3] = -(ptTempUp.x * eyex + ptTempUp.y * eyey + ptTempUp.z * eyez);
	matTemp.matrix[2][0] = -ptTemp.x;
	matTemp.matrix[2][1] = -ptTemp.y;
	matTemp.matrix[2][2] = -ptTemp.z;
	matTemp.matrix[2][3] = ptTemp.x * eyex + ptTemp.y * eyey + ptTemp.z * eyez;
	compose(matTemp);
      }

        //
        // 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 matTemp = new Matrix3D();
	matTemp.matrix[0][0] = 2.0f*near/(right - left);
	matTemp.matrix[1][1] = 2.0f*near/(bottom - top);
	matTemp.matrix[2][2] = (far + near)/(far - near);
	matTemp.matrix[0][2] = -(right + left)/(right-left);
	matTemp.matrix[1][2] = -(bottom + top)/(bottom-top);
	matTemp.matrix[2][3] = -far*near*2.0f/(far-near);
	matTemp.matrix[3][3] = 0.0f;
	matTemp.matrix[3][2] = 1.0f;
	compose(matTemp);
      }

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