#include <GL/glut.h>

#define RECURSTION_DEPTH 4

void myInit () {
	
	glClearColor (1.0, 1.0, 1.0, 1.0);
	glPointSize (1);
	
	glEnable (GL_DEPTH_TEST);
	
	glMatrixMode (GL_PROJECTION);
	glLoadIdentity();
	glOrtho (-5.0, 5.0, 0.0, 5.0, -5.0, 5.0);
	glMatrixMode (GL_MODELVIEW);
}

void sierpinski3D (float vertex1[], float vertex2[], float vertex3[], float vertex4[], int current_depth) {
	
	float new_vertex1[3], new_vertex2[3], new_vertex3[3], new_vertex4[3];
	float new_vertex5[3], new_vertex6[3];
	
	if (current_depth > RECURSTION_DEPTH)
		return;
	
	/* Plot the triangles only in the deeptest of the depths ;-) */	
	if (current_depth == RECURSTION_DEPTH) {
		glBegin (GL_TRIANGLES);
			glColor3f (1.0, 0.0, 0.0);
			glVertex3fv (vertex1);
			glVertex3fv (vertex2);
			glVertex3fv (vertex3);
	
			glColor3f (0.0, 1.0, 0.0);
			glVertex3fv (vertex1);
			glVertex3fv (vertex2);
			glVertex3fv (vertex4);
			
			glColor3f (0.0, 0.0, 1.0);
			glVertex3fv (vertex1);
			glVertex3fv (vertex3);
			glVertex3fv (vertex4);
			
			glColor3f (0.0, 0.0, 0.0);
			glVertex3fv (vertex2);
			glVertex3fv (vertex3);
			glVertex3fv (vertex4);
		glEnd();
	}

	
	/* The tetrahedron has 6 edges, find the 6 midpoints of each of them */
	
	/* first mid-point b/w vertices 1 and 2 */
	new_vertex1[0] = (vertex1[0]+vertex2[0])/2;
	new_vertex1[1] = (vertex1[1]+vertex2[1])/2;
	new_vertex1[2] = (vertex1[2]+vertex2[2])/2;
	
	/* second mid-point b/w vertices 1 and 3 */
	new_vertex2[0] = (vertex1[0]+vertex3[0])/2;
	new_vertex2[1] = (vertex1[1]+vertex3[1])/2;
	new_vertex2[2] = (vertex1[2]+vertex3[2])/2;
	
	/* third mid-point b/w vertices 1 and 4 */
	new_vertex3[0] = (vertex1[0]+vertex4[0])/2;
	new_vertex3[1] = (vertex1[1]+vertex4[1])/2;
	new_vertex3[2] = (vertex1[2]+vertex4[2])/2;
	
	/* fourth mid-point b/w vertices 2 and 3 */
	new_vertex4[0] = (vertex2[0]+vertex3[0])/2;
	new_vertex4[1] = (vertex2[1]+vertex3[1])/2;
	new_vertex4[2] = (vertex2[2]+vertex3[2])/2;
	
	/* fifth mid-point b/w vertices 2 and 4 */
	new_vertex5[0] = (vertex2[0]+vertex4[0])/2;
	new_vertex5[1] = (vertex2[1]+vertex4[1])/2;
	new_vertex5[2] = (vertex2[2]+vertex4[2])/2;
	
	/* sixth mid-point b/w vertices 3 and 4 */
	new_vertex6[0] = (vertex3[0]+vertex4[0])/2;
	new_vertex6[1] = (vertex3[1]+vertex4[1])/2;
	new_vertex6[2] = (vertex3[2]+vertex4[2])/2;
	
	/* Now call the sierpinski3D() for each of the 4 sub-tetrahedrons formed.  
	Each tetrahedron consists of one of the initial four vertices and the
	mid-points of the three corresponding edges */

	sierpinski3D (vertex1, new_vertex1, new_vertex2, new_vertex3, current_depth+1);
	sierpinski3D (vertex2, new_vertex1, new_vertex4, new_vertex5, current_depth+1);
	sierpinski3D (vertex3, new_vertex2, new_vertex4, new_vertex6, current_depth+1);
	sierpinski3D (vertex4, new_vertex3, new_vertex5, new_vertex6, current_depth+1);
}

void display () {
	
	float vertex1[3] = {0.0, 2.5, 5.0};
	float vertex2[3] = {5.0, 0.0, 0.0};
	float vertex3[3] = {-5.0, 0.0, 0.0};
	float vertex4[3] = {0.0, 5.0, 0.0};
	
	glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	sierpinski3D (vertex1, vertex2, vertex3, vertex4, 0);
	glFlush();
}

int main (int argc, char **argv) {
	
	glutInit (&argc, argv);
	glutCreateWindow ("Sierpinski3D");
	glutInitDisplayMode (GLUT_RGB | GLUT_DEPTH | GLUT_SINGLE);
	glutInitWindowPosition (400, 400);
	glutInitWindowSize (800, 600);
	glutDisplayFunc (display);
	myInit();
	glutMainLoop ();
}