#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "coords.h"		// fvector arithmetic (c++ only)

////////////////////////// james, you'll need to change all these defines

#define SX 256	// size of subset of R3.... (ie, for 3d density array defined as density[SZ][SY][SX])
#define SY 256
#define SZ 256

#define BLOBSIZE 1	// size of each cube in space = scale factor
#define MAXBV 50000	// max number of vertices. DO NOT EXCEED THIS OR CORE DUMP
#define MAXBF 100000 // max number of faces. DO NOT EXCEED THIS OR CORE DUMP
#define MAXHASH 1024

#define potential(x,y,z) (x+y+z)	// define this to return sloath potential at integer coords x,y,z. returns a float. 
												// if the surface comes out "inside out", negate this and the threshold values you pass to "tracesurface"

void tracesurface(float thrsh); // this is THE MAIN ROUTINE! just call it with the float value of the isosurface
								// you want to trace. quite slow. (order of seconds on a p2)
								// see main() below for a demo of how to use it. the output format is as follows:
// you get points and faces (doh). the structures are defined just below here (BLOBV and BLOBF). the arrays are in
// global variables: bv, for the vertices, and bf. the number of points is numbv, number of faces (triangles) is 
// numbf. there's also normal information there (how kind). they might be pointing in the wrong direction though
// (swap the cross product over in tracesurface if so). its up to you to write this to a pov friendly format.
// see "output_to_pov" below for a demo of what you might do.
// oh and one note - DONT CHANGE global variables bv and bf (the array pointers) between calls to tracesurface,
// it expects the arrays to stay put.

#define error(x) {printf("OH DEAR too many points/tris. increase maxbv/maxbf in james.cpp\n");exit(x);} 

//////////////////////// end of nasty customise bit

struct BLOBV	// output vertices
{	
	FVector p,n;	// position and normal	
	BLOBV *next;	// for hash table	

};

struct BLOBF	// output faces
{
	BLOBV *a,*b,*c;	
};

float pval[8];
float tempcoord[3];
float thresh;

typedef struct
{
	SBYTE nedge;
	SBYTE edge[12];
	SBYTE ntri;
	SBYTE tri[8][3];
} SUBCUBE;
// big table!
#include "jamestab.cpp"

BLOBV *bv=NULL;
BLOBF *bf=NULL;
BLOBV *bvhash[MAXHASH];
int numbv,numbf;

int hashp(const FVector &pp)
{
	return (int(pp.X)^int(pp.Y*pp.Z))&(MAXHASH-1);
}

int addpoint(const FVector &pp)
{
	int h=hashp(pp);
	BLOBV *v=bvhash[h];
	while (v)
	{
		if (v->p==pp) return v-bv;
		v=v->next;
	}
	v=&bv[numbv++];
	if (numbv>=MAXBV) error(1);
	v->p=pp;
	v->next=bvhash[h];
	v->n=FVector(0,0,0);
	bvhash[h]=v;
	return v-bv;
}

void addface(int a,int b, int c)
{
	if (a==b || a==c || b==c) return;
	BLOBF *f=&bf[numbf++];
	if (numbf>=MAXBF) error(1);
	f->a=bv+a;
	f->b=bv+b;
	f->c=bv+c;
}



void calcblobpoint(SLONG i)
{
	switch (i)
	{
	case 0:

		if (pval[0]!=pval[4]) tempcoord[0]=(BLOBSIZE*pval[0])/(pval[0]-pval[4]);
		tempcoord[1]=0;
		tempcoord[2]=0;
		break;
	case 1:
		if (pval[1]!=pval[5]) tempcoord[0]=(BLOBSIZE*pval[1])/(pval[1]-pval[5]);
		tempcoord[1]=0;
		tempcoord[2]=BLOBSIZE;
		break;
	case 2:
		if (pval[2]!=pval[6]) tempcoord[0]=(BLOBSIZE*pval[2])/(pval[2]-pval[6]);
		tempcoord[1]=BLOBSIZE;
		tempcoord[2]=0;
		break;
	case 3:
		if (pval[7]!=pval[3]) tempcoord[0]=(BLOBSIZE*pval[3])/(pval[3]-pval[7]);
		tempcoord[1]=BLOBSIZE;
		tempcoord[2]=BLOBSIZE;
		break;
	case 4:
		tempcoord[0]=0;
		if (pval[0]!=pval[2]) tempcoord[1]=(BLOBSIZE*pval[0])/(pval[0]-pval[2]);
		tempcoord[2]=0;
		break;
	case 5:
		tempcoord[0]=0;
		if (pval[1]!=pval[3]) tempcoord[1]=(BLOBSIZE*pval[1])/(pval[1]-pval[3]);
		tempcoord[2]=BLOBSIZE;
		break;
	case 6:
		tempcoord[0]=BLOBSIZE;
		if (pval[6]!=pval[4]) tempcoord[1]=(BLOBSIZE*pval[4])/(pval[4]-pval[6]);
		tempcoord[2]=0;
		break;
	case 7:
		tempcoord[0]=BLOBSIZE;
		if (pval[5]!=pval[7]) tempcoord[1]=(BLOBSIZE*pval[5])/(pval[5]-pval[7]);
		tempcoord[2]=BLOBSIZE;
		break;
	case 8:
		tempcoord[0]=0;
		tempcoord[1]=0;
		if (pval[0]!=pval[1]) tempcoord[2]=(BLOBSIZE*pval[0])/(pval[0]-pval[1]);
		break;
	case 9:
		tempcoord[0]=0;
		tempcoord[1]=BLOBSIZE;
		if (pval[2]!=pval[3]) tempcoord[2]=(BLOBSIZE*pval[2])/(pval[2]-pval[3]);
		break;
	case 10:
		tempcoord[0]=BLOBSIZE;
		tempcoord[1]=0;
		if (pval[5]!=pval[4]) tempcoord[2]=(BLOBSIZE*pval[4])/(pval[4]-pval[5]);
		break;
	case 11:
		tempcoord[0]=BLOBSIZE;
		tempcoord[1]=BLOBSIZE;
		if (pval[6]!=pval[7]) tempcoord[2]=(BLOBSIZE*pval[6])/(pval[6]-pval[7]);
		break;
	}

}



void addcubetris(SLONG x,SLONG y,SLONG z)
{
	SLONG in;
	in  = ((pval[0]=potential(x  ,y  ,z  )-thresh)>0?  1:0);
	in |= ((pval[1]=potential(x  ,y  ,z+1)-thresh)>0?  2:0);
	in |= ((pval[2]=potential(x  ,y+1,z  )-thresh)>0?  4:0);
	in |= ((pval[3]=potential(x  ,y+1,z+1)-thresh)>0?  8:0);
	in |= ((pval[4]=potential(x+1,y  ,z  )-thresh)>0? 16:0);
	in |= ((pval[5]=potential(x+1,y  ,z+1)-thresh)>0? 32:0);
	in |= ((pval[6]=potential(x+1,y+1,z  )-thresh)>0? 64:0);
	in |= ((pval[7]=potential(x+1,y+1,z+1)-thresh)>0?128:0);

	SLONG i,j;
	int idx[12];
	for (i=0;i<cubeinfo[in].nedge;i++)
	{
		calcblobpoint(j=cubeinfo[in].edge[i]);
		if (tempcoord[0]<0) tempcoord[0]=0;
		if (tempcoord[0]>BLOBSIZE) tempcoord[0]=BLOBSIZE;
		if (tempcoord[1]<0) tempcoord[1]=0;
		if (tempcoord[1]>BLOBSIZE) tempcoord[1]=BLOBSIZE;
		if (tempcoord[2]<0) tempcoord[2]=0;
		if (tempcoord[2]>BLOBSIZE) tempcoord[2]=BLOBSIZE;

		idx[j]=addpoint(FVector(tempcoord[0]+(x-SX/2)*BLOBSIZE,tempcoord[1]+(y-SY/2)*BLOBSIZE,tempcoord[2]+(z-SZ/2)*BLOBSIZE));
	}
	for (i=0;i<cubeinfo[in].ntri;i++)
	{
		addface(idx[cubeinfo[in].tri[i][0]],idx[cubeinfo[in].tri[i][1]],idx[cubeinfo[in].tri[i][2]]);
	}
}

void tracesurface(float thrsh)
{
	thresh=thrsh;
	int x,y,z;
	numbf=numbv=0;
	memset(bvhash,0,sizeof(bvhash));
	if (!bv) bv=new BLOBV[MAXBV];
	if (!bf) bf=new BLOBF[MAXBF];
	printf("tracing surface %0.3f",thresh);
	
	for (z=0;z<SZ-1;z++)
	{
		for (y=0;y<SY-1;y++)
		{
			for (x=0;x<SX-1;x++)
			{
				addcubetris(x,y,z);
			}
		}
		if ((z&7)==0) {printf(".");flushall();}
	}

	// normals - works by averaging the
	// face normals shared by a point. to give a nice curved shading.
	
	BLOBF *f=bf;	// list of faces. each face has 3 indices
	int c1;
	for (c1=0;c1<numbf;c1++,f++)
	{
		FVector a=f->a->p-f->b->p;a.Normalise();
		FVector b=f->c->p-f->b->p;b.Normalise();
		FVector n=a^b;n.Normalise();
		f->a->n+=n;
		f->b->n+=n;
		f->c->n+=n;
	}
	BLOBV *v=bv;
	for (c1=0;c1<numbv;c1++,v++) v->n.Normalise();

	printf("found surface - %d points, %d faces\n",numbv,numbf);
}



//////////////////// DEMO of how to use it (I hope)

FILE *povf;

void output_to_pov()
{
	// go through each triangle in turn
	printf("writing to pov file\n");
	BLOBF *f=bf;
	for (int c1=0;c1<numbf;c1++)
	{		
		fprintf(povf,"triangle: (3 points, then 3 normals) (%f,%f,%f),(%f,%f,%f),(%f,%f,%f) normals: (%f,%f,%f),(%f,%f,%f),(%f,%f,%f)\n",
			f->a->p.x,f->a->p.y,f->a->p.z,
			f->b->p.x,f->b->p.y,f->b->p.z,
			f->c->p.x,f->c->p.y,f->c->p.z,
			f->a->n.x,f->a->n.y,f->a->n.z,
			f->b->n.x,f->b->n.y,f->b->n.z,
			f->c->n.x,f->c->n.y,f->c->n.z);

		f++;
	}
}

int main(int argc, char *argv[])
{
	// this example traces 3 isosurfaces at value 10, 20 and 30 and "outputs" them
	// set up sloath data here
	// .....
	// open povray output file here. no idea about format...
	// .......
	povf=fopen("out.pov","w");
	
	// now do the work like this:
	tracesurface(10);	// trace surface with value 10
	output_to_pov();	// output to povray file here...
	tracesurface(20);	// and 20
	output_to_pov();	// output to povray file here...
	tracesurface(30);	// and 30
	output_to_pov();	// output to povray file here...
	
	// close povray file
	// .......
	fclose(povf);
	// the end.
	delete [] bv;
	delete [] bf;
}
