#include<sys/types.h>
#include<sys/stat.h>
#include<fcntl.h>
#include<unistd.h>
#include<stdlib.h>
#include<stdio.h>
#include<malloc.h>
#include<string.h>

#define Epsilon 0.001

typedef struct {
	float x,y,z;
	float r,g,b,a;
} Vertex;

typedef struct {
	int p1, p2, p3;
} Triangle ;

typedef struct ls {
	 int vert;
	 struct ls *next;
} listsom ;

typedef struct {
	 listsom *points;
	 int newpoint;
} grid ;

grid *Grille;
Vertex *v, *sv;
Triangle *t, *st;
int nbv,nbt, nbsv, nbst;
float *W;
grid **Cell;
char basefilename[256];

void AjouteAListe( grid *g, int i )
{
	listsom *ls;

	ls = (listsom *)malloc( sizeof(listsom) );
	if ( ls == NULL ) fprintf(stderr,"Erreur malloc pour une liste\n");

	ls->next = g->points;
	ls->vert = i;
	g->points = ls;
}

int ReadGeomFile( char *filename , int *nbv, int *nbt, Vertex **V, Triangle **T )
{
	FILE *fp;
	char line[512];
	char *s;
	int lineno = 1;
	int p1,p2,p3,p4,tri,nbpol,nbpts,nbtri;
	int dummy,tt,lus;
	int i;
	float x,y,z;
	float r,g,b,alpha;
	Vertex *v; Triangle *t;

	if ( (fp = fopen( filename, "r" )) == NULL ) {
		fprintf(stderr,"Could not open file %s.\n",filename);
		return(0);
	}

	if ( (s = fgets( line , 512 , fp )) == NULL ) {
		fprintf(stderr,"Error in file %s near line %d\n",filename,lineno);
		return(0);
	}

	while ( sscanf( line , "%d %d %d" , &nbpts, &nbpol, &dummy ) != 3 ) {
		lineno++;
		if ( (s = fgets( line , 512 , fp )) == NULL ) {
			fprintf(stderr,"Error: wrong format in file %s\n",filename);
			return(0);
		}
	}
	
	v = *V = (Vertex *)malloc(nbpts * sizeof(Vertex) );

	for ( i=0 ; i < nbpts ; ++i ) {
		if ( (s = fgets( line , 512 , fp )) == NULL ) {
			fprintf(stderr,"Error: could not read point %d in file %s\n",
				i+1, filename);
			return(0);
		}

	 	if ( (lus = sscanf( line, "%f %f %f %f %f %f %f", &x, &y, &z, &r, &g, &b, &alpha)) < 3 ) {
			fprintf(stderr,"Error: wrong format for point %d in file %s\n",
				i+1, filename);
			return(0);
		}

		v[i].x = x;
		v[i].y = y;
		v[i].z = z;
		if ( lus > 3 ) {
			v[i].r = r;
			v[i].g = g;
			v[i].b = b;
			v[i].a = alpha;
		} else {
			v[i].r = v[i].g = v[i].b = 0.9;
		}
	}


	t = *T = (Triangle *)malloc(2*nbpol * sizeof(Triangle) );
	nbtri = 0;

	for ( i=0 ; i < nbpol ; ++i ) {
		if ( (s = fgets( line , 512 , fp )) == NULL ) {
			fprintf(stderr,"Error: could not read pol %d in file %s\n",
				i+1, filename);
			return(0);
		}

	 	tt = sscanf( line, "%d %d %d %d %d", &tri, &p1, &p2, &p3, &p4 );

		if ( tt < 4 ) {
			fprintf(stderr,"Error: wrong format for pol %d in file %s\n",
				i+1, filename);
			return(0);
		}
		
		if ( tri == 3 ) {
			t[nbtri].p1 = p1;
			t[nbtri].p2 = p2;
			t[nbtri].p3 = p3;
			nbtri ++;
		}  else
		if ( tri == 4 ) {

			/* Creer DEUX triangles a partir des points p1--p4 */
			/* et les ranger dans le tableau aux emplacements nbtri et nbtri+1 */

			t[nbtri].p1 = p1;
			t[nbtri].p2 = p2;
			t[nbtri].p3 = p3;

			t[nbtri+1].p1 = p1;
			t[nbtri+1].p2 = p3;
			t[nbtri+1].p3 = p4;

			nbtri += 2;

		} else {
			fprintf(stderr,"Warning: ignoring polygon with %d vertices\n", tri);
		}
	}
	fclose(fp);

	printf("Successfully read file %s with %d vertices and %d triangles\n",
		filename, nbpts, nbtri);

	*nbv = nbpts;
	*nbt = nbtri;
	return(1);
}

int WriteGeomFile( char *filename , int nbv, int nbt, Vertex *v, Triangle *t )
{
	FILE *fp;
	int i;

	if ( (fp = fopen( filename, "w" )) == NULL ) {
		fprintf(stderr,"Could not open file %s for writing.\n",filename);
		return(0);
	}

	fprintf(fp,"COFF\n%d %d %d\n",nbv,nbt,0);

	for ( i=0 ; i < nbv ; ++i ) {
		fprintf(fp,"\t%9.6f\t%9.6f\t%9.6f\t\t%5.3f\t%5.3f\t%5.3f\t%5.3f\n",
			v[i].x, v[i].y, v[i].z,
			v[i].r, v[i].g, v[i].b, v[i].a
			);
	}

	for ( i=0 ; i < nbt ; ++i ) {
		fprintf(fp,"3\t%3d\t%3d\t%3d\n",
/*			t[i].p1-1, t[i].p2-1, t[i].p3-1 */
			t[i].p1, t[i].p2, t[i].p3 
			);
	}
	fclose(fp);
	printf("Successfully wrote file %s with %d vertices and %d triangles\n",
		filename, nbv, nbt);
	return 1;
}

void Simplify( int numcells )
{
	float min_X = 1e6;
	float max_X = -1e6;
	float min_Y = 1e6;
	float max_Y = -1e6;
	float min_Z = 1e6;
	float max_Z = -1e6;
	char filename[256];

	int nx,ny,nz;
	int i;
	grid *g;
	float weight;
	float X,Y,Z;
	float R,G,B,A;
	listsom *l;

	/* D'abord on regarde l'etendue du modele en 3D */
	for ( i=0 ; i < nbv ; ++i ) {
		if ( v[i].x < min_X ) min_X = v[i].x;
		if ( v[i].x > max_X ) max_X = v[i].x;
		if ( v[i].y < min_Y ) min_Y = v[i].y;
		if ( v[i].y > max_Y ) max_Y = v[i].y;
		if ( v[i].z < min_Z ) min_Z = v[i].z;
		if ( v[i].z > max_Z ) max_Z = v[i].z;
	}
	min_X -= Epsilon;
	min_Y -= Epsilon;
	min_Z -= Epsilon;
	max_X += Epsilon;
	max_Y += Epsilon;
	max_Z += Epsilon;

/*	printf("Bounding box %8.4f %8.4f %8.4f ->  %8.4f %8.4f %8.4f\n",
		min_X, min_Y, min_Z, max_X, max_Y, max_Z ); */

	/* Puis on cree la grille... */
	Grille = (grid *)malloc( numcells*numcells*numcells*sizeof(grid));
	if ( Grille == NULL ) {
		fprintf(stderr,"Erreur de malloc\n");
		exit(-1);
	}
	/* et on l'initialise... */
	for ( i=0 ; i < numcells*numcells*numcells ; ++i ) {
		Grille[i].newpoint = 0;
		Grille[i].points = NULL;
	}

	/* Et on alloue les tableaux annexes */
	W = (float *)malloc( nbv*sizeof(float) );
	Cell = (grid **)malloc( nbv*sizeof(grid *) );
	if ( (!W) || (!Cell) ) {
		fprintf(stderr,"Erreur de malloc\n");
		exit(-1);
	}

	/* Maintenant, on remplit le tableau Cell et les listes de sommets de
	   chaque case de la grille */
	for ( i=0 ; i < nbv ; ++i ) {
		/* on obtient les coordonnees de notre case dans la grille par
		   division entiere : */
		nx = (int)(numcells*( v[i].x - min_X ) / (max_X - min_X));
		ny = (int)(numcells*( v[i].y - min_Y ) / (max_Y - min_Y));
		nz = (int)(numcells*( v[i].z - min_Z ) / (max_Z - min_Z));

		g = Grille + nx + numcells*( ny + numcells*nz );
		Cell[i] = g;

		AjouteAListe(g,i);

/*		printf("Sommet %f %f %f -> %d %d %d\n",v[i].x, v[i].y, v[i].z,nx,ny,nz); */
		/* on calcule aussi le poids de ce sommet */
		/* Pour l'instant on ne reflechit pas, tous les poids sont egaux */
		W[i] = 1.0;
	}

	/* Pour chaque case de la grille, on elit un representant */
	g = Grille; 
	nbsv = 0;
	sv = (Vertex *)malloc( nbv * sizeof( Vertex ) );

	for ( i=0 ; i < numcells*numcells*numcells ; ++i ) {
		/* Calcule le barycentre des points de cette case */
		weight = 0;
/*		printf("Grille %d: sommets\n",i); */
		X=Y=Z=0;
		R=G=B=A=0;
		for ( l = g->points ; l ; l = l->next ) {
			weight += W[ l->vert ];
			X +=  W[ l->vert ]* v[ l->vert ].x;
			Y +=  W[ l->vert ]* v[ l->vert ].y;
			Z +=  W[ l->vert ]* v[ l->vert ].z;
			R +=  W[ l->vert ]* v[ l->vert ].r;
			G +=  W[ l->vert ]* v[ l->vert ].g;
			B +=  W[ l->vert ]* v[ l->vert ].b;
			A +=  W[ l->vert ]* v[ l->vert ].a;
/*			printf("\t %8.4f %8.4f %8.4f\n",v[ l->vert ].x, v[ l->vert ].y, v[ l->vert ].z); */
		}
		if ( weight > 0 ) {
			/* Le poids total est non nul ssi il y avait au moins un point */
			X /= weight; Y /= weight; Z /= weight;
			R /= weight; G /= weight; B /= weight; A /= weight;
			sv[nbsv].x = X;
			sv[nbsv].y = Y;
			sv[nbsv].z = Z;
			sv[nbsv].r = R;
			sv[nbsv].g = G;
			sv[nbsv].b = B;
			sv[nbsv].a = A;
/*			printf("\t nouveau point %8.4f %8.4f %8.4f (%d)\n",X,Y,Z,nbsv); */
			g->newpoint = nbsv;
			nbsv++;
		}
		g++;
	}

	/* OK maintenant on cree les nouveaux triangles... */

	nbst = 0;
	st = (Triangle *)malloc( nbt * sizeof( Triangle ) );
	for ( i=0 ; i < nbt ; ++i ) {
		/* On regarde si les trois points sont dans des cellules distinctes */
		if ( (Cell[t[i].p1] != Cell[t[i].p2]) && 
			 (Cell[t[i].p1] != Cell[t[i].p3]) && (Cell[t[i].p3] != Cell[t[i].p2]) ) {

/*			 printf("%8.4f %8.4f %8.4f -> %8.4f %8.4f %8.4f\n",
				v[ t[i].p1 ].x,  v[ t[i].p1 ].y,  v[ t[i].p1 ].z,
				sv[ Cell[ t[i].p1 ]->newpoint ].x,
				sv[ Cell[ t[i].p1 ]->newpoint ].y,
				sv[ Cell[ t[i].p1 ]->newpoint ].z ); */

			st[nbst].p1 = Cell[ t[i].p1 ]->newpoint;
			st[nbst].p2 = Cell[ t[i].p2 ]->newpoint;
			st[nbst].p3 = Cell[ t[i].p3 ]->newpoint;
			nbst++;
		}
	}

	sprintf(filename,"%s.S%d", basefilename, numcells);
	WriteGeomFile(filename, nbsv, nbst, sv, st );
}

int main( int argc, char **argv )
{
	int i;
	if ( argc > 1 ) {
		/* on elimine la partie inutile du nom de fichier */
		for ( i=strlen(argv[1])-1; i && argv[1][i] != '/' ; i-- );
		if ( argv[1][i] == '/' ) i++;
		strcpy(basefilename, &(argv[1][i]));
		printf("Extracted base filename %s\n",basefilename);

		if ( ! ReadGeomFile(argv[1],&nbv,&nbt,&v,&t) )
			exit(-1);
		Simplify( argc > 2 ? atoi(argv[2]) : 10 );
	} else
		printf("Usage: %s Geomfile [grid_size]\n",argv[0]);
}