////////////// OpenGL CGCore for 2D Functions ver 1 /////////////////
#include <stdio.h>
#include <math.h>
#include <GL/glut.h>
#include "CGCore3D.h"

void C_Point::draw()
{
	color.set();
	glBegin(GL_POINTS);
		glVertex2f(x, y);
	glEnd();
}

void C_Vertex::resize() //resize to unit size
{
	float length=sqrt(x*x+y*y+z*z);
	float invlength;
	if (length!=0.) {
		invlength=1./length;
		x*=invlength;
		y*=invlength;
		z*=invlength;
	}
}

void C_Vertex::homogenize()
{
	float inv=1./t;
	x*=inv;
	y*=inv;
	z*=inv;
	t=1.;
}

C_Line::C_Line(float x1,float y1,float z1,float x2,float y2,float z2,C_Color col)
{
	p.x=x1;p.y=y1;p.z=z1;
	q.x=x2;q.y=y2;q.z=z2;
	color=col;
	p.t=1.;q.t=1.;
}

void C_Line::draw()
{
	color.set();
	glBegin(GL_LINES);
		glVertex2f(p.x, p.y);
		glVertex2f(q.x, q.y);
	glEnd();
}

void C_Line::homogenize()
{
	p.homogenize();
	q.homogenize();
}

C_Polyline::C_Polyline(int numberPoints0)
{
	numberPoints=numberPoints0;
	p=new C_Vertex[numberPoints];
}


C_Polyline::C_Polyline(C_Vertex *p0,int numberPoints0)
{
	int i;
	numberPoints=numberPoints0;
	p=new C_Vertex[numberPoints];
	for (i=0;i<numberPoints;i++) p[i]=p0[i];
}

C_Polyline::C_Polyline(C_Vertex *p0,int numberPoints0,C_Color color0)
{
	int i;
	numberPoints=numberPoints0;
	p=new C_Vertex[numberPoints];
	for (i=0;i<numberPoints;i++) p[i]=p0[i];
	color=color0;
}

C_Polyline::C_Polyline(const C_Polyline &pl)
{
	int i;
	numberPoints=pl.numberPoints;
	p=new C_Vertex[numberPoints];
	for (i=0;i<numberPoints;i++) p[i]=pl.p[i];
	color=pl.color;
}

C_Polyline::~C_Polyline()
{
	if (numberPoints) delete[] p;
	numberPoints=0;
}

C_Polyline C_Polyline::operator =(const C_Polyline &pl)
{
	int i;
	if (numberPoints) delete[] p;
	numberPoints=pl.numberPoints;
	p=new C_Vertex[numberPoints];
	for (i=0;i<numberPoints;i++) p[i]=pl.p[i];
	color=pl.color;
	return *this;
}

void C_Polyline::draw()
{
	int i;
	color.set();
	glBegin(GL_LINE_STRIP);
		for (i=0;i<numberPoints;i++) {
			glVertex2f(p[i].x, p[i].y);
		}
	glEnd();
}

void C_Polyline::homogenize()
{
	int i;
	for (i=0;i<numberPoints;i++) {
		p[i].homogenize();
	}
}

void C_Polyline::resize()
{
	int i;
	for (i=0;i<numberPoints;i++) {
		p[i].resize();
	}
}

C_Polygon::C_Polygon(int numberPoints0,int enableColorArray)
{
	numberPoints=numberPoints0;
	p=new C_Vertex[numberPoints];
	if (enableColorArray) colorarray=new C_Color[numberPoints];
	else colorarray=NULL;
}

C_Polygon::C_Polygon(C_Vertex *p0,int numberPoints0)
{
	int i;
	numberPoints=numberPoints0;
	p=new C_Vertex[numberPoints];
	for (i=0;i<numberPoints;i++) p[i]=p0[i];
	colorarray=NULL;
	C_Color white(1.,1.,1.);
	color=white;
}

C_Polygon::C_Polygon(C_Vertex *p0,int numberPoints0,C_Color color0)
{
	int i;
	numberPoints=numberPoints0;
	p=new C_Vertex[numberPoints];
	for (i=0;i<numberPoints;i++) p[i]=p0[i];
	colorarray=NULL;
	color=color0;
}

C_Polygon::C_Polygon(C_Vertex *p0,C_Color *colorarray0,int numberPoints0)
{
	int i;
	numberPoints=numberPoints0;
	p=new C_Vertex[numberPoints];
	colorarray=new C_Color[numberPoints];
	for (i=0;i<numberPoints;i++) { p[i]=p0[i]; colorarray[i]=colorarray0[i]; }
	C_Color white(1.,1.,1.);
	color=white;
}

C_Polygon::C_Polygon(const C_Polygon &pg)
{
	int i;
	color=pg.color;
	if (pg.colorarray!=NULL) {
		numberPoints=pg.numberPoints;
		p=new C_Vertex[numberPoints];
		colorarray=new C_Color[numberPoints];
		for (i=0;i<numberPoints;i++) { p[i]=pg.p[i]; colorarray[i]=pg.colorarray[i]; }
	} else {
		numberPoints=pg.numberPoints;
		p=new C_Vertex[numberPoints];
		for (i=0;i<numberPoints;i++) p[i]=pg.p[i];
		colorarray=NULL;
	}
}

C_Polygon::~C_Polygon()
{
	if (numberPoints) delete[] p;
	if (colorarray!=NULL) delete[] colorarray;
	numberPoints=0;
}

C_Polygon C_Polygon::operator =(const C_Polygon &pg)
{
	int i;
	if (numberPoints) delete[] p;
	if (colorarray!=NULL) delete[] colorarray;
	color=pg.color;
	if (pg.colorarray!=NULL) {
		numberPoints=pg.numberPoints;
		p=new C_Vertex[numberPoints];
		colorarray=new C_Color[numberPoints];
		for (i=0;i<numberPoints;i++) { p[i]=pg.p[i]; colorarray[i]=pg.colorarray[i]; }
	} else {
		numberPoints=pg.numberPoints;
		p=new C_Vertex[numberPoints];
		for (i=0;i<numberPoints;i++) p[i]=pg.p[i];
		colorarray=NULL;
	}
	return *this;
}

void C_Polygon::draw()
{
	int i;
	color.set();
	glBegin(GL_LINE_LOOP);
		for (i=0;i<numberPoints;i++) {
			glVertex2f(p[i].x, p[i].y);
		}
	glEnd();
}

void C_Polygon::fill()
{
	int i;
	color.set();
	glBegin(GL_POLYGON);
		for (i=0;i<numberPoints;i++) {
			glVertex2f(p[i].x, p[i].y);
		}
	glEnd();
}

void C_Polygon::gradate()
{
	int i;
	glBegin(GL_POLYGON);
		for (i=0;i<numberPoints;i++) {
			colorarray[i].set();
			glVertex2f(p[i].x, p[i].y);
		}
	glEnd();
}

void C_Polygon::homogenize()
{
	int i;
	for (i=0;i<numberPoints;i++) {
		p[i].homogenize();
	}
}

C_Matrix::C_Matrix()
{
	int i,j;
	for (i=0;i<4;i++) {
		for(j=0;j<4;j++) m[i][j]=0.;
		m[i][i]=1.;
	}
}

C_Matrix::C_Matrix(const C_Matrix &mx)
{
	int i,j;
	for (i=0;i<4;i++) for(j=0;j<4;j++) m[i][j]=mx.m[i][j];
}

C_Matrix::C_Matrix(const float m0[4][4])
{
	int i,j;
	for (i=0;i<4;i++) for(j=0;j<4;j++) m[i][j]=m0[i][j];
}

C_Matrix C_Matrix::operator =(const C_Matrix &mx)
{
	int i,j;
	for (i=0;i<4;i++) for(j=0;j<4;j++) m[i][j]=mx.m[i][j];
	return *this;
}

void C_Matrix::createIdentity()
{
	int i,j;
	for (i=0;i<4;i++) {
		for(j=0;j<4;j++) m[i][j]=0.;
		m[i][i]=1.;
	}
}

void C_Matrix::setTranslation(C_Vector v)
{
	createIdentity();
	m[0][3]=v.x;
	m[1][3]=v.y;
	m[2][3]=v.z;
}

void C_Matrix::setTranslation(float dx,float dy,float dz)
{
	createIdentity();
	m[0][3]=dx;
	m[1][3]=dy;
	m[2][3]=dz;
}

void C_Matrix::setRotationX(float theta)
{
	createIdentity();
	float cs=cos(theta);
	float sn=sin(theta);
	m[1][1]=cs; m[1][2]=-sn;
	m[2][1]=sn; m[2][2]=cs;
}

void C_Matrix::setRotationY(float theta)
{
	createIdentity();
	float cs=cos(theta);
	float sn=sin(theta);
	m[0][0]=cs; m[0][2]=sn;
	m[2][0]=-sn; m[2][2]=cs;
}

void C_Matrix::setRotationZ(float theta)
{
	createIdentity();
	float cs=cos(theta);
	float sn=sin(theta);
	m[0][0]=cs; m[0][1]=-sn;
	m[1][0]=sn; m[1][1]=cs;
}

void C_Matrix::setScale(float fx,float fy,float fz)
{
	createIdentity();
	m[0][0]=fx;
	m[1][1]=fy;
	m[2][2]=fz;
}

void C_Matrix::setPerspective(float zeye)
{
	createIdentity();
	m[3][2]=-1./zeye;
}

void C_Face::createFace(int VertexNumber0[],int numberV)
{
	int i;
	numberVertices=numberV;
	VertexNumber=new int[numberV];
	for (i=0;i<numberV;i++) VertexNumber[i]=VertexNumber0[i];
}

C_Face::C_Face(const C_Face &fc)
{
	int i;
	numberVertices=fc.numberVertices;
	VertexNumber=new int[numberVertices];
	for (i=0;i<numberVertices;i++) VertexNumber[i]=fc.VertexNumber[i];
}

C_Face::~C_Face()
{
	if (numberVertices) {
		delete[] VertexNumber;
		numberVertices=0;
	}
}

C_Face C_Face::operator =(const C_Face &fc)
{
	int i;
	if (numberVertices) delete[] VertexNumber;
	numberVertices=fc.numberVertices;
	VertexNumber=new int[numberVertices];
	for (i=0;i<numberVertices;i++) VertexNumber[i]=fc.VertexNumber[i];
	return *this;
}

void C_Polyhedron::createPolyhedron(C_Vertex Vertex0[],int numberV,C_Face Face0[],int numberF,C_Color color0,opticalAttribute_t &attr0)
{
	int i;
	numberVertices=numberV;
	Vertex=new C_Vertex[numberV];
	for (i=0;i<numberV;i++) Vertex[i]=Vertex0[i];
	numberFaces=numberF;
	Face=new C_Face[numberF];
	for (i=0;i<numberF;i++) {
		Face[i]=Face0[i];
	}
	color=color0;
	attr=attr0;
}

C_Polyhedron::C_Polyhedron(const C_Polyhedron &ph)
{
	int i;
	numberVertices=ph.numberVertices;
	Vertex=new C_Vertex[numberVertices];
	for (i=0;i<numberVertices;i++) Vertex[i]=ph.Vertex[i];
	numberFaces=ph.numberFaces;
	Face=new C_Face[numberFaces];
	for (i=0;i<numberFaces;i++) Face[i]=ph.Face[i];
	color=ph.color;
	attr=ph.attr;
}

C_Polyhedron::~C_Polyhedron()
{
	int i;
	if (numberVertices) {
		delete[] Vertex;
		numberVertices=0;
	}
	if (numberFaces) {
		for (i=0;i<numberFaces;i++) {
			if (Face[i].numberVertices) {
				Face[i].~C_Face();
			}
		}
		delete[] Face;
		numberFaces=0;
	}
}

C_Polyhedron C_Polyhedron::operator =(const C_Polyhedron &ph)
{
	int i;
	if (numberVertices) delete[] Vertex;
	if (numberFaces) {
		for (i=0;i<numberFaces;i++) {
			Face[i].~C_Face();
		}
		delete[] Face;
	}
	numberVertices=ph.numberVertices;
	Vertex=new C_Vertex[numberVertices];
	for (i=0;i<numberVertices;i++) Vertex[i]=ph.Vertex[i];
	numberFaces=ph.numberFaces;
	Face=new C_Face[numberFaces];
	for (i=0;i<numberFaces;i++) Face[i]=ph.Face[i];
	color=ph.color;
	attr=ph.attr;
	return *this;
}

void C_Polyhedron::homogenize()
{
	int i;
	for (i=0;i<numberVertices;i++) {
		Vertex[i].homogenize();
	}
}

void C_Smoothpolyhedron::createSmoothPolyhedron(C_Vertex Vertex0[],C_Vector NormalVector0[],int numberV,C_Face Face0[],int numberF,C_Color color0,opticalAttribute_t &attr0)
{
	int i;
	numberVertices=numberV;
	Vertex=new C_Vertex[numberV];
	NormalVector=new C_Vector[numberV];
	for (i=0;i<numberV;i++) {
		Vertex[i]=Vertex0[i];
		NormalVector[i]=NormalVector0[i];
	}
	numberFaces=numberF;
	Face=new C_Face[numberF];
	for (i=0;i<numberF;i++) {
		Face[i]=Face0[i];
	}
	color=color0;
	attr=attr0;
}

C_Smoothpolyhedron::C_Smoothpolyhedron(const C_Smoothpolyhedron &ph)
{
	int i;
	numberVertices=ph.numberVertices;
	Vertex=new C_Vertex[numberVertices];
	NormalVector=new C_Vector[numberVertices];
	for (i=0;i<numberVertices;i++) {
		Vertex[i]=ph.Vertex[i];
		NormalVector[i]=ph.NormalVector[i];
	}
	numberFaces=ph.numberFaces;
	Face=new C_Face[numberFaces];
	for (i=0;i<numberFaces;i++) Face[i]=ph.Face[i];
	color=ph.color;
	attr=ph.attr;
}

C_Smoothpolyhedron::~C_Smoothpolyhedron()
{
	int i;
	if (numberVertices) {
		delete[] Vertex;
		delete[] NormalVector;
		numberVertices=0;
	}
	if (numberFaces) {
		for (i=0;i<numberFaces;i++) {
			if (Face[i].numberVertices) {
				Face[i].~C_Face();
			}
		}
		delete[] Face;
		numberFaces=0;
	}
}

C_Smoothpolyhedron C_Smoothpolyhedron::operator =(const C_Smoothpolyhedron &ph)
{
	int i;
	if (numberVertices) {
		delete[] Vertex;
		delete[] NormalVector;
	}
	if (numberFaces) {
		for (i=0;i<numberFaces;i++) {
			Face[i].~C_Face();
		}
		delete[] Face;
	}
	numberVertices=ph.numberVertices;
	Vertex=new C_Vertex[numberVertices];
	NormalVector=new C_Vector[numberVertices];
	for (i=0;i<numberVertices;i++) {
		Vertex[i]=ph.Vertex[i];
		NormalVector[i]=ph.NormalVector[i];
	}
	numberFaces=ph.numberFaces;
	Face=new C_Face[numberFaces];
	for (i=0;i<numberFaces;i++) Face[i]=ph.Face[i];
	color=ph.color;
	attr=ph.attr;
	return *this;
}

void C_Smoothpolyhedron::resize()
{
	int i;
	for (i=0;i<numberVertices;i++) {
		NormalVector[i].resize();
	}
}

C_Scene::C_Scene()
{
	C_Vector unitz(0,0,1);
	C_Matrix roty,rotx;
	ViewDirection=unitz;
	zeye=1000;
	theta=-0.5;
	phi=0.5;
	rotx.setRotationX(phi);
	roty.setRotationY(theta);
	Tilting=rotx*roty;
	rotx.setRotationX(-1.2);
	roty.setRotationY(0.2);
	LightDirection=Tilting*roty*rotx*unitz;
	Perspective.setPerspective(zeye);
	stackpointer=0;
}

C_Scene::C_Scene(float theta0,float phi0,float zeye0,C_Vector light)
{
	C_Vector unitz(0,0,1);
	C_Matrix roty,rotx;
	theta=theta0;
	phi=phi0;
	zeye=zeye0;
	ViewDirection=unitz;
	rotx.setRotationX(phi);
	roty.setRotationY(theta);
	Tilting=rotx*roty;
	LightDirection=Tilting*light;
	Perspective.setPerspective(zeye);
	stackpointer=0;
}

void C_Scene::initializeTMatrix()
{
	C_Matrix mat;
	Transformation=TransformationN=mat;
}

void C_Scene::setTranslation(C_Vector v)
{
	C_Matrix mat;
	mat.setTranslation(v);
	Transformation=Transformation*mat;
}

void C_Scene::setTranslation(float dx,float dy,float dz)
{
	C_Matrix mat;
	mat.setTranslation(dx,dy,dz);
	Transformation=Transformation*mat;
}

void C_Scene::setRotationX(float theta)
{
	C_Matrix mat;
	mat.setRotationX(theta);
	Transformation=Transformation*mat;
	TransformationN=TransformationN*mat;
}

void C_Scene::setRotationY(float theta)
{
	C_Matrix mat;
	mat.setRotationY(theta);
	Transformation=Transformation*mat;
	TransformationN=TransformationN*mat;
}

void C_Scene::setRotationZ(float theta)
{
	C_Matrix mat;
	mat.setRotationZ(theta);
	Transformation=Transformation*mat;
	TransformationN=TransformationN*mat;
}

void C_Scene::setScale(float fx,float fy,float fz)
{
	C_Matrix mat;
	mat.setScale(fx,fy,fz);
	Transformation=Transformation*mat;
	mat.setScale(1/fx,1/fy,1/fz);
	TransformationN=TransformationN*mat;
}

void C_Scene::pushMatrix()
{
	if (stackpointer==32) return;
	Stack[stackpointer]=Transformation;
	StackN[stackpointer]=TransformationN;
	stackpointer++;
}

void C_Scene::popMatrix()
{
	stackpointer--;
	if (stackpointer<0) stackpointer=0;
	Transformation=Stack[stackpointer];
	TransformationN=StackN[stackpointer];
}

static void drawcharX(float x,float y,float z,C_Color col)
{
	C_Line Line1(x,y,z,x+10,y+12,z,col);Line1.draw();
	C_Line Line2(x,y+12,z,x+10,y,z,col);Line2.draw();
}

static void drawcharY(float x,float y,float z,C_Color col)
{
	C_Line Line1(x+5,y,z,x+5,y+7,z,col);Line1.draw();
	C_Line Line2(x,y+12,z,x+5,y+7,z,col);Line2.draw();
	C_Line Line3(x+10,y+12,z,x+5,y+7,z,col);Line3.draw();
}

static void drawcharZ(float x,float y,float z,C_Color col)
{
	C_Line Line1(x,y+12,z,x+10,y+12,z,col);Line1.draw();
	C_Line Line2(x+10,y+12,z,x,y,z,col);Line2.draw();
	C_Line Line3(x,y,z,x+10,y,z,col);Line3.draw();
}

void C_Scene::drawAxes()
{
#define HALFAXIS  220
#define HALFAXIS1 (HALFAXIS-10)
	C_Color White(1.,1.,1.);
	float Axes00[14][3]={
		{-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}
	};
	C_Vertex Axes0[14];
	int i;
	for (i=0;i<14;i++) {
		Axes0[i].x=Axes00[i][0];
		Axes0[i].y=Axes00[i][1];
		Axes0[i].z=Axes00[i][2];
	}
	C_Polyline Axes(Axes0,14,White);
	C_Polyline Axes1=Perspective*Tilting*Axes;
	Axes1.homogenize();

	Axes1.draw();
	drawcharX(Axes1.p[1].x,Axes1.p[1].y,Axes1.p[1].z,White);
	drawcharY(Axes1.p[6].x,Axes1.p[6].y,Axes1.p[6].z,White);
	drawcharZ(Axes1.p[11].x-14,Axes1.p[11].y,Axes1.p[11].z,White);
}

void C_Scene::drawPolyhedronWire(C_Polyhedron &ph)
{

	C_Matrix mat=Perspective*Tilting*Transformation;
	C_Color visiblecolor=ph.color;
	C_Color hiddencolor=0.5*ph.color;
	C_Vector NormalVector;
	int i,j,num;

	C_Polyline vertex(ph.numberVertices);
	for (i=0;i<ph.numberVertices;i++) vertex.p[i]=mat*ph.Vertex[i];
	vertex.homogenize();
	num=0;
	for (i=0;i<ph.numberFaces;i++) if (num<ph.Face[i].numberVertices) num=ph.Face[i].numberVertices;
	C_Polygon pg(num,0);
	for (i=0;i<ph.numberFaces;i++) {
		num=ph.Face[i].numberVertices;
		for (j=0;j<num;j++) {
			pg.p[j]=vertex.p[ph.Face[i].VertexNumber[j]];
		}
		NormalVector=(pg.p[1]-pg.p[0])^(pg.p[2]-pg.p[0]);
		if (NormalVector.z<0.) {
			pg.numberPoints=num;
			pg.setColor(hiddencolor);
			pg.draw();
		}
	}
	for (i=0;i<ph.numberFaces;i++) {
		num=ph.Face[i].numberVertices;
		for (j=0;j<num;j++) {
			pg.p[j]=vertex.p[ph.Face[i].VertexNumber[j]];
		}
		NormalVector=(pg.p[1]-pg.p[0])^(pg.p[2]-pg.p[0]);
		if (0.<NormalVector.z) {
			pg.numberPoints=num;
			pg.setColor(visiblecolor);
			pg.draw();
		}
	}
}

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

// x to yth power
static 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));
}

//PhongModel color calculation
// LightVector, NormalVector, ViewVector, ColorofObject
C_Color C_Scene::PhongModel(C_Vector Normal,C_Color col,opticalAttribute_t &attr)
{
	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;
	C_Color ColWhite(1,1,1);
	C_Vector ReflectionVector=(2.*(LightDirection*Normal)*Normal)-LightDirection;
	tmp=Normal*LightDirection;
	NL=funcPositive(tmp);
	tmp=ReflectionVector*ViewDirection;
	RV=funcPositive(tmp);
	return attr.kdif*NL*col+attr.kspe*power(RV,4)*ColWhite+attr.kamb*col;
}

void C_Scene::drawPolyhedronSurface(C_Polyhedron &ph)
{
	C_Matrix mat1=Tilting*Transformation;  //for color calcularion
	C_Matrix mat2=Perspective*mat1;        //for draw
	C_Vector NormalVector1;
	C_Vector NormalVector2;
	int i,j,num;

	C_Polyline vertex1(ph.numberVertices);
	C_Polyline vertex2(ph.numberVertices);
	for (i=0;i<ph.numberVertices;i++) {
		vertex1.p[i]=mat1*ph.Vertex[i];
		vertex2.p[i]=mat2*ph.Vertex[i];
	}

	vertex2.homogenize();
	num=0;
	for (i=0;i<ph.numberFaces;i++) if (num<ph.Face[i].numberVertices) num=ph.Face[i].numberVertices;
	C_Polygon pg1(num,0);
	C_Polygon pg2(num,0);
	for (i=0;i<ph.numberFaces;i++) {
		int num=ph.Face[i].numberVertices;
		for (j=0;j<num;j++) {
			pg1.p[j]=vertex1.p[ph.Face[i].VertexNumber[j]];
			pg2.p[j]=vertex2.p[ph.Face[i].VertexNumber[j]];
		}
		NormalVector1=(pg1.p[1]-pg1.p[0])^(pg1.p[2]-pg1.p[0]);
		NormalVector2=(pg2.p[1]-pg2.p[0])^(pg2.p[2]-pg2.p[0]);
		if (0.<NormalVector2.z) {
			NormalVector1.resize();
			C_Color facecolor=PhongModel(NormalVector1,ph.color,ph.attr);
			pg2.numberPoints=num;
			pg2.setColor(facecolor);
			pg2.fill();
		}
	}
}

void C_Scene::drawPolyhedronSurfaceZSort(C_Polyhedron &ph)
{
	C_Matrix mat1=Tilting*Transformation;  //for color calcularion
	C_Matrix mat2=Perspective*mat1;        //for draw
	C_Vector NormalVector1;
	C_Vector NormalVector2;
	int i,j,k,num,numFaces;
	C_Polygon *polygon=new C_Polygon[ph.numberFaces];
	int *order=new int[ph.numberFaces];
	numFaces=0;

	C_Polyline vertex1(ph.numberVertices);
	C_Polyline vertex2(ph.numberVertices);
	for (i=0;i<ph.numberVertices;i++) {
		vertex1.p[i]=mat1*ph.Vertex[i];
		vertex2.p[i]=mat2*ph.Vertex[i];
	}

	vertex2.homogenize();
	num=0;
	for (i=0;i<ph.numberFaces;i++) if (num<ph.Face[i].numberVertices) num=ph.Face[i].numberVertices;
	C_Polygon pg1(num,0);
	C_Polygon pg2(num,0);
	for (i=0;i<ph.numberFaces;i++) {
		int num=ph.Face[i].numberVertices;
		for (j=0;j<num;j++) {
			pg1.p[j]=vertex1.p[ph.Face[i].VertexNumber[j]];
			pg2.p[j]=vertex2.p[ph.Face[i].VertexNumber[j]];
		}
		NormalVector1=(pg1.p[1]-pg1.p[0])^(pg1.p[2]-pg1.p[0]);
		NormalVector2=(pg2.p[1]-pg2.p[0])^(pg2.p[2]-pg2.p[0]);
		if (0.<NormalVector2.z) {
			NormalVector1.resize();
			C_Color facecolor=PhongModel(NormalVector1,ph.color,ph.attr);
			pg2.numberPoints=num;
			pg2.setColor(facecolor);
			polygon[numFaces++]=pg2;
		}
	}
	// z-sort
	for (i=0;i<numFaces;i++) order[i]=i;
	for (i=0;i<numFaces;i++) {
		for (j=numFaces-2;i<=j;j--) {
			k=j+1;
			if (polygon[order[k]].p[0].z<polygon[order[j]].p[0].z) {
				int tmp=order[k];
				order[k]=order[j];
				order[j]=tmp;
			}
		}
	}
	for (i=0;i<numFaces;i++) polygon[order[i]].fill();
	delete[] polygon;
	delete[] order;
}

void C_Scene::drawPolyhedronSmoothSurface(C_Smoothpolyhedron &ph)
{
	C_Matrix mat1=Tilting*TransformationN;  //for color calcularion
	C_Matrix mat2=Perspective*Tilting*Transformation;        //for draw
	C_Vector NormalVector;
	int i,j,num;

	C_Polyline vertex1(ph.numberVertices);
	C_Polyline vertex2(ph.numberVertices);
	for (i=0;i<ph.numberVertices;i++) {
		vertex1.p[i]=mat1*ph.NormalVector[i];
		vertex2.p[i]=mat2*ph.Vertex[i];
	}
	vertex1.resize();
	vertex2.homogenize();
	num=0;
	for (i=0;i<ph.numberFaces;i++) if (num<ph.Face[i].numberVertices) num=ph.Face[i].numberVertices;
	C_Polygon pg(num,1);
	for (i=0;i<ph.numberFaces;i++) {
		int num=ph.Face[i].numberVertices;
		for (j=0;j<num;j++) {
			pg.p[j]=vertex2.p[ph.Face[i].VertexNumber[j]];
		}
		NormalVector=(pg.p[1]-pg.p[0])^(pg.p[2]-pg.p[0]);
		if (0.<NormalVector.z) {
			for (j=0;j<num;j++) {
				NormalVector=vertex1.p[ph.Face[i].VertexNumber[j]];
				pg.colorarray[j]=PhongModel(NormalVector,ph.color,ph.attr);
			}
			pg.numberPoints=num;
			pg.gradate();
		}
	}
}

void C_Scene::drawPolyhedronSmoothSurfaceZSort(C_Smoothpolyhedron &ph)
{
	C_Matrix mat1=Tilting*TransformationN;  //for color calcularion
	C_Matrix mat2=Perspective*Tilting*Transformation;        //for draw
	C_Vector NormalVector;
	int i,j,k,num,numFaces;
	C_Polygon *polygon=new C_Polygon[ph.numberFaces];
	int *order=new int[ph.numberFaces];
	numFaces=0;

	C_Polyline vertex1(ph.numberVertices);
	C_Polyline vertex2(ph.numberVertices);
	for (i=0;i<ph.numberVertices;i++) {
		vertex1.p[i]=mat1*ph.NormalVector[i];
		vertex2.p[i]=mat2*ph.Vertex[i];
	}
	vertex1.resize();
	vertex2.homogenize();
	num=0;
	for (i=0;i<ph.numberFaces;i++) if (num<ph.Face[i].numberVertices) num=ph.Face[i].numberVertices;
	C_Polygon pg(num,1);
	for (i=0;i<ph.numberFaces;i++) {
		int num=ph.Face[i].numberVertices;
		for (j=0;j<num;j++) {
			pg.p[j]=vertex2.p[ph.Face[i].VertexNumber[j]];
		}
		NormalVector=(pg.p[1]-pg.p[0])^(pg.p[2]-pg.p[0]);
		if (0.<NormalVector.z) {
			for (j=0;j<num;j++) {
				NormalVector=vertex1.p[ph.Face[i].VertexNumber[j]];
				pg.colorarray[j]=PhongModel(NormalVector,ph.color,ph.attr);
			}
			pg.numberPoints=num;
			polygon[numFaces++]=pg;
		}
	}
	// z-sort
	for (i=0;i<numFaces;i++) order[i]=i;
	for (i=0;i<numFaces;i++) {
		for (j=numFaces-2;i<=j;j--) {
			k=j+1;
			if (polygon[order[k]].p[0].z<polygon[order[j]].p[0].z) {
				int tmp=order[k];
				order[k]=order[j];
				order[j]=tmp;
			}
		}
	}
	for (i=0;i<numFaces;i++) polygon[order[i]].gradate();
	delete[] polygon;
	delete[] order;
}

////////////////////////////////////////////////////////////////////////////////

C_Matrix operator * (const C_Matrix &a, const C_Matrix &b)
{
	C_Matrix 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;
}

C_Vertex operator * (const C_Matrix &a, const C_Vertex &b)
{
	C_Vertex c;//c=a*b
	c.x=a.m[0][0]*b.x+a.m[0][1]*b.y+a.m[0][2]*b.z+a.m[0][3]*b.t;
	c.y=a.m[1][0]*b.x+a.m[1][1]*b.y+a.m[1][2]*b.z+a.m[1][3]*b.t;
	c.z=a.m[2][0]*b.x+a.m[2][1]*b.y+a.m[2][2]*b.z+a.m[2][3]*b.t;
	c.t=a.m[3][0]*b.x+a.m[3][1]*b.y+a.m[3][2]*b.z+a.m[3][3]*b.t;
	return c;
}

C_Line operator * (const C_Matrix &a, const C_Line &b)
{
	C_Line c;//c=a*b
	c.p=a*b.p;
	c.q=a*b.q;
	c.color=b.color;
	return c;
}

C_Polyline operator * (const C_Matrix &a, const C_Polyline &b)
{
	C_Polyline c;//c=a*b
	int i;
	c=b;
	for (i=0;i<b.numberPoints;i++) c.p[i]=a*b.p[i];
	return c;
}

C_Polygon operator * (const C_Matrix &a, const C_Polygon &b)
{
	C_Polygon c;//c=a*b
	int i;
	c=b;
	for (i=0;i<b.numberPoints;i++) c.p[i]=a*b.p[i];
	return c;
}

C_Polyhedron operator * (const C_Matrix &a, const C_Polyhedron &b)
{
	C_Polyhedron c;//c=a*b
	int i;
	c=b;
	for (i=0;i<b.numberVertices;i++) c.Vertex[i]=a*b.Vertex[i];
	return c;
}

C_Color operator + (const C_Color &a, const C_Color &b)
{
	C_Color col;
	col.r=a.r+b.r;
	col.g=a.g+b.g;
	col.b=a.b+b.b;
	return col;
}

C_Color operator * (const float &a, const C_Color &b)
{
	C_Color col;
	col.r=a*b.r;
	col.g=a*b.g;
	col.b=a*b.b;
	return col;
}

C_Color operator * (const C_Color &a, const float &b)
{
	C_Color col;
	col.r=b*a.r;
	col.g=b*a.g;
	col.b=b*a.b;
	return col;
}

float operator * (const C_Vector &a, const C_Vector &b) //“à�Ï
{
	float c;//c=a*b
	c=a.x*b.x+a.y*b.y+a.z*b.z;
	return c;
}

C_Vector operator * (const float &a, const C_Vector &b)
{
	C_Vector c;//c=a*b
	c.x=a*b.x;
	c.y=a*b.y;
	c.z=a*b.z;
	c.t=1.;
	return c;
}

C_Vector operator * (const C_Vector &a, const float &b)
{
	C_Vector c;//c=a*b
	c.x=a.x*b;
	c.y=a.y*b;
	c.z=a.z*b;
	c.t=1.;
	return c;
}

// outer product (cross product ) of homogeneous vector
//       i         j         k
//       a0       a1        a2
//       b0       b1        b2
C_Vector operator ^ (const C_Vector &a, const C_Vector &b) //ŠO�Ï
{
	C_Vector c;//c=a*b
	c.x=a.y*b.z-a.z*b.y;
	c.y=a.z*b.x-a.x*b.z;
	c.z=a.x*b.y-a.y*b.x;
	c.t=1.;
	return c;
}

C_Vector operator + (const C_Vector &a, const C_Vector &b)
{
	C_Vector c;//c=a+b
	c.x=a.x+b.x;
	c.y=a.y+b.y;
	c.z=a.z+b.z;
	c.t=1.;
	return c;
}

C_Vector operator - (const C_Vector &a, const C_Vector &b)
{
	C_Vector c;//c=a-b
	c.x=a.x-b.x;
	c.y=a.y-b.y;
	c.z=a.z-b.z;
	c.t=1.;
	return c;
}

C_Vector operator - (const C_Vector &a)
{
	C_Vector c;//c=-a
	c.x=-a.x;
	c.y=-a.y;
	c.z=-a.z;
	c.t=1.;
	return c;
}