// COMP 236 Programming Assignment 3
// Adrian Ilie

#include <string.h>
#include <GL/glut.h>
#include <GL/glut.h>
#include <fstream.h>
#include <iostream.h>
#include <stdlib.h>
#include <math.h>

#define NUM_PATCHES 32
#define C_X 0
#define C_Y 1
#define C_Z 2

int do_draw_normals = 1;
int do_wire_frame = 0;
int main_window_id;
float scale = 1;
const float material[]={0.3f, 0.2f, 0.4f, 100.0f};
double max_coord[3]; // max coordinates for view volume
double**** patches=NULL; // patches read from file
double**** points=NULL;	// curve points computed from patches using the Bezier formulas
double**** normals=NULL; // normals at each point
double** mat_bezier;// Bezier matrix, equal to its transpose

const char* file_name = "teapot.txt";
int num_curves = 10;

float xy_ratio=1.3;	// x/y aspect ratio

float light0_intensity = 1.0;
float light1_intensity = 0.4;
float light2_intensity = 0.5;
float pos_teapot[] = { 0.0, 0.0, 120 };
float rotX = 0,rotY = 0;
float pi = 3.141529;

// lights - from Paul Rademacher's examples
GLfloat light0_ambient[] =  {0.1f, 0.1f, 0.3f, 1.0f};
GLfloat light0_diffuse[] =  {.6f, .6f, 1.0f, 1.0f};
GLfloat light0_position[] = {1.0f, 1.0f, -1.0f, 0.0f};

GLfloat light1_ambient[] =  {0.1f, 0.1f, 0.3f, 1.0f};
GLfloat light1_diffuse[] =  {.9f, .6f, 0.0f, 1.0f};
GLfloat light1_position[] = {-1.0f, 1.0f, 1.0f, 0.0f};

GLfloat light2_ambient[] =  {1.0f, 1.0f, 1.0f, 1.0f};
GLfloat light2_diffuse[] =  {1.0f, 1.0f, 1.0f, 1.0f};
GLfloat light2_position[] = {1.0f, 1.0f, 1.0f, 0.0f};

GLfloat light_specular[] = {1.0f, 1.0f, 1.0f, 1.0f};


double** mat_multiply(double** mat1, double** mat2, int m1, int n, int m2)
{
	int i,j,k;
	double** r = new double*[m1];
	for (i=0; i<m1; i++) r[i] = new double[m2];
	for (i=0; i<m1; i++)
		for (j=0; j<m2; j++)
		{	r[i][j] = 0;
			for (k=0; k<n; k++)
				r[i][j] += mat1[i][k]*mat2[k][j];
		}
	return r;
}

void alloc_matrices()
{
	patches = new double***[NUM_PATCHES];
	points = new double***[NUM_PATCHES];
	normals = new double***[NUM_PATCHES];
	int i,j,k;
	for (i=0; i<NUM_PATCHES; i++) //patches
	{
		patches[i] = new double**[3];
		points[i] = new double**[3];
		normals[i] = new double**[3];
		for (j=0; j<3; j++) //coordinates
		{
			patches[i][j] = new double*[4]; //4*4 control points in a patch
			points[i][j] = new double*[num_curves]; //curve points
			normals[i][j] = new double*[num_curves]; //normals
			for (k=0; k<4; k++) patches[i][j][k] = new double[4]; //4*4 control points in a patch
			for (k=0; k<num_curves; k++)
			{
				points[i][j][k] = new double[num_curves];
				normals[i][j][k] = new double[num_curves];
			}
		}
	}
}

void delete_matrices()
{
	int i,j,k;
	for (i=0; i<NUM_PATCHES; i++) //patches
	{
		for (j=0; j<3; j++) //coordinates
		{
			for (k=0; k<4; k++) delete[] patches[i][j][k]; //control points
			for (k=0; k<num_curves; k++) //curve points
			{
				delete[] points[i][j][k];
				delete[] normals[i][j][k];
			}
			delete[] patches[i][j];
			delete[] points[i][j];
			delete[] normals[i][j];
		}
		delete[] patches[i];
		delete[] points[i];
		delete[] normals[i];
	}
	delete[] patches;		patches=NULL;
	delete[] points;	points=NULL;
	delete[] normals;	normals=NULL;
}

void read_teapot()
{
	ifstream f;
	f.open(file_name);
	if (!f.is_open())
	{	cerr << "Cannot read file " << file_name << endl;
		exit(1);
	}
	for (int i=0; i<NUM_PATCHES; i++) //patches
	{
		for (int j=0; j<4; j++) //control points
		{
			for (int k=0; k<4; k++) //edges
			{
				for (int c=0; c<3; c++) //coordinates
				{
					f >> patches[i][c][j][k];
					max_coord[c] = (abs(patches[i][c][j][k]) < max_coord[c])
						? max_coord[c] : abs(patches[i][c][j][k]);
				}
			}
		}
	}
	f.close();
}

double one_coord(double s, double v, int coord, int p)
{
	double** c;
	double** S = new double*[1];
	S[0] = new double[4];
	S[0][0] = s*s*s;
	S[0][1] = s*s;
	S[0][2] = s;
	S[0][3] = 1;
	double** T = new double*[4];
	for (int i=0; i<4; i++)	T[i] = new double[1];
	T[0][0] = v*v*v;
	T[1][0] = v*v;
	T[2][0] = v;
	T[3][0] = 1;
	c = patches[p][coord];//G
	c = mat_multiply(mat_bezier, c, 4,4,4);//MG
	c = mat_multiply(c, mat_bezier, 4,4,4);//MGM (M=MT)
	c = mat_multiply(S, c, 1,4,4);//UMGM
	c = mat_multiply(c, T, 1,4,1);//UMGMV
	double coordinate = c[0][0];
	for (int j=0; i<4; i++) delete[] T[i];
	delete[] S[0];
	delete[] S;
	delete[] T;
	delete[] c[0];
	delete[] c;
	return coordinate;
}

void compute_patches()
{
	double delta = 1.0/(num_curves-1);
	double s,v;
	for (int p=0; p<NUM_PATCHES; p++) //patches
	{
		s = 0;
		for (int cs = 0; cs < num_curves; cs++) //curves in S
		{
			v = 0;
			for (int ct = 0; ct < num_curves; ct++) //curves in T
			{
				points[p][C_X][cs][ct] = one_coord(s, v, 0, p);
				points[p][C_Y][cs][ct] = one_coord(s, v, 1, p);
				points[p][C_Z][cs][ct] = one_coord(s, v, 2, p);
				v += delta;
			}
			s += delta;
		}
	}
}

void one_normal(double s, double v, int coord, int p, double &vds, double &vdt)
{
	double** c;
	double** MGM;
	double** ds;
	double** dt;
	int i;
	double** S = new double*[1];
	S[0] = new double[4];
	S[0][0] = s*s*s;
	S[0][1] = s*s;
	S[0][2] = s;
	S[0][3] = 1;
	double** T = new double*[4];
	for (i=0; i<4; i++)	T[i] = new double[1];
	T[0][0] = v*v*v;
	T[1][0] = v*v;
	T[2][0] = v;
	T[3][0] = 1;
	double** dS = new double*[1];
	dS[0] = new double[4];
	dS[0][0] = 3*s*s;
	dS[0][1] = 2*s;
	dS[0][2] = 1;
	dS[0][3] = 0;
	double** dT = new double*[4];
	for (i=0; i<4; i++)	dT[i] = new double[1];
	dT[0][0] = 3*v*v;
	dT[1][0] = 2*v;
	dT[2][0] = 1;
	dT[3][0] = 0;
	c = patches[p][coord]; //G
	c = mat_multiply(mat_bezier, c, 4,4,4); //MG
	MGM = mat_multiply(c, mat_bezier, 4,4,4); //MGM (M=MT)
	c = mat_multiply(dS, MGM, 1,4,4); //UMGM
	ds = mat_multiply(c, T, 1,4,1); //UMGMdV
	vds = ds[0][0]; //tangent in S
	c = mat_multiply(S, MGM, 1,4,4); //dUMGM
	dt = mat_multiply(c, dT, 1,4,1); //dUMGMV
	vdt = dt[0][0]; //tangent in T
	for (i=0; i<4; i++)	delete[] T[i];
	delete[] S[0];
	delete[] S;
	delete[] T;
	for (i=0; i<4; i++)	delete[] dT[i];
	delete[] dS[0];
	delete[] dS;
	delete[] dT;
}

double* cross_product(double x1, double x2, double x3,
					 double y1, double y2, double y3)
{
	double* v = new double[3];
	double l=0;
	int i;
	v[0] = x2*y3 - x3*y2;
	v[1] = x3*y1 - x1*y3;
	v[2] = x1*y2 - x2*y1;
	for (i=0; i<3; i++) l += v[i] * v[i];
	l = sqrt(l);
	for (i=0; i<3; i++) v[i] /= l;//normalize
	return v;
}

void compute_normals()
{
	double s,v,x0,x1,x2,y0,y1,y2;
	double *norm;
	double dsv = 1.0/(double)num_curves;

	for (int p=0; p<NUM_PATCHES; p++) //patches
	{
		s = 0;
		for (int cs = 0; cs < num_curves; cs++) //curves in S
		{
			v = 0;
			for (int ct = 0; ct < num_curves; ct++) //curves in T
			{
				one_normal(s,v,C_X,p,x0,y0); //X
				one_normal(s,v,C_Y,p,x1,y1); //Y
				one_normal(s,v,C_Z,p,x2,y2); //Z

				norm = cross_product(x0, x1, x2, y0, y1, y2);
				normals[p][C_X][cs][ct] = norm[C_X]; //nX
				normals[p][C_Y][cs][ct] = norm[C_Y]; //nY
				normals[p][C_Z][cs][ct] = norm[C_Z]; //nZ
				v += dsv;
			}
			s += dsv;
		}
	}
}

void fix_one_normal(int p1, int p2, int c11, int c12, int c21, int c22)
{
	double nx,ny,nz,s;
	if ((abs(floor(points[p1][C_X][c11][c12])-floor(points[p2][C_X][c21][c22]))<5)&&
	    (abs(floor(points[p1][C_Y][c11][c12])-floor(points[p2][C_Y][c21][c22]))<5)&&
	    (abs(floor(points[p1][C_Z][c11][c12])-floor(points[p2][C_Z][c21][c22]))<5)) //if points coincide
	{
		nx = (normals[p1][C_X][c11][c12] + normals[p2][C_X][c21][c22])/2; //compute the mean
		ny = (normals[p1][C_Y][c11][c12] + normals[p2][C_Y][c21][c22])/2;
		nz = (normals[p1][C_Z][c11][c12] + normals[p2][C_Z][c21][c22])/2;
//		s = sqrt(nx * nx + ny * ny + nz * nz); //normalize - looks bad ;)
//		nx /= s;
//		ny /= s;
//		nz /= s;
		normals[p1][C_X][c11][c12] = nx;
		normals[p2][C_X][c21][c22] = nx;
		normals[p1][C_Y][c11][c12] = ny;
		normals[p2][C_Y][c21][c22] = ny;
		normals[p1][C_Z][c11][c12] = nz;
		normals[p2][C_Z][c21][c22] = nz;
	}
}

void fix_edges()
{
	int p1,p2,pt;
	for (p1=0; p1<NUM_PATCHES; p1++) //cycle through pairs of patches
		for (p2=0; p2<NUM_PATCHES; p2++)
		{
			if (p1!=p2) //if the patch is not the same
			for (pt=0; pt<num_curves; pt++)
			{
				fix_one_normal(p1,p2,0,pt,num_curves-1,pt);
				fix_one_normal(p1,p2,pt,0,pt,num_curves-1);
				fix_one_normal(p1,p2,num_curves-1,pt,0,pt);
				fix_one_normal(p1,p2,pt,num_curves-1,pt,0);
				fix_one_normal(p1,p2,pt,num_curves-1,pt,num_curves-1);
				fix_one_normal(p1,p2,num_curves-1,pt,num_curves-1,pt);
				fix_one_normal(p1,p2,pt,0,pt,0);
				fix_one_normal(p1,p2,0,pt,0,pt);
			}
		}
}

void draw_normals()
{
	double x,y,z,nx,ny,nz;
	int p,cs,ct;
	glDisable(GL_LIGHTING);
	for (p=0; p<NUM_PATCHES; p++) //parches
	{
		for (cs = 0; cs < num_curves; cs++) //curves in S
		{
			for (ct = 0; ct < num_curves; ct++) //curves in T
			{
				nx = normals[p][C_X][cs][ct]; //normal
				ny = normals[p][C_Y][cs][ct];
				nz = normals[p][C_Z][cs][ct];
				x = points[p][C_X][cs][ct]; //point
				y = points[p][C_Y][cs][ct];
				z = points[p][C_Z][cs][ct];
				glColor3f(1,0,0);
				glBegin(GL_LINES);
				glVertex3f(x,y,z);
				glVertex3f(x+4*nx,y+4*ny,z+4*nz);
				glEnd();
			}
		}
	}
	glEnable(GL_LIGHTING);
}

void draw_patches()
{
	double x,y,z,nx,ny,nz;
	int p,cs,ct;
	for (p=0; p<NUM_PATCHES; p++) //patches
	{
		for (cs = 0; cs < num_curves-1; cs++) //curves in S
		{
			for (ct = 0; ct < num_curves-1; ct++)  //curves in T
			{
				if(do_wire_frame) glBegin(GL_LINE_STRIP);
					else glBegin(GL_QUADS);
				x = points[p][C_X][cs][ct]; //point 00
				y = points[p][C_Y][cs][ct];
				z = points[p][C_Z][cs][ct];
				nx = normals[p][C_X][cs][ct];
				ny = normals[p][C_Y][cs][ct];
				nz = normals[p][C_Z][cs][ct];
  				glNormal3f(nx,ny,nz);
				glVertex3f(x,y,z);
				x = points[p][C_X][cs+1][ct]; //point 10
				y = points[p][C_Y][cs+1][ct];
				z = points[p][C_Z][cs+1][ct];
				nx = normals[p][C_X][cs+1][ct];
				ny = normals[p][C_Y][cs+1][ct];
				nz = normals[p][C_Z][cs+1][ct];
				glNormal3f(nx,ny,nz);
				glVertex3f(x,y,z);
				x = points[p][C_X][cs+1][ct+1]; //point 11
				y = points[p][C_Y][cs+1][ct+1];
				z = points[p][C_Z][cs+1][ct+1];
				nx = normals[p][C_X][cs+1][ct+1];
				ny = normals[p][C_Y][cs+1][ct+1];
				nz = normals[p][C_Z][cs+1][ct+1];
				glNormal3f(nx,ny,nz);
				glVertex3f(x,y,z);
				x = points[p][C_X][cs][ct+1]; //point 01
				y = points[p][C_Y][cs][ct+1];
				z = points[p][C_Z][cs][ct+1];
				nx = normals[p][C_X][cs][ct+1];
				ny = normals[p][C_Y][cs][ct+1];
				nz = normals[p][C_Z][cs][ct+1];
				glNormal3f(nx,ny,nz);
				glVertex3f(x,y,z);
				glEnd();
			}
		}
	}
}

void Keyboard(unsigned char Key, int x, int y)
{
  switch(Key)
  {
  case 'q':exit(0);break;
  case '1':pos_teapot[0]++;break;
  case '2':pos_teapot[0]--;break;
  case '3':pos_teapot[1]++;break;
  case '4':pos_teapot[1]--;break;
  case '5':pos_teapot[2]++;break;
  case '6':pos_teapot[2]--;break;
  case '7':rotX++;break;
  case '8':rotX--;break;
  case '9':rotY++;break;
  case '0':rotY--;break;
  case 'w':do_wire_frame=!do_wire_frame;break;
  case 'n':do_draw_normals=!do_draw_normals;break;
  };
  
  glutPostRedisplay();
}


void Idle( void )
{
  if ( glutGetWindow() != main_window_id ) glutSetWindow(main_window_id);  
  glutPostRedisplay();
}

void Motion(int x, int y )
{
  glutPostRedisplay(); 
}

void Reshape(int x, int y )
{
  xy_ratio=x/y;
  glutPostRedisplay(); 
}

void Display( void )
{
	glClearColor( 1.0f, 1.0f, 1.0f, 1.0f );
    glClear( GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT );
    glMatrixMode( GL_PROJECTION );
    glLoadIdentity();
    gluPerspective(90, xy_ratio, 0.1, 2*max_coord[2]);
    glMatrixMode( GL_MODELVIEW );
    glLoadIdentity();
    glLightfv(GL_LIGHT0, GL_POSITION, light0_position);
    glLightfv(GL_LIGHT1, GL_POSITION, light1_position);
    glLightfv(GL_LIGHT2, GL_POSITION, light2_position);
    glLoadIdentity();
    glTranslatef( pos_teapot[0], pos_teapot[1], -pos_teapot[2] ); 
	glRotatef(rotX,0,1,1);
	glRotatef(rotY,1,0,1);
    glScalef( scale, scale, scale );
    glPushMatrix();
    glTranslatef(0,0,-max_coord[2]); 
    draw_patches();
    if (do_draw_normals) draw_normals();
    glPopMatrix();
    glutSwapBuffers(); 
}

void my_exit()
{
	delete_matrices();
	exit(1);
}

void main(int argc, char* argv[])
{

	if (argc == 2) num_curves = atoi(argv[1]) + 1;
	mat_bezier = new double*[4];
	for (int i=0; i<4; i++)	mat_bezier[i] = new double[4];
	mat_bezier[0][0] = -1;	mat_bezier[0][1] = 3;	mat_bezier[0][2] = -3;	mat_bezier[0][3] = 1;
	mat_bezier[1][0] = 3;	mat_bezier[1][1] = -6;	mat_bezier[1][2] = 3;	mat_bezier[1][3] = 0;
	mat_bezier[2][0] = -3;	mat_bezier[2][1] = 3;	mat_bezier[2][2] = 0;	mat_bezier[2][3] = 0;
	mat_bezier[3][0] = 1;	mat_bezier[3][1] = 0;	mat_bezier[3][2] = 0;	mat_bezier[3][3] = 0;
	alloc_matrices();
	read_teapot();
	compute_patches();
	compute_normals();
	fix_edges();
	glutInitDisplayMode( GLUT_RGB | GLUT_DOUBLE | GLUT_DEPTH );
	glutInitWindowPosition( 50, 50 );
	glutInitWindowSize( 400, 300 );
	main_window_id = glutCreateWindow( "The Utah Teapot" );
	glutDisplayFunc( Display );
	glutMotionFunc( Motion );
	glutKeyboardFunc( Keyboard );
	glutReshapeFunc( Reshape );

	glEnable(GL_LIGHTING);
	glEnable( GL_NORMALIZE );
	glEnable(GL_LIGHT0);
	glLightfv(GL_LIGHT0, GL_AMBIENT, light0_ambient);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, light0_diffuse);
	glLightfv(GL_LIGHT0, GL_POSITION, light_specular);
	glLightfv(GL_LIGHT0, GL_POSITION, light0_position);
	glEnable(GL_LIGHT1);
	glLightfv(GL_LIGHT1, GL_AMBIENT, light1_ambient);
	glLightfv(GL_LIGHT1, GL_DIFFUSE, light1_diffuse);
	glLightfv(GL_LIGHT1, GL_POSITION, light_specular);
	glLightfv(GL_LIGHT1, GL_POSITION, light1_position);
	glEnable(GL_LIGHT2);
	glLightfv(GL_LIGHT2, GL_AMBIENT, light2_ambient);
	glLightfv(GL_LIGHT2, GL_DIFFUSE, light2_diffuse);
	glLightfv(GL_LIGHT2, GL_POSITION, light_specular);
	glLightfv(GL_LIGHT2, GL_POSITION, light2_position);

	glMaterialfv(GL_FRONT, GL_AMBIENT,  &material[0]);
	glMaterialfv(GL_FRONT, GL_DIFFUSE,  &material[1]);
	glMaterialfv(GL_FRONT, GL_SPECULAR, &material[2]);
	glMaterialf(GL_FRONT, GL_SHININESS, material[3]);
	glEnable(GL_DEPTH_TEST);
	glutMainLoop();
}