#include <stdio.h>
#include <math.h>
#include <exec/types.h>
#include <graphics/gfx.h>
#include <graphics/rastport.h>
#include <graphics/copper.h>
#include <graphics/view.h>
#include <graphics/gels.h>
#include <graphics/regions.h>
#include <graphics/clip.h>
#include <exec/exec.h>
#include <graphics/text.h>
#include <graphics/gfxbase.h>
#include <hardware/dmabits.h>
#include <hardware/custom.h>
#include <hardware/blit.h>
#include <devices/keymap.h>

#define CARRE(x) ((x) * (x))

#define DB1 0
#define DB2 1
#define PLAN 1
#define COLOR 2
#define HORIZONT 640
#define VERT 512
#define MODE HIRES|LACE

UWORD dbuff[2] = {DB1,DB2};

int bascule;
int i;

struct View view; 
struct ViewPort viewport;

struct RasInfo rasinfo[2];

struct BitMap bitmap[2];
struct RastPort rastport[2];
struct GfxBase *GfxBase;

struct View *ecran_sauvegarde;

UWORD colortable[COLOR] = { 0x000, 0x77F };

char *bouton_gauche = (char *) 0xbfe001;

float sinus[360],cosinus[360]; 

double r;

struct cprlist *LOF[2];
struct cprlist *SHF[2];

void init(),ouvreecran(),chargepoint(),chargeligne(),libere(),cree_ecran();
void init_table(),elimine(),trace(),chargeface(),echange();

double point3D[100][3];
double sommet[100][3];
int ligne[100][2];
struct { int A,B,C;
		 double a,b,c,h;
} face[100];

int nbrpoint,nbrligne,nbrface ,nbrligneE;

double peps = 1e-5, meps = -1e-5, unmeps = 1-1e-5, unpeps = 1+1e-5;
int alpha = 0 ,teta = 0;
int va,vt;
char nomligne[100],nompoint[100],nomface[100];



/* --------------------*/
/* Programme principal */
/* --------------------*/

void main()
{
	init();
	printf("nom du fichier de point : ");
	scanf("%s",nompoint);
	printf("nom du fichier de ligne : ");
	scanf("%s",nomligne);
	printf("nom du fichier de face : ");
	scanf("%s",nomface);
	printf("valeur de la vitesse V(alpha) : ");
	scanf("%d",&va);
	printf("valeur de la vitesse V(teta) : ");
	scanf("%d",&vt);
	printf("distance r : ");
	scanf("%lf",&r);
	chargepoint();
	chargeligne();
	chargeface();
	init_table();
	ouvreecran();
	SetAPen(&rastport[dbuff[DB1]],1);
	SetAPen(&rastport[dbuff[DB2]],1);
	bascule = 0;
	while(!((*bouton_gauche & 0x40) - 64)) {
		SetRast(&rastport[dbuff[!bascule]],0);

		alpha += va;
		if(alpha >= 360)
			alpha = 0;

		teta += vt;
		if(teta >= 360)
			teta = 0;

		cree_ecran(&view,bascule);
		elimine();
		bascule ^= 1;
	}
	LoadView(ecran_sauvegarde); /* récupération de l'écran sauvegardé */
	libere();
}

/* ------------------------------- */
/* init() Ouvre la bibliothèque et */
/* sauvegarde l'ancine écran       */
/* ------------------------------- */

void init()
{
	if((GfxBase = (struct GfxBase *)OpenLibrary("graphics.library",0)) == NULL)
		exit(1);

	ecran_sauvegarde = GfxBase->ActiView; /* sauvegarde de l'écran courant */
}

/* ---------------------------------------------------- */
/* chargepoints() charge la table des points sur disque */
/* ---------------------------------------------------- */

void chargepoint()
{
	FILE *fich,*fopen();
	
	fich = fopen(nompoint, "r");
	if(fich == NULL) {
		printf("Problème d'ouverture du ficher Point.dat\n");
		exit(FALSE);
	}
	else
	{
		fscanf(fich,"%d",&nbrpoint);
		printf("Chargement de la structure POINT en cours...\n");

		for(i = 0; i < nbrpoint;i++)
		{
			fscanf(fich, "%lf %lf %lf",&point3D[i][0],&point3D[i][1],&point3D[i][2]);
		}
	}
	fclose(fich);
}

/* --------------------------------------------- */
/* chargeligne() charge les lignes sur le disque */
/* --------------------------------------------- */

void chargeligne()
{
	FILE *fich,*fopen();
	int origine,extreme;

	fich = fopen(nomligne,"r");
	if(fich == NULL) {
		printf("Problème d'ouverture du fichier ligne.dat\n");
		exit(FALSE);
	}
	else
	{
		fscanf(fich,"%d",&nbrligne);
		printf("Chargement de la structure LIGNE en cours...\n");

		for(i = 0; i < nbrligne;i++) {
			fscanf(fich,"%d %d",&origine,&extreme);
			ligne[i][0] = origine;
			ligne[i][1] = extreme;
		}
	}
	fclose(fich);
}

/* ------------------------------------------- */
/* chargeface() charge les faces sur le disque */
/* ------------------------------------------- */

void chargeface()
{
	FILE *fich,*fopen();
	int i;
	int A, B, C;
	fich = fopen(nomface,"r");
	if(fich == NULL) {
		printf("Problème d'ouverture du fichier Face.dat\n");
		exit(FALSE);
	}
	else
	{
		fscanf(fich,"%d",&nbrface);
		printf("Chargement de la structure FACE en cours...\n");
		printf("nbrface = %d\n",nbrface);
		for(i = 0; i < nbrface;i++)	{
			fscanf(fich,"%d %d %d",&A,&B,&C);
			face[i].A = A;
			face[i].B = B;
			face[i].C = C;
		}
	}
	fclose(fich);
}

/* ---------------------------------------------------------- */
/* ouvreecran() initialise le double buffer et ouvre un écran */
/* ---------------------------------------------------------- */

void ouvreecran()
{
	InitView(&view);
	InitVPort(&viewport); 
	view.ViewPort = &viewport;
	
	view.Modes = MODE;
	
	InitBitMap(&bitmap[DB1],PLAN,HORIZONT,VERT);
	InitBitMap(&bitmap[DB2],PLAN,HORIZONT,VERT);

	for(i=0; i<2; i++)
	{
		rasinfo[i].BitMap = &bitmap[i];
		rasinfo[i].RxOffset = 0;
		rasinfo[i].RyOffset = 0;
		rasinfo[i].Next = NULL;
	}
	
	viewport.DxOffset = 0;
	viewport.DyOffset = 0;
	viewport.DWidth = HORIZONT;
	viewport.DHeight = VERT;
	viewport.RasInfo = &rasinfo[DB1];
	viewport.Modes = MODE;
	viewport.Next = NULL;
	viewport.ColorMap = (struct ColorMap *)GetColorMap(COLOR);
	LoadRGB4(&viewport,&colortable[0],COLOR);

	for( i=0;i<PLAN;i++) {
		if ((bitmap[DB1].Planes[i] = (PLANEPTR)AllocRaster(HORIZONT,VERT)) == NULL)
			exit(2);
		BltClear ((UBYTE *)bitmap[DB1].Planes[i],RASSIZE(HORIZONT,VERT),0);
	}

	for( i=0;i<PLAN;i++) {
		if ((bitmap[DB2].Planes[i] = (PLANEPTR)AllocRaster(HORIZONT,VERT)) == NULL)
			exit(2);
		BltClear ((UBYTE *)bitmap[DB2].Planes[i],RASSIZE(HORIZONT,VERT),0);
	}

	for( i=0;i<2;i++)
	{
		InitRastPort(&rastport[i]);
		rastport[i].BitMap = &bitmap[i];
	}
	
	MakeVPort(&view, &viewport);
	
	MrgCop(&view);
	
	LOF[DB1] = view.LOFCprList;
	SHF[DB1] = view.SHFCprList;

	viewport.RasInfo = &rasinfo[DB2];
	rasinfo[DB2].Next = &rasinfo[DB1];

	view.LOFCprList = NULL;
	view.SHFCprList = NULL;
	
	MakeVPort(&view, &viewport);
	
	MrgCop(&view);
	
	LOF[DB2] = view.LOFCprList;
	SHF[DB2] = view.SHFCprList;
}

/* ----------------------------------- */
/* libere() libère la mémoire utilisée */
/* ----------------------------------- */

void libere()
{
	for(i = 0;i<PLAN;i++)
		
		FreeRaster(bitmap[DB1].Planes[i],HORIZONT,VERT);

	for(i = 0;i<PLAN;i++)
		FreeRaster(bitmap[DB2].Planes[i],HORIZONT,VERT);

	FreeColorMap(viewport.ColorMap);
	
	FreeVPortCopLists(&viewport);
	
	FreeCprList(LOF[DB1]);
	FreeCprList(SHF[DB1]);
	FreeCprList(LOF[DB2]);
	FreeCprList(SHF[DB2]);
	
	CloseLibrary(GfxBase);
}

/* ---------------------------- */
/* cree_ecran() affiche l'écran */
/* ---------------------------- */

void cree_ecran(view,bascule)
struct View *view;
int bascule;
{
	view->LOFCprList = LOF[bascule];
	view->SHFCprList = SHF[bascule];
	LoadView(view);
	WaitTOF();
}

/* ----------------------------------------------------- */
/* init_table() initialise la table des sinus et cosinus */
/* ----------------------------------------------------- */

void init_table()
{
	double rad;

	rad = (double)(2.*PI/360.);
	for (i=0;i<360;i++)
	{
		sinus[i]   = sin((rad*(double)i));
		cosinus[i] = cos((rad*(double)i));
	}
}

/* ------------------------------------------ */
/* Elimination des parties cachées et traçage */
/* ------------------------------------------ */

void elimine()
{
	double v11,v12,v13,v21,v22,v23,v32,v33;
	double unit;

	/* Coefficients de la matrice de transformation perspective */
	
	v11 = -sinus[alpha];
	v12 = -cosinus[teta]*cosinus[alpha];
	v13 = -sinus[teta]*cosinus[alpha];
	v21 = cosinus[alpha];
	v22 = -cosinus[teta]*sinus[teta];
	v23 = -sinus[teta]*sinus[alpha];
	v32 = sinus[teta];
	v33 = -cosinus[teta];
	
	/* Passage des points dans le repère observateur */
	
	for(i=0; i<nbrpoint; i++) {
		sommet[i][0]=v11*point3D[i][0]+v21*point3D[i][1];
		sommet[i][1]=v12*point3D[i][0]+v22*point3D[i][1]+v32*point3D[i][2];
		sommet[i][2]=v13*point3D[i][0]+v23*point3D[i][1]+v33*point3D[i][2]+r;
	}
	
	/* Initialisation des coefficients des faces */

	for(i=0; i<nbrface; i++) {
		face[i].a= sommet[face[i].A][1]*(sommet[face[i].B][2] - sommet[face[i].C][2])
			-sommet[face[i].B][1]*(sommet[face[i].A][2] - sommet[face[i].C][2])
			+sommet[face[i].C][1]*(sommet[face[i].A][2] - sommet[face[i].B][2]);

		face[i].b= -(sommet[face[i].A][0]*(sommet[face[i].B][2] - sommet[face[i].C][2])
			- sommet[face[i].B][0]*(sommet[face[i].A][2] - sommet[face[i].C][2])
			+ sommet[face[i].C][0]*(sommet[face[i].A][2] - sommet[face[i].B][2]));

		face[i].c= sommet[face[i].A][0]*(sommet[face[i].B][1] - sommet[face[i].C][1])
			-sommet[face[i].B][0]*(sommet[face[i].A][1] - sommet[face[i].C][1])
			+sommet[face[i].C][0]*(sommet[face[i].A][1] - sommet[face[i].B][1]);

		face[i].h = face[i].a*sommet[face[i].A][0]+
					face[i].b*sommet[face[i].A][1]+
					face[i].c*sommet[face[i].A][2];
		
		unit=sqrt(CARRE(face[i].a) + CARRE(face[i].b) + CARRE(face[i].c));
		
		/* Normalisation des coefficients du plan'appuyé sur la face */
		
		face[i].a /= unit;
		face[i].b /= unit;
		face[i].c /= unit;
		face[i].h /= unit;
	}
	
	/* Elimination et traçage du segment PQ */

	for(i=0; i<nbrligne; i++) {
		int l1= ligne[i][0],l2 = ligne[i][1];
		trace(0,sommet[l1][0],sommet[l1][1],sommet[l1][2],
			  sommet[l2][0],sommet[l2][1],sommet[l2][2]);
	}
}

/* --------------------------------------------- */
/* Traçage d'un segment après test de visibilité */
/* --------------------------------------------- */

void trace(jj,xP,yP,zP,xQ,yQ,zQ)
double xP,yP,zP,xQ,yQ,zQ;
int jj;
{
	int j = jj, traceln =1,A,B,C,i,Pb,Qb,dehors,Pdehors,Qdehors,eA,eB,eC,somme;
	int xe1,ye1,xe2,ye2;
	
	double Cpos,Ppos,Qpos,eps1;
	double a,b,c,h,hP,hQ,r1,r2,r3,xA,yA,zA,xB,yB,zB,xC,yC,zC,dA,dB,dC,MIN,MAX;
	double lambda,mu,xmin,ymin,zmin,xmax,ymax,zmax,C1,C2,C3,K1,K2,K3,denom1,denom2;
	
	/* Pour toutes les faces */
	while(j<nbrface) {
	
		a = face[j].a; /* Coefficients du plan associé à la face */
		b = face[j].b;
		c = face[j].c;
		h = face[j].h;
		eps1 = peps + (peps*h); /* Calcul de l'erreur relative */
	
		if (h <= 0) { /* si h = 0 la face passe par l'observateur E et ne cache rien */
			j++; /* si h < 0 la face est une face arrière et ne cache rien */
			continue;
		}
	
		/* test 1*/
	
		hP = a*xP + b*yP + c*zP;
		hQ = a*xQ + b*yQ + c*zQ;
	
		if((hP <= (h+eps1)) && (hQ <= (h+eps1))) {
			j++;
			continue;
		} /* PQ n'est pas derrière ABC */
	
		/* test 2*/
	
		K1 = yP*zQ - yQ*zP;
		K2 = zP*xQ - zQ*xP;
		K3 = xP*yQ - xQ*yP;
		
		A = face[j].A;
		B = face[j].B;
		C = face[j].C;
		
		xA = sommet[A][0]; yA = sommet[A][1]; zA = sommet[A][2];
		xB = sommet[B][0]; yB = sommet[B][1]; zB = sommet[B][2];
		xC = sommet[C][0]; yC = sommet[C][1]; zC = sommet[C][2];
		
		dA = K1*xA + K2*yA + K3*zA;
		dB = K1*xB + K2*yB + K3*zB;
		dC = K1*xC + K2*yC + K3*zC;
		
		/* si dA,dB,dC sont de même signe les points A,B,C sont du même coté que EPQ */
		
		eA = (dA > peps) ? 1 : (dA < meps) ? -1 : 0;
		eB = (dB > peps) ? 1 : (dB < meps) ? -1 : 0;
		eC = (dC > peps) ? 1 : (dC < meps) ? -1 : 0;
		
		somme = eA+eB+eC;
		
		if(abs(somme) >= 2) {
			j++;
			continue;
		}
		
		/* Si ce test est un succès alors la ligne infinie PQ est extérieure à la pyramide EABC */
		/* dans le cas contraire il y a intersection */
		
		/* test 3 */
		
		Pdehors = Qdehors = 0;
		MIN = 1;
		MAX = 0;
		
		for(i=0; i<3; i++) {
			C1 = yA * zB - yB * zA;
			C2 = zA * xB - zB * xA; /* Cl*x + C2*y + C3*z = 0 est le plan EAB */
			C3 = xA * yB - xB * yA;
			
			Cpos = C1*xC + C2*yC + C3*zC;
			Ppos = C1*xP + C2*yP + C3*zP;
			Qpos = C1*xQ + C2*yQ + C3*zQ;
			
			denom1 = Qpos-Ppos;
			
			if (Cpos>peps) {
				Pb = (Ppos < meps);
				Qb = (Qpos < meps);
				dehors = (Pb && (Qpos <= peps)) || (Qb && (Ppos <= peps));
			}
			else
				if (Cpos<meps) {
					Pb = (Ppos > peps);
					Qb = (Qpos > peps);
					dehors = (Pb && (Qpos >= meps)) || (Qb && (Ppos >= meps));
				}
				else dehors = 1;
		
			if (dehors) break;
		
			lambda = (abs(denom1) <= peps) ? 1.e7 : -Ppos/denom1;
			
			Pdehors |= Pb;
			Qdehors |= Qb;
			denom2 = dB - dA;
			mu = (abs(denom2) <= peps) ? 1.e7 : -dA/denom2;
			
			if ((mu >= meps) && (mu <= unpeps) && (lambda >= meps) && (lambda <= unpeps)) {
				if (lambda < MIN) MIN = lambda;
				if (lambda > MAX) MAX = lambda;
			}
			
			echange(&xA, &xB, &xC);
			echange(&yA, &yB, &yC);
			echange(&zA, &zB, &zC);
			echange(&dA, &dB, &dC);
		}
		if (dehors) {
			j++;
			continue;
		}
		
		/* test4 */
		
		if (!(Pdehors || Qdehors)) {
			traceln = 0; /* PQ est invisible */
			break;
		}
		
		/* test5 */
		
		r1 = xQ - xP;
		r2 = yQ - yP;
		r3 = zQ - zP;
		
		xmin = xP + MIN*r1;
		ymin = yP + MIN*r2;
		zmin = zP + MIN*r3;
		
		if ((a*xmin + b*ymin + c*zmin) < (h-eps1)) {
			j++;
			continue;
		}
		
		xmax = xP + MAX*r1;
		ymax = yP + MAX*r2;
		zmax = zP + MAX*r3;
		
		if ((a*xmax + b*ymax + c*zmax) < (h-eps1)) {
			j++;
			continue;
		}

		/* si ce test est un succès il existe un point d'intersection entre PQ et ABC */
		
		/* test6 */
		
		if (Pdehors || (hP < (h-eps1))) 
			trace(j+1,xP, yP, zP, xmin, ymin, zmin);
		
		if (Qdehors || (hQ < (h-eps1)))
			trace(j+1,xQ, yQ, zQ, xmax, ymax, zmax);
		
		traceln = 0;
		break;
	}
	
	if (traceln) {
		xe1 = (int)(320 + (xP/zP)*100);
		ye1 = (int)(256 + (yP/zP)*100); /* Calcul des coordonnées écran */
		xe2 = (int)(320 + (xQ/zQ)*100);
		ye2 = (int)(256 + (yQ/zQ)*100);
		
		/* Attention il n'y a pas de test de débordement */
		
		Move(&rastport[dbuff[!bascule]],xe1,ye1);
		Draw(&rastport[dbuff[!bascule]],xe2,ye2);
	}
}

void echange(a,b,c)
double *a,*b,*c;
{
	double aux = *a;
	*a = *b;
	*b = *c;
	*c = aux;
}