// How To Make sphere
// explanation of sphere mesh modeling
// Gouraud shading
// T. Kosaka CS TNCT 2001

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

#ifndef M_PI
#define M_PI 3.141592653589793
#endif

typedef struct {
	float m[4][4];
} matrix3D_t;

typedef struct {
	float v[4];
} vector3D_t;

typedef struct {
	float x;
	float y;
	float z;
} point3D_t;

typedef struct {
	float x;
	float y;
} point2D_t;

typedef struct {
	float r;
	float g;
	float b;
} color_t;

////////////////// matrices and vectors 3D ver 2 /////////////////
matrix3D_t createIdentity(void)
{
	matrix3D_t u;
	int i,j;
	for (i=0;i<4;i++) {
		for(j=0;j<4;j++) u.m[i][j]=0.;
		u.m[i][i]=1.;
	}
	return u;
}

matrix3D_t operator * (matrix3D_t a,matrix3D_t b)
{
	matrix3D_t c;//c=a*b
	int i,j,k;
	for (i=0;i<4;i++) for (j=0;j<4;j++) {
		c.m[i][j]=0;
		for (k=0;k<4;k++) c.m[i][j]+=a.m[i][k]*b.m[k][j];
	}
	return c;
}

vector3D_t operator * (matrix3D_t a, vector3D_t b)
{
	vector3D_t c;//c=a*b
	int i,j;
	for (i=0;i<4;i++) {
		c.v[i]=0;
		for (j=0;j<4;j++) c.v[i]+=a.m[i][j]*b.v[j];
	}
	return c;
}

matrix3D_t translationMTX(float dx,float dy,float dz)
{
	matrix3D_t trans=createIdentity();
	trans.m[0][3]=dx;
	trans.m[1][3]=dy;
	trans.m[2][3]=dz;
	return trans;
}

matrix3D_t rotationXMTX(float theta)
{
	matrix3D_t rotate=createIdentity();
	float cs=cos(theta);
	float sn=sin(theta);
	rotate.m[1][1]=cs; rotate.m[1][2]=-sn;
	rotate.m[2][1]=sn; rotate.m[2][2]=cs;
	return rotate;
}

matrix3D_t rotationYMTX(float theta)
{
	matrix3D_t rotate=createIdentity();
	float cs=cos(theta);
	float sn=sin(theta);
	rotate.m[0][0]=cs; rotate.m[0][2]=sn;
	rotate.m[2][0]=-sn; rotate.m[2][2]=cs;
	return rotate;
}

matrix3D_t rotationZMTX(float theta)
{
	matrix3D_t rotate=createIdentity();
	float cs=cos(theta);
	float sn=sin(theta);
	rotate.m[0][0]=cs; rotate.m[0][1]=-sn;
	rotate.m[1][0]=sn; rotate.m[1][1]=cs;
	return rotate;
}

matrix3D_t scalingMTX(float factorx,float factory,float factorz)
{
	matrix3D_t scale=createIdentity();
	scale.m[0][0]=factorx;
	scale.m[1][1]=factory;
	scale.m[2][2]=factorz;
	return scale;
}

matrix3D_t perspectiveMTX(float eyelength)
{
	matrix3D_t perspective=createIdentity();
	perspective.m[3][2]=-1./eyelength;
	return perspective;
}

point2D_t Vector2Point2D(vector3D_t vec)
{
	point2D_t pnt;
	pnt.x=vec.v[0];
	pnt.y=vec.v[1];
	return pnt;
}

point3D_t Vector2Point3D(vector3D_t vec)
{
	point3D_t pnt;
	pnt.x=vec.v[0];
	pnt.y=vec.v[1];
	pnt.z=vec.v[2];
	return pnt;
}

vector3D_t Point2Vector(point3D_t pnt)
{
	vector3D_t vec;
	vec.v[0]=pnt.x;
	vec.v[1]=pnt.y;
	vec.v[2]=pnt.z;
	vec.v[3]=1.;
	return vec;
}

vector3D_t homogenizeVector(vector3D_t vec)
{
	int i;
	for (i=0;i<3;i++) {
		vec.v[i]/=vec.v[3];
	}
	vec.v[3]=1.;
	return vec;
}

vector3D_t unitVector(vector3D_t vec)
{
	int i;
	float vec2=0.;
	float vec1,invvec1;
	for (i=0;i<3;i++) {
		vec2+=vec.v[i]*vec.v[i];
	}
	vec1=sqrt(vec2);
	if (vec1!=0.) {
		invvec1=1./vec1;
		for (i=0;i<3;i++) {
			vec.v[i]*=invvec1;
		}
	}
	vec.v[3]=1.;
	return vec;
}

// inner product (dot product) of homogeneous vector
float operator * (vector3D_t a, vector3D_t b)
{
	float c;//c=a*b
	int i;
	c=0;
	for (i=0;i<3;i++) {
		c+=a.v[i]*b.v[i];
	}
	return c;
}

// outer product (cross product ) of homogeneous vector
//       i         j         k
//       a0       a1        a2
//       b0       b1        b2
vector3D_t operator ^ (vector3D_t a, vector3D_t b)
{
	vector3D_t c;//c=a*b
	c.v[0]=a.v[1]*b.v[2]-a.v[2]*b.v[1];
	c.v[1]=a.v[2]*b.v[0]-a.v[0]*b.v[2];
	c.v[2]=a.v[0]*b.v[1]-a.v[1]*b.v[0];
	c.v[3]=1.;
	return c;
}

vector3D_t operator - (vector3D_t v1,vector3D_t v0)
{
	vector3D_t c;//c=v1-v0
	c.v[0]=v1.v[0]-v0.v[0];
	c.v[1]=v1.v[1]-v0.v[1];
	c.v[2]=v1.v[2]-v0.v[2];
	c.v[3]=1.;
	return c;
}

vector3D_t operator - (vector3D_t v)
{
	vector3D_t c;//c=-v
	c.v[0]=-v.v[0];
	c.v[1]=-v.v[1];
	c.v[2]=-v.v[2];
	c.v[3]=1.;
	return c;
}

vector3D_t operator * (float r, vector3D_t b)
{
	vector3D_t c;//c=r*b
	int i;
	for (i=0;i<3;i++) {
		c.v[i]=r*b.v[i];
	}
	c.v[3]=1.;
	return c;
}

vector3D_t operator * (vector3D_t b, float r)
{
	vector3D_t c;//c=r*b
	int i;
	for (i=0;i<3;i++) {
		c.v[i]=r*b.v[i];
	}
	c.v[3]=1.;
	return c;
}

float funcPositive(float x)
{
	if (0.<x) return x;
	else return 0.;
}

// x to yth power
float power(float x,float y)
{
	//ln z = y ln x        z = exp (y ln x)
	if (x==0.) return 0;
	return exp(y*log(x));
}

color_t operator + (color_t c1, color_t c2)
{
	color_t col;
	col.r=c1.r+c2.r;
	col.g=c1.g+c2.g;
	col.b=c1.b+c2.b;
	return col;
}

color_t operator * (float r, color_t c)
{
	color_t col;
	col.r=r*c.r;
	col.g=r*c.g;
	col.b=r*c.b;
	return col;
}

color_t operator * (color_t c, float r)
{
	color_t col;
	col.r=r*c.r;
	col.g=r*c.g;
	col.b=r*c.b;
	return col;
}

//PhongModel color calculation
// LightVector, NormalVector, ViewVector, ColorofObject
color_t PhongModel(vector3D_t Light,vector3D_t Normal,vector3D_t View,color_t col)
{
	float kspe=0.7; // specular reflection coefficient
	float kdif=0.6; // diffuse reflection coefficient
	float kamb=0.4; // ambient light coefficient
	float tmp,NL,RV;
	color_t ColWhite={1,1,1};
	vector3D_t ReflectionVector=(2.*(Light*Normal)*Normal)-Light;
	tmp=Normal*Light;
	NL=funcPositive(tmp);
	tmp=ReflectionVector*View;
	RV=funcPositive(tmp);
	return kdif*NL*col+kspe*power(RV,4)*ColWhite+kamb*col;
}

////////////// End of matrices and vectors 3D ver 2 //////////////

////////////// OpenGL drawShape Functions ver 1 /////////////////
void setColor(float red,float green,float blue)
{
	glColor3f(red, green, blue);	
}

void setColor(color_t col)
{
	glColor3f(col.r, col.g, col.b);	
}

void drawDot(float x,float y)
{
	glBegin(GL_POINTS);
		glVertex2f(x, y);
	glEnd();
}

void drawLine(float x1, float y1, float x2, float y2)
{
	glBegin(GL_LINES);
		glVertex2f(x1, y1);
		glVertex2f(x2, y2);
	glEnd();
}

void drawLine(point2D_t p1,point2D_t p2)
{
	drawLine(p1.x,p1.y,p2.x,p2.y);
}

//n: number of points
void drawPolyline(point2D_t pnt[],int n)
{
	int i;
	glBegin(GL_LINE_STRIP);
		for (i=0;i<n;i++) {
			glVertex2f(pnt[i].x, pnt[i].y);
		}
	glEnd();
}

//n: number of vertices
void drawPolygon(point2D_t pnt[],int n)
{
	int i;
	glBegin(GL_LINE_LOOP);
		for (i=0;i<n;i++) {
			glVertex2f(pnt[i].x, pnt[i].y);
		}
	glEnd();
}

// The function fillPolygon can fills only convex polygons
//n: number of vertices
void fillPolygon(point2D_t pnt[],int n,color_t color)
{
	int i;
	setColor(color);
	glBegin(GL_POLYGON);
		for (i=0;i<n;i++) {
			glVertex2f(pnt[i].x, pnt[i].y);
		}
	glEnd();
}

// The function gradatePolygon can fills only convex polygons
// The vertices will be painted with corresponding given colors.
// The points inside the polygon will be painted with the mixed color.
//n: number of vertices
void gradatePolygon(point2D_t pnt[],color_t col[],int num)
{
	int i;
	glBegin(GL_POLYGON);
		for (i=0;i<num;i++) {
			setColor(col[i]);
			glVertex2f(pnt[i].x, pnt[i].y);
		}
	glEnd();
}

//////////// End of OpenGL drawShape Functions ver 1 ////////////

void userdraw(void);

void display(void)
{
	glClear( GL_COLOR_BUFFER_BIT);
	userdraw();
	glutSwapBuffers();
}

//////////////////////////////////////////////////////////////////
void drawcharX(float x,float y)
{
	drawLine(x,y,x+10,y+12);drawLine(x,y+12,x+10,y);
}

void drawcharY(float x,float y)
{
	drawLine(x+5,y,x+5,y+7);drawLine(x,y+12,x+5,y+7);drawLine(x+10,y+12,x+5,y+7);
}

void drawcharZ(float x,float y)
{
	drawLine(x,y+12,x+10,y+12);drawLine(x+10,y+12,x,y);drawLine(x,y,x+10,y);
}

void drawAxes(matrix3D_t view)
{
#define HALFAXIS  220
#define HALFAXIS1 (HALFAXIS-10)
	point3D_t axes[14]={
		{-HALFAXIS,0,0},{HALFAXIS,0,0},{HALFAXIS1,5,0},{HALFAXIS1,0,0},{0,0,0},
		{0,-HALFAXIS,0},{0,HALFAXIS,0},{0,HALFAXIS1,5},{0,HALFAXIS1,0},{0,0,0},
		{0,0,-HALFAXIS},{0,0,HALFAXIS},{5,0,HALFAXIS1},{0,0,HALFAXIS1}
	};
	vector3D_t vec[14];
	point2D_t buff[14];
	int i;
	for (i=0;i<14;i++) {
		vec[i]=Point2Vector(axes[i]);
		vec[i]=view*vec[i];
		buff[i]=Vector2Point2D(vec[i]);
	}
	drawPolyline(buff,14);
	drawcharX(buff[1].x,buff[1].y);
	drawcharY(buff[6].x,buff[6].y);
	drawcharZ(buff[11].x-14,buff[11].y);
}

//////////////////////////////////////////////////////////////////
typedef struct {
	int NumberofVertices; //in the face
	short int pnt[16];
} face_t;
typedef struct {
	int NumberofVertices; //of the object
	point3D_t pnt[1024];
	vector3D_t NormalVector[1024]; // at the pnt
	int NumberofFaces; //of the object
	face_t fc[1024];
} object3D_t;

int NumLatEdge=16;
int NumLngEdge=32;

static void makeSpherePointData(float r,object3D_t &sphere)
{
	// r : radius of the sphere
	float latitude;
	float longitude;
	int i,j,k;
	vector3D_t unitz={0,0,1,1};
	vector3D_t unityn={0,1,0,1};
	vector3D_t unitys={0,-1,0,1};
	matrix3D_t mat;
	
	sphere.NumberofVertices=1+(NumLatEdge-1)*NumLngEdge+1;

	sphere.NormalVector[0]=unityn;
	sphere.pnt[0]=Vector2Point3D(r*sphere.NormalVector[0]);
	for(i=1;i<NumLatEdge;i++) {
		latitude=(NumLatEdge/2.-i)/(float)NumLatEdge*M_PI;
		mat=rotationXMTX(-latitude);
		for(j=0;j<NumLngEdge;j++) {
			k=NumLngEdge*(i-1)+j+1;
			longitude=(float)j/NumLngEdge*2.*M_PI;
			sphere.NormalVector[k]=rotationYMTX(longitude)*mat*unitz;
			sphere.pnt[k]=Vector2Point3D(r*sphere.NormalVector[k]);
		}
	}
	k=(NumLatEdge-1)*NumLngEdge+1;
	sphere.NormalVector[k]=unitys;
	sphere.pnt[k]=Vector2Point3D(r*sphere.NormalVector[k]);

	sphere.NumberofFaces=NumLatEdge*NumLngEdge;

	for(i=0;i<NumLatEdge;i++) {
		if (i==0) {
			for(j=0;j<NumLngEdge;j++) {
				k=j;
				sphere.fc[k].NumberofVertices=3;
				sphere.fc[k].pnt[0]=0;
				sphere.fc[k].pnt[1]=k+1;
				if (j!=NumLngEdge-1) {
					sphere.fc[k].pnt[2]=k+2;
				} else {
					sphere.fc[k].pnt[2]=1;
				}
			}
		} else if (i<NumLatEdge-1) {
			for(j=0;j<NumLngEdge;j++) {
				k=NumLngEdge*i+j;
				sphere.fc[k].NumberofVertices=4;
				sphere.fc[k].pnt[0]=k+1-NumLngEdge;
				sphere.fc[k].pnt[1]=k+1;
				if (j!=NumLngEdge-1) {
					sphere.fc[k].pnt[2]=k+2;
					sphere.fc[k].pnt[3]=k+2-NumLngEdge;
				} else {
					sphere.fc[k].pnt[2]=NumLngEdge*i+1;
					sphere.fc[k].pnt[3]=NumLngEdge*i+1-NumLngEdge;
				}
			}
		} else {
			for(j=0;j<NumLngEdge;j++) {
				k=NumLngEdge*i+j;
				sphere.fc[k].NumberofVertices=3;
				sphere.fc[k].pnt[0]=k+1-NumLngEdge;
				sphere.fc[k].pnt[1]=(NumLatEdge-1)*NumLngEdge+1;
				if (j!=NumLngEdge-1) {
					sphere.fc[k].pnt[2]=k+2-NumLngEdge;
				} else {
					sphere.fc[k].pnt[2]=NumLngEdge*i+1-NumLngEdge;
				}
			}
		}
	}
}

void drawSquare(float x,float y,color_t fillcolor)
{
	point2D_t p[4];
	float d=3;
	p[0].x=x-d,p[0].y=y-d;
	p[1].x=x+d,p[1].y=y-d;
	p[2].x=x+d,p[2].y=y+d;
	p[3].x=x-d,p[3].y=y+d;
	fillPolygon(p,4,fillcolor);
}

void userdraw(void)
{
	static int tick=0;
	float theta=0.5;
	matrix3D_t tilting=rotationXMTX(theta)*rotationYMTX(-theta);
	setColor(1,1,1);
	drawAxes(tilting);

	int tick1=500; // 0-tick1 setting vertices with increasing
	int tick2=600; // tick1-tick2 showing vertices
	int tick3=700; // tick2-tick3 showing mesh
	int tick4=1200; // tick3-tick4 putting flat face
	int tick5=1300; // tick4-tick5 showing flat face
	int tick6=1800; // tick5-tick6 putting smooth face
	int tick7=1900; // tick6-tick7 showing smooth face
	static object3D_t sphere0;
	static object3D_t sphere;
	vector3D_t unitz={0,0,1,1};
	vector3D_t LightVector0=rotationYMTX(0.2)*rotationXMTX(-1.2)*unitz;
	vector3D_t LightVector=tilting*LightVector0;
	vector3D_t ViewVector={0,0,1,1};
	vector3D_t ZeroVector={0,0,0,1};
	color_t ColCyan={0,1,1};
	static vector3D_t vec[1024];
	static vector3D_t nrmvec[1024];
	static vector3D_t vecbuff[8];
	static vector3D_t nrmvecbuff[1024];
	vector3D_t NormalVector;
	static point2D_t buff[1024];
	static point2D_t pbuff[1024];
	static color_t colbuff[1024];
	float normalzi;//z of normal vector of the face
	static float normalz[1024];//z of normal vector of the face
	color_t gray={0.7,0.7,0.7};
	color_t green={0,1,0};
	color_t fillcolor;

	int i,j;
	if (tick==0) {
		makeSpherePointData(150,sphere0);
	}
	sphere=sphere0;
	setColor(0,1,1);

	for (i=0;i<sphere.NumberofVertices;i++) {
		vec[i]=Point2Vector(sphere.pnt[i]);
		vec[i]=tilting*vec[i];
		buff[i]=Vector2Point2D(vec[i]);
	}
	for(i=0,j=0;i<sphere.NumberofVertices;j++){
		if ((i<2)||((i%NumLngEdge)!=1)) {
			pbuff[j]=buff[i];
			i++;
		} else{
			pbuff[j]=buff[i-NumLngEdge];
			j++;
			pbuff[j]=buff[i];
			i++;
		}
	}
	int nDispPoint=j;

	if (tick<tick2) {
		glutSetWindowTitle("sphere - making semicircular point data");
		int disp=(float)tick/tick1*nDispPoint;
		if (nDispPoint<disp) disp=nDispPoint;
		drawPolyline(pbuff,disp);
		for (i=0;i<disp;i++) {
			fillcolor=gray;
			if (i%(NumLngEdge+1)==0||i%(NumLngEdge+1)==1) fillcolor=green;
			drawSquare(pbuff[i].x,pbuff[i].y,fillcolor);
		}
	} else if (((tick2<tick)&&(tick<tick4))||((tick5<tick)&&(tick<tick6))){
		if (tick2<tick && tick<tick3) {
			glutSetWindowTitle("sphere - showing wire frame model with blue hidden faces");
		}
		setColor(0,0,1);
		for (i=0;i<sphere.NumberofFaces;i++) {
			for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
				vecbuff[j]=vec[sphere.fc[i].pnt[j]];
			}
			NormalVector=(vecbuff[1]-vecbuff[0])^(vecbuff[2]-vecbuff[0]);
			normalz[i]=NormalVector.v[2];
			if (normalz[i]<0.) {
				for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
					buff[j]=Vector2Point2D(vecbuff[j]);
				}
				drawPolygon(buff,sphere.fc[i].NumberofVertices);
			}
		}
		if (tick2<tick && tick<tick3) {
			setColor(0,1,1);
			for (i=0;i<sphere.NumberofFaces;i++) {
				if (0.<=normalz[i]) {
					for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
						vecbuff[j]=vec[sphere.fc[i].pnt[j]];
						buff[j]=Vector2Point2D(vecbuff[j]);
					}
					drawPolygon(buff,sphere.fc[i].NumberofVertices);
				}
			}
		}
	}

	if (tick3<tick && tick<tick5) {
	//removing hidden faces and Phong model Flat shading
		int disp=(float)(tick-tick3)/(tick4-tick3)*sphere.NumberofFaces;
		if (sphere.NumberofFaces<disp) disp=sphere.NumberofFaces;
		glutSetWindowTitle("sphere - removing hidden faces and Phong model Flat shading");
		for (i=0;i<disp;i++) {
			for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
				vecbuff[j]=vec[sphere.fc[i].pnt[j]];
			}
			NormalVector=(vecbuff[1]-vecbuff[0])^(vecbuff[2]-vecbuff[0]);
			normalzi=NormalVector.v[2];
			if (0.<normalzi) {
				for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
					buff[j]=Vector2Point2D(vecbuff[j]);
				}
				NormalVector=unitVector(NormalVector);
				fillcolor=PhongModel(LightVector,NormalVector,ViewVector,ColCyan);
				fillPolygon(buff,sphere.fc[i].NumberofVertices,fillcolor);
			} else {
				disp++;
				if (sphere.NumberofFaces<disp) disp=sphere.NumberofFaces;
			}
		}
		setColor(0,1,1);
		for (;i<sphere.NumberofFaces;i++) {
			if (0.<=normalz[i]) {
				for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
					vecbuff[j]=vec[sphere.fc[i].pnt[j]];
					buff[j]=Vector2Point2D(vecbuff[j]);
				}
				drawPolygon(buff,sphere.fc[i].NumberofVertices);
			}
		}
	}
	
	if (tick5<tick) {
		//removing hidden faces and Phong model Gouraud shading
		int disp=(float)(tick-tick5)/(tick6-tick5)*sphere.NumberofFaces;
		if (sphere.NumberofFaces<disp) disp=sphere.NumberofFaces;
		glutSetWindowTitle("sphere - removing hidden faces and Phong model Gouraud shading");
		ZeroVector=tilting*ZeroVector;
		for (i=0;i<sphere.NumberofFaces;i++) {
			vec[i]=Point2Vector(sphere.pnt[i]);
			vec[i]=tilting*vec[i];
			nrmvec[i]=tilting*sphere.NormalVector[i]-ZeroVector;
		}
		for (i=0;i<disp;i++) {
			for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
				vecbuff[j]=vec[sphere.fc[i].pnt[j]];
				nrmvecbuff[j]=nrmvec[sphere.fc[i].pnt[j]];
			}
			NormalVector=(vecbuff[1]-vecbuff[0])^(vecbuff[2]-vecbuff[0]);
			normalzi=NormalVector.v[2];
			if (0.<normalzi) {
				for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
					buff[j]=Vector2Point2D(vecbuff[j]);
					NormalVector=unitVector(nrmvecbuff[j]);
					colbuff[j]=PhongModel(LightVector,NormalVector,ViewVector,ColCyan);
				}
				gradatePolygon(buff,colbuff,sphere.fc[i].NumberofVertices);
			} else {
				disp++;
				if (sphere.NumberofFaces<disp) disp=sphere.NumberofFaces;
			}
		}
		setColor(0,1,1);
		for (;i<sphere.NumberofFaces;i++) {
			if (0.<=normalz[i]) {
				for (j=0;j<sphere.fc[i].NumberofVertices;j++) {
					vecbuff[j]=vec[sphere.fc[i].pnt[j]];
					buff[j]=Vector2Point2D(vecbuff[j]);
				}
				drawPolygon(buff,sphere.fc[i].NumberofVertices);
			}
		}
	}
	if (tick7<tick) tick=0;
	tick++;
}

int main(int argc, char **argv)
{
	glutInit(&argc,argv);
	glutInitDisplayMode ( GLUT_DOUBLE | GLUT_RGB );
	glutInitWindowPosition(100,100);
	glutInitWindowSize(640,480);
	glutCreateWindow ("sphere");
	glClearColor(0.0, 0.0, 0.0, 0.0);
	gluOrtho2D(-320., 320., -240.0, 240.0);
	  // Define the dimensions of the Orthographic Viewing Volume
	glutIdleFunc(display); // idle event call back
	glutDisplayFunc(display);
	glutMainLoop();
	return 0;
}