/*
	This program draws Bˇzier surfaces using evaluators, with an emphasis
	on extending control points from a single patch to a draw a full smooth
	surface.  This is adapted from work of Ben Eadington, CSU Stanislaus
*/

#include "glut.h"
#include <stdlib.h>
#include <stdio.h>
#include <math.h>

#define ROT_ANG 5    // angle of rotation for each keystroke
#define GRIDSIZE 16  // control point array is of size GRIDSIZE x GRIDSIZE
#define MAXHITS 200  // size of array to hold hit records for selection

typedef GLfloat point3[3];
typedef GLfloat point4[4];

static GLfloat Angle = 0.0;
static int RotAxis = -1;  // axis model for rotation

// initial size of window
GLint windW = 500, windH = 500;
GLint vp[4];

point4 white={1.0, 1.0, 1.0, 1.0};
point4 red={1.0, 0.0, 0.0, 1.0};
point4 green={0.0, 1.0, 0.0, 1.0};
point4 light_pos={30.0, 30.0, 30.0, 1.0};
GLfloat shiny[] = {25};

static GLfloat RotMatrix[16] =
 {1.0, 0.0, 0.0, 0.0,
  0.0, 1.0, 0.0, 0.0,
  0.0, 0.0, 1.0, 0.0,
  0.0, 0.0, 0.0, 1.0};

// control point array for pool surface
point3 ctrlpts[GRIDSIZE][GRIDSIZE];

/* prototypes for functions */
GLint DoSelect(GLint x, GLint y);
void myinit(void);
void kybd(unsigned char, int, int);
void genPoints(void);
void drawpoints(void);
void drawpatch(point3 patch[4][4]);
void surface(point3 ctrlpts[GRIDSIZE][GRIDSIZE]);
void render(void);
void display( void );
void reshape(int,int);

void myinit(void)
{
	glClearColor(0.0, 0.0, 0.0, 1.0);
	glLightfv(GL_LIGHT0, GL_AMBIENT, white);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, white);
	glLightfv(GL_LIGHT0, GL_SPECULAR, white);
	glLightfv(GL_LIGHT0, GL_POSITION, light_pos);

	glShadeModel(GL_SMOOTH);
	glEnable(GL_LIGHTING);
	glEnable(GL_LIGHT0);
	glEnable(GL_DEPTH_TEST);
	glEnable(GL_BLEND);
	glBlendFunc(GL_SRC_ALPHA, GL_ONE_MINUS_SRC_ALPHA);
	glEnable(GL_AUTO_NORMAL); // automatic generation of surface normals
	glEnable(GL_MAP2_VERTEX_3);

	genPoints();
}

void kybd(unsigned char key, int x, int y)
{
	switch(key) {
		case 'z': {
			Angle = ROT_ANG; RotAxis=0; break; }
		case 'x': {
			Angle = -ROT_ANG; RotAxis=0; break; }
		case 'a': {
			Angle = ROT_ANG; RotAxis=1; break; }
		case 's': {
			Angle = -ROT_ANG; RotAxis=1; break; }
		case 'q': {
			Angle = ROT_ANG; RotAxis=2; break; }
		case 'w': {
			Angle = -ROT_ANG; RotAxis=2; break; }
		default: {return;}
	}
	glutPostRedisplay();
}

void genPoints(void)
{
#define PI 3.14159
#define R1 6.0
#define R2 3.0
	int i, j;
	GLfloat alpha, beta, step;
	
	alpha = -PI;
	step = PI/(GLfloat)(GRIDSIZE-1);
	for (i=0; i<GRIDSIZE; i++) {
		beta = -PI;
		for (j=0; j<GRIDSIZE; j++) {
			ctrlpts[i][j][0] = (R1 + R2*cos(beta))*cos(alpha);
			ctrlpts[i][j][1] = (R1 + R2*cos(beta))*sin(alpha);
			ctrlpts[i][j][2] = R2*sin(beta);
			beta -= step;
		}
		alpha += step;
	}
}

void drawpoints(void)
{
	int i, j;
	int name=0;
	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, green);
	// iterate through control point array
	for(i=0; i<GRIDSIZE; i++)
		for(j=0; j<GRIDSIZE; j++) {
			glPushMatrix();
			// move point to appropriate location
			glTranslatef(ctrlpts[i][j][0],ctrlpts[i][j][1],ctrlpts[i][j][2]);
			// scale it to match scale of pool ( so it is still a cube )
			glScalef(1.0, 1.5, 1.5);
			glutSolidCube(0.10);
			glPopMatrix();
		}
}

// this function draws a patch defined by a 4 x 4 array
void drawpatch(point3 patch[4][4])
{
#define NUM 10

	glMap2f(GL_MAP2_VERTEX_3, 0.0, 1.0, 3, 4, 0.0, 1.0, 12, 4, &patch[0][0][0]);
	glMapGrid2f(NUM, 0.0, 1.0, NUM, 0.0, 1.0);
	glEvalMesh2(GL_FILL, 0, NUM, 0, NUM);
}

/*
	This function draws the entire surface defined by the control point
	array. Special cases must be defined for each corner and each edge, 
	for a total of 9 cases in all.  In each case the appropriate midpoints 
	are calculated, and the points are copied to a 4 x 4 array to be drawn.
*/
void surface(point3 ctrlpts[GRIDSIZE][GRIDSIZE])
{
#define LAST_STEP (GRIDSIZE/2)-1

	point3 patch[4][4];
	int xstep, ystep, i, j, k;

	for(xstep=0; xstep<LAST_STEP; xstep++)
		for(ystep=0; ystep<LAST_STEP; ystep++) {
			if(xstep==0 && ystep==0) { // front left corner of grid
				for(i=0; i<3; i++)
					for(j=0; j<3; j++)
						for(k=0; k<3; k++) patch[i][j][k] = ctrlpts[i][j][k];
				for(i=0; i<3; i++)
					for(k=0; k<3; k++)  {
						patch[i][3][k]=(ctrlpts[i][2][k]+ctrlpts[i][3][k])/2.0;
						patch[3][i][k]=(ctrlpts[2][i][k]+ctrlpts[3][i][k])/2.0;
					}
				for(k=0; k<3; k++)
					patch[3][3][k]=(ctrlpts[2][2][k]+ctrlpts[2][3][k]
										+ctrlpts[3][2][k]+ctrlpts[3][3][k])/4.0;
			}
			else if(xstep==0 && ystep==LAST_STEP-1) { // back left corner of grid
				for(i=0; i<3; i++)
					for(j=1; j<4; j++)
						for(k=0; k<3; k++) patch[i][j][k]=ctrlpts[i][GRIDSIZE-4+j][k];
				for(i=0; i<3; i++)
					for(k=0; k<3; k++)
						patch[i][0][k]=(ctrlpts[i][GRIDSIZE-4][k]
											+ctrlpts[i][GRIDSIZE-3][k])/2.0;
				for(i=1; i<4; i++)
					for(k=0; k<3; k++)
						patch[3][i][k]=(ctrlpts[2][GRIDSIZE-4+i][k]
											+ctrlpts[3][GRIDSIZE-4+i][k])/2.0;
				for(k=0; k<3; k++)
					patch[3][0][k]=(ctrlpts[2][GRIDSIZE-4][k]
										+ctrlpts[2][GRIDSIZE-3][k]
										+ctrlpts[3][GRIDSIZE-4][k]
										+ctrlpts[3][GRIDSIZE-3][k])/4.0;
			}
			else if(xstep==0) { /* patch is along left edge of grid */
				for(i=0; i<3; i++)
					for(j=1; j<3; j++)
						for(k=0; k<3; k++) patch[i][j][k]=ctrlpts[i][2*ystep+j][k];
				for(i=0; i<3; i++)
					for(k=0; k<3; k++) {
						patch[i][0][k]=(ctrlpts[i][2*ystep][k]
											+ctrlpts[i][2*ystep+1][k])/2.0;
						patch[i][3][k]=(ctrlpts[i][2*ystep+2][k]
											+ctrlpts[i][2*ystep+3][k])/2.0;
					}
				for(k=0; k<3; k++) {
					patch[3][0][k]=(ctrlpts[2][2*ystep][k]+ctrlpts[2][2*ystep+1][k]
										+ctrlpts[3][2*ystep][k]
										+ctrlpts[3][2*ystep+1][k])/4.0;
					patch[3][1][k]=(ctrlpts[2][2*ystep+1][k]
										+ctrlpts[3][2*ystep+1][k])/2.0;
					patch[3][2][k]=(ctrlpts[2][2*ystep+2][k]
										+ctrlpts[3][2*ystep+2][k])/2.0;
					patch[3][3][k]=(ctrlpts[2][2*ystep+2][k]+ctrlpts[2][2*ystep+3][k]
										+ctrlpts[3][2*ystep+2][k]
										+ctrlpts[3][2*ystep+3][k])/4.0;
				}
			}
			else if(xstep==LAST_STEP-1 && ystep==0) { // front right corner
				for(i=1; i<4; i++)
					for(j=0; j<3; j++)
						for(k=0; k<3; k++) patch[i][j][k]=ctrlpts[GRIDSIZE-4+i][j][k];
				for(j=0; j<3; j++)
					for(k=0; k<3; k++)
						patch[0][j][k]=(ctrlpts[GRIDSIZE-4][j][k]
											+ctrlpts[GRIDSIZE-3][j][k])/2.0;
				for(i=1; i<4; i++)
					for(k=0; k<3; k++)
						patch[i][3][k]=(ctrlpts[GRIDSIZE-4+i][3][k]
											+ctrlpts[GRIDSIZE-4+i][4][k])/2.0;
				for(k=0; k<3; k++)
					patch[0][3][k]=(ctrlpts[GRIDSIZE-4][2][k]
										+ctrlpts[GRIDSIZE-3][2][k]
										+ctrlpts[GRIDSIZE-4][3][k]
										+ctrlpts[GRIDSIZE-3][3][k])/4.0;
			}
			else if(ystep==0) { // patch on front edge of grid
				for(i=1; i<3; i++)
					for(j=0; j<3; j++)
						for(k=0; k<3; k++) patch[i][j][k]=ctrlpts[2*xstep+i][j][k];
				for(j=0; j<3; j++)
					for(k=0; k<3; k++) {
						patch[0][j][k]=(ctrlpts[2*xstep][j][k]
											+ctrlpts[2*xstep+1][j][k])/2.0;
						patch[3][j][k]=(ctrlpts[2*xstep+2][j][k]
											+ctrlpts[2*xstep+3][j][k])/2.0;
					}
				for(k=0; k<3; k++) {
					patch[0][3][k]=(ctrlpts[2*xstep][2][k]+ctrlpts[2*xstep+1][2][k]
										+ctrlpts[2*xstep][3][k]
										+ctrlpts[2*xstep+1][3][k])/4.0;
					patch[1][3][k]=(ctrlpts[2*xstep+1][2][k]
										+ctrlpts[2*xstep+1][3][k])/2.0;
					patch[2][3][k]=(ctrlpts[2*xstep+2][2][k]
										+ctrlpts[2*xstep+2][3][k])/2.0;
					patch[3][3][k]=(ctrlpts[2*xstep+2][2][k]+ctrlpts[2*xstep+3][2][k]
										+ctrlpts[2*xstep+2][3][k]
										+ctrlpts[2*xstep+3][3][k])/4.0;
				}
			}
			else if(xstep==LAST_STEP-1&&ystep==LAST_STEP-1) { //back right corner
				for(i=1; i<4; i++)
					for(j=1; j<4; j++)
						for(k=0; k<3; k++)
							patch[i][j][k]=ctrlpts[GRIDSIZE-4+i][GRIDSIZE-4+j][k];
				for(i=1; i<4; i++)
					for(k=0; k<3; k++) {
						patch[i][0][k]=(ctrlpts[GRIDSIZE-4+i][GRIDSIZE-4][k]
											+ctrlpts[GRIDSIZE-4+i][GRIDSIZE-3][k])/2.0;
						patch[0][i][k]=(ctrlpts[GRIDSIZE-4][GRIDSIZE-4+i][k]
											+ctrlpts[GRIDSIZE-3][GRIDSIZE-4+i][k])/2.0;
					}
				for(k=0; k<3; k++)
					patch[0][0][k]=(ctrlpts[GRIDSIZE-4][GRIDSIZE-4][k]
										+ctrlpts[GRIDSIZE-3][GRIDSIZE-4][k]
										+ctrlpts[GRIDSIZE-4][GRIDSIZE-3][k]
										+ctrlpts[GRIDSIZE-3][GRIDSIZE-3][k])/4.0;
			}
			else if(xstep==LAST_STEP-1) { // patch on right edge of grid
				for(i=1; i<4; i++)
					for(j=1; j<3; j++)
						for(k=0; k<3; k++)
							patch[i][j][k]=ctrlpts[GRIDSIZE-4+i][2*ystep+j][k];
				for(i=1; i<4; i++)
					for(k=0; k<3; k++) {
						patch[i][0][k]=(ctrlpts[GRIDSIZE-4+i][2*ystep][k]
											+ctrlpts[GRIDSIZE-4+i][2*ystep+1][k])/2.0;
						patch[i][3][k]=(ctrlpts[GRIDSIZE-4+i][2*ystep+2][k]
											+ctrlpts[GRIDSIZE-4+i][2*ystep+3][k])/2.0;
					}
				for(k=0; k<3; k++) {
					patch[0][0][k]=(ctrlpts[GRIDSIZE-4][2*ystep][k]
										+ctrlpts[GRIDSIZE-4][2*ystep+1][k]
										+ctrlpts[GRIDSIZE-3][2*ystep][k]
										+ctrlpts[GRIDSIZE-3][2*ystep+1][k])/4.0;
					patch[0][1][k]=(ctrlpts[GRIDSIZE-4][2*ystep+1][k]
										+ctrlpts[GRIDSIZE-3][2*ystep+1][k])/2.0;
					patch[0][2][k]=(ctrlpts[GRIDSIZE-4][2*ystep+2][k]
										+ctrlpts[GRIDSIZE-3][2*ystep+2][k])/2.0;
					patch[0][3][k]=(ctrlpts[GRIDSIZE-4][2*ystep+2][k]
										+ctrlpts[GRIDSIZE-4][2*ystep+3][k]
										+ctrlpts[GRIDSIZE-3][2*ystep+2][k]
										+ctrlpts[GRIDSIZE-3][2*ystep+3][k])/4.0;
				}
			}
			else if(ystep==LAST_STEP-1) { // patch on back edge of grid
				for(i=1; i<3; i++)
					for(j=1; j<4; j++)
						for(k=0; k<3; k++)
							patch[i][j][k]=ctrlpts[2*xstep+i][GRIDSIZE-4+j][k];
					for(j=1; j<4; j++)
						for(k=0; k<3; k++) {
							patch[0][j][k]=(ctrlpts[2*xstep][GRIDSIZE-4+j][k]
												+ctrlpts[2*xstep+1][GRIDSIZE-4+j][k])/2.0;
							patch[3][j][k]=(ctrlpts[2*xstep+2][GRIDSIZE-4+j][k]
												+ctrlpts[2*xstep+3][GRIDSIZE-4+j][k])/2.0;
					}
				for(k=0; k<3; k++) {
					patch[0][0][k]=(ctrlpts[2*xstep][GRIDSIZE-4][k]
										+ctrlpts[2*xstep+1][GRIDSIZE-4][k]
										+ctrlpts[2*xstep][GRIDSIZE-3][k]
										+ctrlpts[2*xstep+1][GRIDSIZE-3][k])/4.0;
					patch[1][0][k]=(ctrlpts[2*xstep+1][GRIDSIZE-4][k]
										+ctrlpts[2*xstep+1][GRIDSIZE-3][k])/2.0;
					patch[2][0][k]=(ctrlpts[2*xstep+2][GRIDSIZE-4][k]
										+ctrlpts[2*xstep+2][GRIDSIZE-3][k])/2.0;
					patch[3][0][k]=(ctrlpts[2*xstep+2][GRIDSIZE-4][k]
										+ctrlpts[2*xstep+3][GRIDSIZE-4][k]
										+ctrlpts[2*xstep+2][GRIDSIZE-3][k]
										+ctrlpts[2*xstep+3][GRIDSIZE-3][k])/4.0;
				}
			}
			else { // general case (internal patch)
				for(i=1; i<3; i++)
					for(j=1; j<3; j++)
						for(k=0; k<3; k++)
							patch[i][j][k]=ctrlpts[2*xstep+i][2*ystep+j][k];
				for(i=1; i<3; i++)
					for(k=0; k<3; k++) {
						patch[i][0][k]=(ctrlpts[2*xstep+i][2*ystep][k]
											+ctrlpts[2*xstep+i][2*ystep+1][k])/2.0;
						patch[i][3][k]=(ctrlpts[2*xstep+i][2*ystep+2][k]
											+ctrlpts[2*xstep+i][2*ystep+3][k])/2.0;
						patch[0][i][k]=(ctrlpts[2*xstep][2*ystep+i][k]
											+ctrlpts[2*xstep+1][2*ystep+i][k])/2.0;
						patch[3][i][k]=(ctrlpts[2*xstep+2][2*ystep+i][k]
											+ctrlpts[2*xstep+3][2*ystep+i][k])/2.0;
					}
				for(k=0; k<3; k++) {
					patch[0][0][k]=(ctrlpts[2*xstep][2*ystep][k]
										+ctrlpts[2*xstep+1][2*ystep][k]
										+ctrlpts[2*xstep][2*ystep+1][k]
										+ctrlpts[2*xstep+1][2*ystep+1][k])/4.0;
					patch[3][0][k]=(ctrlpts[2*xstep+2][2*ystep][k]
										+ctrlpts[2*xstep+3][2*ystep][k]
										+ctrlpts[2*xstep+2][2*ystep+1][k]
										+ctrlpts[2*xstep+3][2*ystep+1][k])/4.0;
					patch[0][3][k]=(ctrlpts[2*xstep][2*ystep+2][k]
										+ctrlpts[2*xstep+1][2*ystep+2][k]
										+ctrlpts[2*xstep][2*ystep+3][k]
										+ctrlpts[2*xstep+1][2*ystep+3][k])/4.0;
					patch[3][3][k]=(ctrlpts[2*xstep+2][2*ystep+2][k]
										+ctrlpts[2*xstep+3][2*ystep+2][k]
										+ctrlpts[2*xstep+2][2*ystep+3][k]
										+ctrlpts[2*xstep+3][2*ystep+3][k])/4.0;
				}
			}
		drawpatch(patch);  // draw the currently calculated patch
		}
}

void render(void)
{
	// perform appropriate rotations
	if (RotAxis==0) { // rotation about x-axis
		glPushMatrix();
		glLoadIdentity();
		glRotatef(Angle, 1.0, 0.0, 0.0);
		glMultMatrixf(RotMatrix);
		glGetFloatv(GL_MODELVIEW_MATRIX, RotMatrix);
		glPopMatrix();
	}
	else if (RotAxis==1) { // rotation about y-axis
		glPushMatrix();
		glLoadIdentity();
		glRotatef(Angle, 0.0, 1.0, 0.0);
		glMultMatrixf(RotMatrix);
		glGetFloatv(GL_MODELVIEW_MATRIX, RotMatrix);
		glPopMatrix();
	}
	else if (RotAxis==2) { // rotation about z-axis
		glPushMatrix();
		glLoadIdentity();
		glRotatef(Angle, 0.0, 0.0, 1.0);
		glMultMatrixf(RotMatrix);
		glGetFloatv(GL_MODELVIEW_MATRIX, RotMatrix);
		glPopMatrix();
	}
	glPushMatrix();
	// multiply modelview matrix (of the world view) by the rotation matrix
	glMultMatrixf(RotMatrix);
	glPushMatrix();
	glScalef(1.5,1.0,1.0);
	drawpoints();
	glMaterialfv(GL_FRONT_AND_BACK, GL_AMBIENT_AND_DIFFUSE, red);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SPECULAR, white);
	glMaterialfv(GL_FRONT_AND_BACK, GL_SHININESS, shiny);
	surface(ctrlpts);
	glPopMatrix();
	glPopMatrix();
	RotAxis=-1;  // don't rotate on next display
}

void display( void )
{
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	render();
	glutSwapBuffers();
}

void reshape(int w,int h)
{
	windW = w;
	windH = h;

	glViewport(0,0,(GLsizei)w,(GLsizei)h);
	glGetIntegerv(GL_VIEWPORT, vp);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(100.0,1.0,1.0,300);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
	gluLookAt(20.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0, 0.0);
}

void main(int argc, char** argv)
{
	glutInit(&argc,argv);
	glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glutInitWindowSize(500,500);
	glutInitWindowPosition(50, 50);
	glutCreateWindow("Spline surface, procedural control points");
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(kybd);

	myinit();
	glutMainLoop();
}