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


struct Raster {
	int size;
	int **rast;
};

struct Light {
	double col[3];
	double p[3];
};

struct Triangle {
	double col[3];
	double reflect;
	double alpha;
	double p[3][3];
	double norm[3];
	double min[3],max[3];
};

struct Sphere {
	double col[3];
	double reflect;
	double alpha;
	double c[3];
	double radius;
};

struct Geometry {
	int numLights;
	int numTris;
	int numSpheres;
	struct Light light[100];
	struct Triangle tri[10000];
	struct Sphere sphere[1000];
	int useReflect;
	int useAlpha;
	double index;
	double size;
};

double distance(double *p1,double *p2) {
	double dx=p1[0]-p2[0];
	double dy=p1[1]-p2[1];
	double dz=p1[2]-p2[2];
	return sqrt(dx*dx+dy*dy+dz*dz);
}

void normalize(double *p) {
	double len=sqrt(p[0]*p[0]+p[1]*p[1]+p[2]*p[2]);
	if(len==0.0) return;
	p[0]/=len;
	p[1]/=len;
	p[2]/=len;
}

double dotProduct(double *p1,double *p2) {
	return p1[0]*p2[0]+p1[1]*p2[1]+p1[2]*p2[2];
}

void crossProduct(double *p1,double *p2,double *prod) {
	prod[0]=p1[1]*p2[2]-p2[1]*p1[2];
	prod[1]=p1[2]*p2[0]-p2[2]*p1[0];
	prod[2]=p1[0]*p2[1]-p2[0]*p1[1];
}

void calcTriNorm(struct Triangle *tri) {
	double v1[3];
	double v2[3];
	int i;
	for(i=0; i<3; i++) {
		v1[i]=tri->p[0][i]-tri->p[1][i];
		v2[i]=tri->p[2][i]-tri->p[1][i];
	}
	crossProduct(v1,v2,tri->norm);
	normalize(tri->norm);
}

void setBounds(struct Triangle *tri) {
	int i,j;
	for(i=0; i<3; i++) {
		tri->min[i]=tri->p[0][i];
		tri->max[i]=tri->p[0][i];
		for(j=1; j<3; j++) {
			if(tri->p[j][i]<tri->min[i]) tri->min[i]=tri->p[j][i];
			if(tri->p[j][i]>tri->max[i]) tri->max[i]=tri->p[j][i];
		}
		tri->min[i]-=1e-4;
		tri->max[i]+=1e-4;
	}
}

void readPoint(FILE *fin,double p[]) {
	int i;
	for(i=0; i<3; i++) {
		fscanf(fin,"%lf",p+i);
	}
}

void readTriangle(FILE *fin,double tp[][3]) {
	readPoint(fin,tp[0]);
	readPoint(fin,tp[1]);
	readPoint(fin,tp[2]);
}

void writeGeometry(FILE *file,struct Geometry *geom) {
	int i;
	int ref,al;

	fprintf(file,"<?xml version=\"1.0\"?>\n");
	fprintf(file,"<geometry>\n");
	fprintf(file,"\t<eye x=\"0.5\" y=\"0.5\" z=\"-1.0\"/>\n");
	fprintf(file,"\t<rasterCorner x=\"0.0\" y=\"1.0\" z=\"0.0\"/>\n");
	fprintf(file,"\t<rasterRight x=\"1.0\" y=\"0.0\" z=\"0.0\"/>\n");
	fprintf(file,"\t<rasterDown x=\"0.0\" y=\"-1.0\" z=\"0.0\"/>\n");
	for(i=0; i<geom->numLights; i++) {
		fprintf(file,"\t<light>\n");
		fprintf(file,"\t\t<point x=\"%f\" y=\"%f\" z=\"%f\"/>\n",geom->light[i].p[0],geom->light[i].p[1],geom->light[i].p[2]);
		fprintf(file,"\t\t<color x=\"%f\" y=\"%f\" z=\"%f\"/>\n",geom->light[i].col[0],geom->light[i].col[1],geom->light[i].col[2]);
		fprintf(file,"\t</light>\n");
	}
	for(i=0; i<geom->numTris; i++) {
		fprintf(file,"\t<triangle>\n");
		fprintf(file,"\t\t<point x=\"%f\" y=\"%f\" z=\"%f\"/>\n",geom->tri[i].p[0][0],geom->tri[i].p[0][1],geom->tri[i].p[0][2]);
		fprintf(file,"\t\t<point x=\"%f\" y=\"%f\" z=\"%f\"/>\n",geom->tri[i].p[1][0],geom->tri[i].p[1][1],geom->tri[i].p[1][2]);
		fprintf(file,"\t\t<point x=\"%f\" y=\"%f\" z=\"%f\"/>\n",geom->tri[i].p[2][0],geom->tri[i].p[2][1],geom->tri[i].p[2][2]);
		fprintf(file,"\t\t<color red=\"%f\" green=\"%f\" blue=\"%f\"/>\n",geom->tri[i].col[0],geom->tri[i].col[1],geom->tri[i].col[2]);
		fprintf(file,"\t\t<reflect value=\"%f\"/>\n",geom->tri[i].reflect);
		fprintf(file,"\t</triangle>\n");
	}
	for(i=0; i<geom->numSpheres; i++) {
		fprintf(file,"\t<sphere>\n");
		fprintf(file,"\t\t<point x=\"%f\" y=\"%f\" z=\"%f\"/>\n",geom->sphere[i].c[0],geom->sphere[i].c[1],geom->sphere[i].c[2]);
		fprintf(file,"\t\t<color red=\"%f\" green=\"%f\" blue=\"%f\"/>\n",geom->sphere[i].col[0],geom->sphere[i].col[1],geom->sphere[i].col[2]);
		fprintf(file,"\t\t<reflect value=\"%f\"/>\n",geom->sphere[i].reflect);
		fprintf(file,"\t</sphere>\n");
	}
	fprintf(file,"</geometry>\n");
}

void linearInterp(double p1[],double p2[],double frac,double pout[]) {
	int i;
	for(i=0; i<3; i++) {
		pout[i]=p1[i]+frac*(p2[i]-p1[i]);
	}
}

void interpolateSheet(double i,double j,double sheet[][2][3],double num,double pout[]) {
	double p[4][3];
	linearInterp(sheet[0][0],sheet[1][0],i/num,p[0]);
	linearInterp(sheet[0][1],sheet[1][1],i/num,p[1]);
	linearInterp(sheet[0][0],sheet[0][1],j/num,p[2]);
	linearInterp(sheet[1][0],sheet[1][1],j/num,p[3]);

	linearInterp(p[0],p[1],j/num,p[0]);
	linearInterp(p[2],p[3],i/num,p[2]);

	linearInterp(p[0],p[2],0.5,pout);
}

void addSheetSquare(double i,double j,double sheet[][2][3],double num,double rgb[],double alpha,double reflect,struct Geometry *geom) {
	int tnum=geom->numTris;
	int k;
	double p[4][3];
	interpolateSheet(i,j,sheet,num,p[0]);
	interpolateSheet(i+1,j,sheet,num,p[1]);
	interpolateSheet(i+1,j+1,sheet,num,p[2]);
	interpolateSheet(i,j+1,sheet,num,p[3]);
	for(k=0; k<3; k++) {
		geom->tri[tnum].p[0][k]=p[0][k];
		geom->tri[tnum].p[1][k]=p[1][k];
		geom->tri[tnum].p[2][k]=p[3][k];
		geom->tri[tnum+1].p[0][k]=p[1][k];
		geom->tri[tnum+1].p[1][k]=p[2][k];
		geom->tri[tnum+1].p[2][k]=p[3][k];
	}
	geom->tri[tnum].col[0]=rgb[0];
	geom->tri[tnum].col[1]=rgb[1];
	geom->tri[tnum].col[2]=rgb[2];
	geom->tri[tnum].alpha=alpha;
	geom->tri[tnum].reflect=reflect;
	geom->tri[tnum+1].col[0]=rgb[0];
	geom->tri[tnum+1].col[1]=rgb[1];
	geom->tri[tnum+1].col[2]=rgb[2];
	geom->tri[tnum+1].alpha=alpha;
	geom->tri[tnum+1].reflect=reflect;
	geom->numTris+=2;
}

void combine(double v1[],double f1,double v2[],double f2,double ans[]) {
	int i;
	for(i=0; i<3; i++) {
		ans[i]=v1[i]*f1+v2[i]*f2;
	}
}

void calcCylinderPoints(double cylinder[][21][3],double numSides) {
	double r[3]={cylinder[0][1][0]-cylinder[0][0][0],cylinder[0][1][1]-cylinder[0][0][1],cylinder[0][1][2]-cylinder[0][0][2]};
	double n[3]={cylinder[1][0][0]-cylinder[0][0][0],cylinder[1][0][1]-cylinder[0][0][1],cylinder[1][0][2]-cylinder[0][0][2]};
	double rp[3];
	double add[3];
	double ang=0,delta=6.28/numSides;
	int i=1;

	normalize(n);
	crossProduct(r,n,rp);
	while(ang<6.27) {
		combine(r,cos(ang),rp,sin(ang),add);
		combine(cylinder[0][0],1.0,add,1.0,cylinder[0][i]);
		combine(cylinder[1][0],1.0,add,1.0,cylinder[1][i]);
		ang+=delta;
		i++;
	}
}

void addCylinderPiece(int i,double cylinder[][21][3],int numSides,double rgb[],double alpha,double reflect,struct Geometry *geom) {
	int tnum=geom->numTris;
	int k;
	int ip=i+1;
	if(ip>numSides) ip=1;
	for(k=0; k<3; k++) {
		geom->tri[tnum].p[0][k]=cylinder[0][0][k];
		geom->tri[tnum].p[1][k]=cylinder[0][i][k];
		geom->tri[tnum].p[2][k]=cylinder[0][ip][k];
		geom->tri[tnum+1].p[0][k]=cylinder[1][i][k];
		geom->tri[tnum+1].p[1][k]=cylinder[0][ip][k];
		geom->tri[tnum+1].p[2][k]=cylinder[0][i][k];
		geom->tri[tnum+2].p[0][k]=cylinder[1][i][k];
		geom->tri[tnum+2].p[1][k]=cylinder[1][ip][k];
		geom->tri[tnum+2].p[2][k]=cylinder[0][ip][k];
		geom->tri[tnum+3].p[0][k]=cylinder[1][0][k];
		geom->tri[tnum+3].p[1][k]=cylinder[1][ip][k];
		geom->tri[tnum+3].p[2][k]=cylinder[1][i][k];
	}
	geom->tri[tnum].col[0]=rgb[0];
	geom->tri[tnum].col[1]=rgb[1];
	geom->tri[tnum].col[2]=rgb[2];
	geom->tri[tnum].alpha=alpha;
	geom->tri[tnum].reflect=reflect;
	geom->tri[tnum+1].col[0]=rgb[0];
	geom->tri[tnum+1].col[1]=rgb[1];
	geom->tri[tnum+1].col[2]=rgb[2];
	geom->tri[tnum+1].alpha=alpha;
	geom->tri[tnum+1].reflect=reflect;
	geom->tri[tnum+2].col[0]=rgb[0];
	geom->tri[tnum+2].col[1]=rgb[1];
	geom->tri[tnum+2].col[2]=rgb[2];
	geom->tri[tnum+2].alpha=alpha;
	geom->tri[tnum+2].reflect=reflect;
	geom->tri[tnum+3].col[0]=rgb[0];
	geom->tri[tnum+3].col[1]=rgb[1];
	geom->tri[tnum+3].col[2]=rgb[2];
	geom->tri[tnum+3].alpha=alpha;
	geom->tri[tnum+3].reflect=reflect;
	geom->numTris+=4;
}

double posRand(double p[],int which) {
	int i=0xffffff*p[0];
	int j=0xffffff*p[1];
	int k=0xffffff*p[2];
	int ret=((~i)<<10)+(j>>4)^((k<<12)-(i>>8)+((~j)<<16))+(i*which*k);
	return (double)ret/0x7fffffff;
}

void copyVect(double dest[],double src[]) {
	dest[0]=src[0];
	dest[1]=src[1];
	dest[2]=src[2];
}

void addFract(double fractTri[][3],double fractOffset,int level,double rgb[],double alpha,double reflect,struct Geometry *geom) {
	if(level==0) {
		int tnum=geom->numTris;
		int k;
		for(k=0; k<3; k++) {
			geom->tri[tnum].p[0][k]=fractTri[0][k];
			geom->tri[tnum].p[1][k]=fractTri[1][k];
			geom->tri[tnum].p[2][k]=fractTri[2][k];
		}
		geom->tri[tnum].col[0]=rgb[0];
		geom->tri[tnum].col[1]=rgb[1];
		geom->tri[tnum].col[2]=rgb[2];
		geom->tri[tnum].alpha=alpha;
		geom->tri[tnum].reflect=reflect;
		geom->numTris+=1;
	} else {
		double tri[3][3];
		double mids[3][3];
		int i,j;

		for(i=0; i<3; i++) {
			linearInterp(fractTri[i],fractTri[(i+1)%3],0.5,mids[i]);
			for(j=0; j<3; j++) {
				mids[i][j]+=posRand(mids[i],j)*fractOffset;
			}
		}
		copyVect(tri[0],fractTri[0]);
		copyVect(tri[1],mids[0]);
		copyVect(tri[2],mids[2]);
		addFract(tri,fractOffset*0.5,level-1,rgb,alpha,reflect,geom);
		copyVect(tri[0],fractTri[1]);
		copyVect(tri[1],mids[1]);
		copyVect(tri[2],mids[0]);
		addFract(tri,fractOffset*0.5,level-1,rgb,alpha,reflect,geom);
		copyVect(tri[0],fractTri[2]);
		copyVect(tri[1],mids[2]);
		copyVect(tri[2],mids[1]);
		addFract(tri,fractOffset*0.5,level-1,rgb,alpha,reflect,geom);
		copyVect(tri[0],mids[0]);
		copyVect(tri[1],mids[1]);
		copyVect(tri[2],mids[2]);
		addFract(tri,fractOffset*0.5,level-1,rgb,alpha,reflect,geom);
	}
}

void createGeometry(struct Geometry *geom) {
	int i,j;
	int opt;
	char fileName[40];
	FILE *file;
	double red,green,blue,rgb[3],rgb2[3];
	char buf[20];
	double sheet[2][2][3];
	double cylinder[2][21][3];
	double fractTri[3][3];
	double fractOffset;
	int numSides;
	int checker;
	double alpha,reflect;

	do {
		printf("1. Add Light\n2. Add Triangle\n3. Add Sphere\n4. Add fractal landscape\n5. Add sheet\n6. Add \"cylinder\"\n7. Write geometry file\n8. Print geometry\n9. Quit\n");
		scanf("%d",&opt);
		switch(opt) {
		case 1:
			printf("Input the x, y, and z for the light. ");
			readPoint(stdin,geom->light[geom->numLights].p);
			printf("Input red, green, and blue values (0.0-1.0)\n");
			scanf("%lf %lf %lf",&red,&green,&blue);
			geom->light[geom->numLights].col[0]=red;
			geom->light[geom->numLights].col[1]=green;
			geom->light[geom->numLights].col[2]=blue;
			geom->numLights++;
			break;
		case 2:
			printf("Input the x, y, and z for the three points. ");
			readTriangle(stdin,geom->tri[geom->numTris].p);
			printf("Input red, green, and blue values (0.0-1.0)\n");
			scanf("%lf %lf %lf",&red,&green,&blue);
			geom->tri[geom->numTris].col[0]=red;
			geom->tri[geom->numTris].col[1]=green;
			geom->tri[geom->numTris].col[2]=blue;
			printf("Enter fractional alpha and reflectivity: ");
			scanf("%lf %lf",&(geom->tri[geom->numTris].alpha),&(geom->tri[geom->numTris].reflect));
			geom->numTris++;
			break;
		case 3:
			printf("Input the x, y, and z for the sphere center. ");
			readPoint(stdin,geom->sphere[geom->numSpheres].c);
			printf("Input a radius. ");
			scanf("%lf",&(geom->sphere[geom->numSpheres].radius));
			printf("Input red, green, and blue values (0.0-1.0)\n");
			scanf("%lf %lf %lf",&red,&green,&blue);
			geom->sphere[geom->numSpheres].col[0]=red;
			geom->sphere[geom->numSpheres].col[1]=green;
			geom->sphere[geom->numSpheres].col[2]=blue;
			printf("Enter fractional alpha and reflectivity: ");
			scanf("%lf %lf",&(geom->sphere[geom->numSpheres].alpha),&(geom->sphere[geom->numSpheres].reflect));
			geom->numSpheres++;
			break;
		case 4:
			printf("Input the x, y, and z for the bounding triangle. ");
			readTriangle(stdin,fractTri);
			printf("Enter the original offset magnitude. ");
			scanf("%lf",&fractOffset);
			printf("How many levels do you want subdivide? ");
			scanf("%d",&numSides);
			printf("Input red, green, and blue values (0.0-1.0)\n");
			scanf("%lf %lf %lf",&red,&green,&blue);
			rgb[0]=red;
			rgb[1]=green;
			rgb[2]=blue;
			printf("Enter fractional alpha and reflectivity: ");
			scanf("%lf %lf",&alpha,&reflect);
			addFract(fractTri,fractOffset,numSides,rgb,alpha,reflect,geom);
			break;
		case 5:
			printf("Enter four corner points for the sheet in \"left hand\" order: ");
			readPoint(stdin,sheet[0][0]);
			readPoint(stdin,sheet[1][0]);
			readPoint(stdin,sheet[1][1]);
			readPoint(stdin,sheet[0][1]);
			printf("How many subsquares do you want? ");
			scanf("%d",&numSides);
			printf("Do you want it checkered (0=no, 1=yes)? ");
			scanf("%d",&checker);
			printf("Input red, green, and blue values (0.0-1.0)\n");
			scanf("%lf %lf %lf",&red,&green,&blue);
			rgb[0]=red;
			rgb[1]=green;
			rgb[2]=blue;
			if(checker) {
				printf("Input second red, green, and blue values (0.0-1.0)\n");
				scanf("%lf %lf %lf",&red,&green,&blue);
				rgb2[0]=red;
				rgb2[1]=green;
				rgb2[2]=blue;
			} else {
				rgb2[0]=rgb[0];
				rgb2[1]=rgb[1];
				rgb2[2]=rgb[2];
			}
			printf("Enter fractional alpha and reflectivity: ");
			scanf("%lf %lf",&alpha,&reflect);
			for(i=0; i<numSides; i++) {
				for(j=0; j<numSides; j++) {
					addSheetSquare(i,j,sheet,numSides,((i+j)%2==0)?rgb:rgb2,alpha,reflect,geom);
				}
			}
			break;
		case 6:
			printf("Enter the top center point: ");
			readPoint(stdin,cylinder[0][0]);
			printf("Enter the bottom center point: ");
			readPoint(stdin,cylinder[1][0]);
			printf("Enter one edge point on top rim: ");
			readPoint(stdin,cylinder[0][1]);
			printf("How many sides do you want (3-20)? ");
			scanf("%d",&numSides);
			printf("Input red, green, and blue values (0.0-1.0)\n");
			scanf("%lf %lf %lf",&red,&green,&blue);
			rgb[0]=red;
			rgb[1]=green;
			rgb[2]=blue;
			printf("Enter fractional alpha and reflectivity: ");
			scanf("%lf %lf",&alpha,&reflect);
			calcCylinderPoints(cylinder,numSides);
			for(i=1; i<=numSides; i++) {
				addCylinderPiece(i,cylinder,numSides,rgb,alpha,reflect,geom);
			}
			break;
		case 7:
			printf("Enter the name of the file. ");
			scanf("%s",fileName);
			file=fopen(fileName,"w");
			writeGeometry(file,geom);
			fclose(file);
			break;
		case 8:
			writeGeometry(stdout,geom);
			break;
		case 9:
			printf("Bye!\n");
			break;
		}
	} while(opt!=9);
}

int main(int argc,char **argv) {
	struct Geometry geom;

	createGeometry(&geom);
	return 0;
}