//This code was created by Jordan Isaak
//Code written in April, 2002
//My website is at members.shaw.ca/jordanisaak
//My email is jordanisaak@hotmail.com

#include <GL/glut.h>
#include <math.h>

#define PI 3.1415926535897932384626433832795
#define SPHERE_RESOLUTION 32

//Structure to hold all information
//for a section of a sphere
struct SPHERE_SECTION
{
	float boundingsphere[4];
	float aabb[24];
	float groupnormal[3];
	float tolerance;
	float radianstolerance;
	float points[3 * SPHERE_RESOLUTION * SPHERE_RESOLUTION];
};

/*Lighting variables*/
float lightpos[4] = {1.0f, 1.0f, 1.0f, 0.0f};
float lightambient[4] = {0.2f, 0.2f, 0.2f, 1.0f};
float lightdiffuse[4] = {1.0f, 1.0f, 1.0f, 1.0f};

float cameradistance = 5.0f;
float cameraelevation = 1.0f;
float camerarotation = 0.0f;
float camerapos[3] = {0.0f, 0.0f, 0.0f};
int lastx, lasty;

//Whether to use a bounding box method for group
//backface culling or to use a bounding sphere method
int cullmethod = 0;

//Whether or not to draw the sphere in wireframe
int wireframedraw = 0;

//6 sections to make up the sphere
SPHERE_SECTION sphere[6];

//Generates a section of a sphere
void generatespheresection(SPHERE_SECTION *section, int direction)
{
	int i, j;
	float veclen;
	float minx, miny, minz;
	float maxx, maxy, maxz;
	float temptolerance;
	float radius;

	//Generate the points that are on the sphere section
	if(direction == 0)
	{
		for(i = 0; i < SPHERE_RESOLUTION; i++)
		{
			for(j = 0; j < SPHERE_RESOLUTION; j++)
			{
				section->points[3 * (i * SPHERE_RESOLUTION + j)] = -0.5f + (float) j / (SPHERE_RESOLUTION - 1);
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 1] = -0.5f + (float) i / (SPHERE_RESOLUTION - 1);
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 2] = 0.5f;
			}
		}
	}
	else if(direction == 1)
	{
		for(i = 0; i < SPHERE_RESOLUTION; i++)
		{
			for(j = 0; j < SPHERE_RESOLUTION; j++)
			{
				section->points[3 * (i * SPHERE_RESOLUTION + j)] = 0.5f - (float) j / (SPHERE_RESOLUTION - 1);
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 1] = -0.5f + (float) i / (SPHERE_RESOLUTION - 1);
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 2] = -0.5f;
			}
		}
	}
	else if(direction == 2)
	{
		for(i = 0; i < SPHERE_RESOLUTION; i++)
		{
			for(j = 0; j < SPHERE_RESOLUTION; j++)
			{
				section->points[3 * (i * SPHERE_RESOLUTION + j)] = -0.5;
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 1] = -0.5f + (float) i / (SPHERE_RESOLUTION - 1);
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 2] = -0.5f + (float) j / (SPHERE_RESOLUTION - 1);
			}
		}
	}
	else if(direction == 3)
	{
		for(i = 0; i < SPHERE_RESOLUTION; i++)
		{
			for(j = 0; j < SPHERE_RESOLUTION; j++)
			{
				section->points[3 * (i * SPHERE_RESOLUTION + j)] = 0.5;
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 1] = 0.5f - (float) i / (SPHERE_RESOLUTION - 1);
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 2] = -0.5f + (float) j / (SPHERE_RESOLUTION - 1);
			}
		}
	}
	else if(direction == 4)
	{
		for(i = 0; i < SPHERE_RESOLUTION; i++)
		{
			for(j = 0; j < SPHERE_RESOLUTION; j++)
			{
				section->points[3 * (i * SPHERE_RESOLUTION + j)] = -0.5f + (float) i / (SPHERE_RESOLUTION - 1);
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 1] = -0.5;
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 2] = 0.5f - (float) j / (SPHERE_RESOLUTION - 1);
			}
		}
	}
	else if(direction == 5)
	{
		for(i = 0; i < SPHERE_RESOLUTION; i++)
		{
			for(j = 0; j < SPHERE_RESOLUTION; j++)
			{
				section->points[3 * (i * SPHERE_RESOLUTION + j)] = -0.5f + (float) i / (SPHERE_RESOLUTION - 1);
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 1] = 0.5;
				section->points[3 * (i * SPHERE_RESOLUTION + j) + 2] = -0.5f + (float) j / (SPHERE_RESOLUTION - 1);
			}
		}
	}

	//Normalize points to lie on a unit sphere
	for(i = 0; i < SPHERE_RESOLUTION * SPHERE_RESOLUTION; i++)
	{
		veclen = (float) sqrt(section->points[3 * i] * section->points[3 * i]
							+ section->points[3 * i + 1] * section->points[3 * i + 1]
							+ section->points[3 * i + 2] * section->points[3 * i + 2]);
		veclen = 1.0f / veclen;
		section->points[3 * i] *= veclen;
		section->points[3 * i + 1] *= veclen;
		section->points[3 * i + 2] *= veclen;
	}

	//Calculate an aabb for the sphere section
	minx = section->points[0];
	maxx = minx;
	miny = section->points[1];
	maxy = miny;
	minz = section->points[2];
	maxz = minz;
	for(i = 1; i < SPHERE_RESOLUTION * SPHERE_RESOLUTION; i++)
	{
		if(section->points[3 * i] < minx)
			minx = section->points[3 * i];
		else if(section->points[3 * i] > maxx)
			maxx = section->points[3 * i];

		if(section->points[3 * i + 1] < miny)
			miny = section->points[3 * i + 1];
		else if(section->points[3 * i + 1] > maxy)
			maxy = section->points[3 * i + 1];

		if(section->points[3 * i + 2] < minz)
			minz = section->points[3 * i + 2];
		else if(section->points[3 * i + 2] > maxz)
			maxz = section->points[3 * i + 2];
	}
	section->aabb[0] = minx;
	section->aabb[1] = miny;
	section->aabb[2] = minz;
	section->aabb[3] = maxx;
	section->aabb[4] = miny;
	section->aabb[5] = minz;
	section->aabb[6] = minx;
	section->aabb[7] = maxy;
	section->aabb[8] = minz;
	section->aabb[9] = maxx;
	section->aabb[10] = maxy;
	section->aabb[11] = minz;
	section->aabb[12] = minx;
	section->aabb[13] = miny;
	section->aabb[14] = minz;
	section->aabb[15] = maxx;
	section->aabb[16] = miny;
	section->aabb[17] = maxz;
	section->aabb[18] = minx;
	section->aabb[19] = maxy;
	section->aabb[20] = maxz;
	section->aabb[21] = maxx;
	section->aabb[22] = maxy;
	section->aabb[23] = maxz;

	//Calculate a normal and tolerance to use for back patch culling
	section->groupnormal[0] = 0.0f;
	section->groupnormal[1] = 0.0f;
	section->groupnormal[2] = 0.0f;
	for(i = 0; i < SPHERE_RESOLUTION * SPHERE_RESOLUTION; i++)
	{
		section->groupnormal[0] += section->points[3 * i];
		section->groupnormal[1] += section->points[3 * i + 1];
		section->groupnormal[2] += section->points[3 * i + 2];
	}

	//Calculate bounding sphere
	section->boundingsphere[0] = section->groupnormal[0] / (SPHERE_RESOLUTION * SPHERE_RESOLUTION);
	section->boundingsphere[1] = section->groupnormal[1] / (SPHERE_RESOLUTION * SPHERE_RESOLUTION);
	section->boundingsphere[2] = section->groupnormal[2] / (SPHERE_RESOLUTION * SPHERE_RESOLUTION);
	section->boundingsphere[3] = 0.0f;
	for(i = 0; i < SPHERE_RESOLUTION * SPHERE_RESOLUTION; i++)
	{
		radius = (float) sqrt((section->boundingsphere[0] - section->points[3 * i]) * (section->boundingsphere[0] - section->points[3 * i])
				+ (section->boundingsphere[1] - section->points[3 * i + 1]) * (section->boundingsphere[1] - section->points[3 * i + 1])
				+ (section->boundingsphere[2] - section->points[3 * i + 2]) * (section->boundingsphere[2] - section->points[3 * i + 2]));
		if(radius < section->boundingsphere[3])
			section->boundingsphere[3] = radius;
	}

	//Normalize group normal of the triangles
	veclen = (float) sqrt(section->groupnormal[0] * section->groupnormal[0]
							+ section->groupnormal[1] * section->groupnormal[1]
							+ section->groupnormal[2] * section->groupnormal[2]);
	veclen = 1.0f / veclen;
	section->groupnormal[0] *= veclen;
	section->groupnormal[1] *= veclen;
	section->groupnormal[2] *= veclen;

	//First of all, find the normal from the group of triangles
	//that differs the most from average group normal, and store
	//how much it differs
	section->tolerance = 1.0f;
	for(i = 0; i < SPHERE_RESOLUTION * SPHERE_RESOLUTION; i++)
	{
		temptolerance = section->groupnormal[0] * section->points[3 * i]
				+ section->groupnormal[1] * section->points[3 * i + 1]
				+ section->groupnormal[2] * section->points[3 * i + 2];
		if(temptolerance < section->tolerance)
			section->tolerance = temptolerance;
	}

	//Set the tolerance value to be equal to the cosine of the arcsine of the current value
	//This value indicates the minimum that the dot product from
	//the viewpoint to a point on the convex hull of the group of triangles
	//and the group normal can be without the group of triangles possibly
	//being visible
	section->tolerance = (float)asin(section->tolerance);
	section->radianstolerance = section->tolerance;
	section->tolerance = (float)cos(section->tolerance);

	//Square this value so that square roots are not needed for each calculation
	section->tolerance = section->tolerance * section->tolerance;
}

void drawspheresection(SPHERE_SECTION *section)
{
	int i, j;
	float viewvector[3];
	int visible = 0;
	float a, b, veclength, tolerance;

	if(cullmethod == 0)
	{
		//Check to see if this section of the sphere is entirely back-facing
		//if it is, don't bother rendering it
		//Check visibility of each vertex of the bounding box to correctly
		//Determine if the group of triangles is back facing
		for(i = 0; i < 8; i++)
		{
			//Calculate vector from camera to vertex
			viewvector[0] = section->aabb[3 * i] - camerapos[0];
			viewvector[1] = section->aabb[3 * i + 1] - camerapos[1];
			viewvector[2] = section->aabb[3 * i + 2] - camerapos[2];

			//Calculate dot product of view vector and the
			//average normal of the group of faces
			a = viewvector[0] * section->groupnormal[0]
				+ viewvector[1] * section->groupnormal[1]
				+ viewvector[2] * section->groupnormal[2];

			//If the average normal is pointing towards the camera,
			//the group of triangles is visible, so the program doesn't
			//need to do any more calculations
			if(a < 0.0f)
			{
				visible = 1;
				break;
			}

			//Square this value
			a = a * a;

			//Calculate the length of the view vector squared
			b = viewvector[0] * viewvector[0]
				+ viewvector[1] * viewvector[1]
				+ viewvector[2] * viewvector[2];

			//Check to see if the group normal is back-facing enough
			//to have all possible normals be back facing
			//If not, then the group of triangles is visible, so no more
			//calculations are needed
			if(a / b < section->tolerance)
			{
				visible = 1;
				break;
			}
		}
	}
	else if(cullmethod == 1)
	{
		//Check to see if this section of the sphere is entirely back-facing
		//if it is, don't bother rendering it
		//Check visibility using the bounding sphere of the section
		//to determine the tolerance to use as a threshold

		//Calculate vector from camera to center of bounding sphere
		viewvector[0] = section->boundingsphere[0] - camerapos[0];
		viewvector[1] = section->boundingsphere[1] - camerapos[1];
		viewvector[2] = section->boundingsphere[2] - camerapos[2];

		//Calculate dot product of view vector and average group normal
		a = viewvector[0] * section->groupnormal[0]
			+ viewvector[1] * section->groupnormal[1]
			+ viewvector[2] * section->groupnormal[2];

		if(a > 0.0f)
		{
			//Calculate length of vector
			veclength = (float) sqrt(viewvector[0] * viewvector[0]
									+ viewvector[1] * viewvector[1]
									+ viewvector[2] * viewvector[2]);

			veclength = 1.0f / veclength;
		
			//Calculate tolerance for determining if the object is back facing or not
			tolerance = cos(section->radianstolerance - asin(section->boundingsphere[3] * veclength));

			//Normalize the dot product
			a *= veclength;

			//Check to see if the dot product is below the tolerance
			if(a < tolerance)
				visible = 1;
		}
		else
			visible = 1;
	}

	//If the result of the visibility calculations is
	//that the group of triangles is not visible, then
	//don't draw this group of triangles
	if(!visible)
		return;

	//Draw a series of triangle strips for the sphere section
	for(i = 0; i < SPHERE_RESOLUTION - 1; i++)
	{
		glBegin(GL_TRIANGLE_STRIP);
		for(j = 0; j < SPHERE_RESOLUTION; j++)
		{
			glTexCoord2f((float)j / (SPHERE_RESOLUTION - 1), (float)(i + 1) / (SPHERE_RESOLUTION - 1));
			glNormal3fv(&(section->points[3 * ((i + 1) * SPHERE_RESOLUTION + j)]));
			glVertex3fv(&(section->points[3 * ((i + 1) * SPHERE_RESOLUTION + j)]));
			glTexCoord2f((float)j / (SPHERE_RESOLUTION - 1), (float)i / (SPHERE_RESOLUTION - 1));
			glNormal3fv(&(section->points[3 * (i * SPHERE_RESOLUTION + j)]));
			glVertex3fv(&(section->points[3 * (i * SPHERE_RESOLUTION + j)]));
		}
		glEnd();
	}
}

void idle(void)
{
	glutPostRedisplay();
}

void init(void)
{
	int i;

	//General OpenGL setup
	glShadeModel(GL_SMOOTH);
	glEnable(GL_CULL_FACE);
	glEnable(GL_DEPTH_TEST);

	//Set up lighting
	glLightfv(GL_LIGHT0, GL_AMBIENT, lightambient);
	glLightfv(GL_LIGHT0, GL_DIFFUSE, lightdiffuse);
	glEnable(GL_LIGHTING);
	glEnable(GL_LIGHT0);
	glLightModeli(GL_LIGHT_MODEL_TWO_SIDE, 0);

	//Generate 6 patches to create a sphere
	for(i = 0; i < 6; i++)
		generatespheresection(&(sphere[i]), i);
}

void display(void)
{
	int i;

	//Set the camera position
	camerapos[0] = sin(camerarotation) * sin(cameraelevation) * cameradistance;
	camerapos[1] = cos(cameraelevation) * cameradistance;
	camerapos[2] = cos(camerarotation) * sin(cameraelevation) * cameradistance;

	//Set up for rendering
	glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	glLoadIdentity();
	gluLookAt(camerapos[0], camerapos[1], camerapos[2],
				0.0f, 0.0f, 0.0f,
				0.0f, 1.0f, 0.0f);

	//Set the light position
	glLightfv(GL_LIGHT0, GL_POSITION, lightpos);

	//Set appropriate draw settings
	if(wireframedraw)
	{
		glDisable(GL_LIGHTING);
		glDisable(GL_CULL_FACE);
		glDisable(GL_DEPTH_TEST);
		glPolygonMode(GL_FRONT_AND_BACK, GL_LINE);
	}
	else
	{
		glEnable(GL_LIGHTING);
		glEnable(GL_CULL_FACE);
		glEnable(GL_DEPTH_TEST);
		glPolygonMode(GL_FRONT_AND_BACK, GL_FILL);
	}

	//Draw each section of the sphere
	for(i = 0; i < 6; i++)
		drawspheresection(&(sphere[i]));

	glutSwapBuffers();
}

void reshape(int w, int h)
{
	glViewport(0, 0, (GLsizei) w, (GLsizei) h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	gluPerspective(60.0, (GLfloat) w/(GLfloat) h, 0.1, 100.0);
	glMatrixMode(GL_MODELVIEW);
}

void mousemove(int x, int y)
{
	camerarotation += (x - lastx) / 100.0f;
	cameraelevation += (y - lasty) / 100.0f;

	if(cameraelevation < 0.1f)
			cameraelevation = 0.1f;
	if(cameraelevation > 3.0f)
			cameraelevation = 3.0f;

	lastx = x;
	lasty = y;
	idle();
}

void mousedown(int button, int state, int x, int y)
{
	lastx = x;
	lasty = y;
}

void keyboard (unsigned char key, int x, int y)
{
	switch (key)
	{
		//Quit the program
		case 27:
			exit(0);
			break;
		//Bring the camera farther away from the sphere
		case 'z':
			cameradistance += 0.5f;
			break;
		//Bring the camera closer to the sphere
		case 'x':
			cameradistance -= 0.5f;
			if(cameradistance < 3.0f)
				cameradistance = 3.0f;
			break;
		//Toggle between wireframe and solid draw modes
		case 'w':
			wireframedraw = -wireframedraw + 1;
			break;
		case 'c':
			cullmethod = -cullmethod + 1;
			break;
		default:
			break;
	}
}

int main(int argc, char** argv)
{
	glutInit(&argc, argv);
	glutInitDisplayMode(GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glutInitWindowSize(800, 600);
	glutInitWindowPosition(100, 100);
	glutCreateWindow(argv[0]);
	init();

	//Set up control functions
	glutIdleFunc(idle);
	glutDisplayFunc(display);
	glutReshapeFunc(reshape);
	glutKeyboardFunc(keyboard);
	glutMotionFunc(mousemove);
	glutMouseFunc(mousedown);
	glutMainLoop();
	return 0; 
}