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

float _eye[3];
float _topleft[3];
float _u[3];
float _v[3];
int _num_tr;

typedef struct {
	float p[3][3];
	float plane_eq[4];
	unsigned char intensity;
} Triangle;
Triangle *_tr = NULL;

void cross_product(float A[3], float B[3], float C[3])
{
	// C shouldn't be A or B
	C[0] = A[1]*B[2]-A[2]*B[1];
	C[1] = A[2]*B[0]-A[0]*B[2];
	C[2] = A[0]*B[1]-A[1]*B[0];
}

float dot_product(float A[3], float B[3])
{
	return A[0]*B[0] + A[1]*B[1] + A[2]*B[2];
}

// find the signed volume of tetrahedron
float tetrahedron_volume(float A[3], float B[3], float C[3], float D[3])
{
float AB[3],AC[3],AD[3], ADxAC[3];

	AB[0] = B[0]-A[0];	AC[0] = C[0]-A[0];	AD[0] = D[0]-A[0];
	AB[1] = B[1]-A[1];	AC[1] = C[1]-A[1];	AD[1] = D[1]-A[1];
	AB[2] = B[2]-A[2];	AC[2] = C[2]-A[2];	AD[2] = D[2]-A[2];
	cross_product(AD,AC,ADxAC);

	return dot_product(AB,ADxAC)/6;
}

// inside test
float tetrahedron_inside_test(float A[3], float B[3], float C[3], float E[3], float P[3])
{
	float a1 = tetrahedron_volume(E,P,C,A);
	float a2 = tetrahedron_volume(E,P,B,C);
	float a3 = tetrahedron_volume(E,P,A,B);
	float a4 = tetrahedron_volume(A,P,C,B);

	return (a1*a2>0) && (a2*a3>0) && (a3*a4>0);
}

void normalize_vector(float v[3])
{
float s = (float)sqrt(v[0]*v[0] + v[1]*v[1] + v[2]*v[2]);

	v[0] /= s;	v[1] /= s;	v[2] /= s;
}

char render_pixel(float P[3])
{
float t, dist2smallest=1e10;
float ray[3] = {P[0] - _eye[0],  P[1] - _eye[1],  P[2] - _eye[2]};
char  color = 0;

	for(int i=0; i<_num_tr; i++) {
		if(tetrahedron_inside_test(_tr[i].p[0],_tr[i].p[1],_tr[i].p[2],_eye,P)) {
			// find the t value, t value of 1 represents the depth |ray vector|
			t = -(_tr[i].plane_eq[3] + dot_product(_eye,_tr[i].plane_eq))
				  / dot_product(ray,_tr[i].plane_eq);
			if(t>0 && t<dist2smallest) {
				dist2smallest = t;
				color = _tr[i].intensity;
			}
		}
	}
	return color;
}

void main(int argc, char *argv[])
{
char image[256][256];
float a[3],b[3];
int i,j;

	if(argc!=2) return;

	// read file
	FILE *fp = fopen(argv[1],"rt");
	fscanf(fp,"%f %f %f",&_eye[0],&_eye[1],&_eye[2]);
	fscanf(fp,"%f %f %f",&_topleft[0],&_topleft[1],&_topleft[2]);
	fscanf(fp,"%f %f %f",&_u[0],&_u[1],&_u[2]);
	fscanf(fp,"%f %f %f",&_v[0],&_v[1],&_v[2]);
	fscanf(fp,"%i",&_num_tr);
	_tr = new Triangle[_num_tr];
	for(i=0;i<_num_tr; i++) {
		fscanf(fp,"%f %f %f   %f %f %f   %f %f %f",
			&_tr[i].p[0][0],&_tr[i].p[0][1],&_tr[i].p[0][2],
			&_tr[i].p[1][0],&_tr[i].p[1][1],&_tr[i].p[1][2],
			&_tr[i].p[2][0],&_tr[i].p[2][1],&_tr[i].p[2][2]);
	}
	fclose(fp);

	// calculate the intensity of the triangle, equation of plane
	float light[3] = { -1/(float)sqrt(3), -1/(float)sqrt(3), -1/(float)sqrt(3)};
	for(i=0; i<_num_tr; i++) {
		a[0] = _tr[i].p[1][0] - _tr[i].p[0][0];		b[0] = _tr[i].p[2][0] - _tr[i].p[0][0];
		a[1] = _tr[i].p[1][1] - _tr[i].p[0][1];		b[1] = _tr[i].p[2][1] - _tr[i].p[0][1];
		a[2] = _tr[i].p[1][2] - _tr[i].p[0][2];		b[2] = _tr[i].p[2][2] - _tr[i].p[0][2];
		cross_product(a,b,_tr[i].plane_eq);
		normalize_vector(_tr[i].plane_eq);
		_tr[i].plane_eq[3] = - _tr[i].p[0][0]*_tr[i].plane_eq[0]
							 - _tr[i].p[0][1]*_tr[i].plane_eq[1]
							 - _tr[i].p[0][2]*_tr[i].plane_eq[2];
		_tr[i].intensity = (char)(-dot_product(_tr[i].plane_eq,light)*127+128);
	}

	// render all pixels
	float P[3];
	for(i=0; i<256; i++) for(j=0; j<256; j++) {
		P[0] = _topleft[0] + i*_u[0] + j*_v[0];
		P[1] = _topleft[1] + i*_u[1] + j*_v[1];
		P[2] = _topleft[2] + i*_u[2] + j*_v[2];
		image[j][i] = render_pixel(P);
	}

	// save image
	fp = fopen("out.raw","wb");
	fwrite(image,1,256*256,fp);
	fclose(fp);

	// dump ??
//	for(i=0; i<256; i++) for(j=0; j<256; j++) {
//		printf("%02X ",image[j][i]);
//	}

	// cleanup the memory
	if(_tr!=NULL) delete _tr;
}