Triangle Triangle::standardize(void) const
{

	//  we assume that the triangle is valid;  i.e. the vertices
	//  distinct -- if not this will BOMB!
	float myMatrix[16];
	int i,j,k;


	Triangle standard=*this;
	glPushMatrix();			// save the current matrix

	// translate so the first vertex is at the origin
	standard[1] -= standard[0];
	standard[2] -= standard[0];
	standard[0] = Tuple(0,0,0);


	// rotate to align last edge with the x axis
	Tuple x(1,0,0);
	Tuple w = standard[2];	// set w to the last edge
	w.normalize();			// normalize it
	Tuple v = w.cross(x);	// find the vector of rotation
	v.normalize();			// normalize it
	float phi = (float) acos((double) w.dot(x));
							// find the angle of rotation in radians
	phi *= 360.f/(2.0f*3.14159f);
							// convert to degrees

	glLoadIdentity();		// load the identity matrix	
	glRotatef(phi,v[0],v[1],v[2]);
							// openGL computes the rotation transform
	glGetFloatv(GL_MODELVIEW_MATRIX,myMatrix);
							// we'll grab the matrix transform
	// update vertices
	for (i=0; i<3; i++) { // for each vertex
		Tuple tmp(0,0,0); 
		for (j=0; j<3; j++) { // for each row of the matrix 
			for (k=0; k<3; k++) {  // for each column
				tmp[j] += myMatrix[j+4*k]*standard[i][k];  // the matrix is column-major
			}
		}
		standard[i]=tmp;
	}


	// rotate to bring second vertex into y=0 plane
	w=standard[1];			// set w to the lie in the first edge
	w[0]=0.0f;				// project into the x=0 plane
	float sigma= (float) acos((double) (w[2]/w.length()));
							// compute the angle of rotation in radians
	sigma *= 360.f/(2.0f*3.14159f);
							// convert to degrees
	if (w[1]<0)				// should we rotate about + or - x axis
		x *= -1.0f;
	glLoadIdentity();		// reset the matrix
	glRotatef(sigma,x[0],x[1],x[2]);
							// openGL will compute the rotation transform
	glGetFloatv(GL_MODELVIEW_MATRIX,myMatrix);
							// grab the matrix
	// update vertices
	for (i=0; i<3; i++) { // for each vertex
		Tuple tmp(0,0,0); 
		for (j=0; j<3; j++) { // for each row of the matrix 
			for (k=0; k<3; k++) {  // for each column
				tmp[j] += myMatrix[j+4*k]*standard[i][k];  // the matrix is column-major
			}
		}
		standard[i]=tmp;
	}


	// now compute the single forward (standardizing) transform
	float forwardTransform[16];
	glLoadIdentity();
	glRotatef(sigma,x[0],x[1],x[2]);
	glRotatef(phi,v[0],v[1],v[2]);
	glTranslatef(this->u[0][0],this->u[0][1], this->u[0][2]);
	glGetFloatv(GL_MODELVIEW_MATRIX,forwardTransform);

	// now compute the single reverse transform
	float reverseTransform[16];
	glLoadIdentity();
	glTranslatef(this->u[0][0],this->u[0][1], this->u[0][2]);
	glRotatef(-phi,v[0],v[1],v[2]);
	glRotatef(-sigma,x[0],x[1],x[2]);
	glGetFloatv(GL_MODELVIEW_MATRIX,reverseTransform);


	glPopMatrix(); // restore current matrix


	return standard;

}