/*****************************************************************************
*   Routines to	test all edges for visibility relative to the global polygon *
* list,	and output the visible ones, to	graphics screen	or as text file	     *
*									     *
* Written by:  Gershon Elber				Ver 2.0, Jan. 1989   *
*****************************************************************************/

#ifdef __MSDOS__
#include <conio.h>
#endif /* __MSDOS__ */

#include <stdio.h>
#include <math.h>
#include <setjmp.h>
#include <time.h>
#include "program.h"
#include "genmat.h"
#include "parser.h"

static struct EdgeStruct *OutEdgeHashTable[EDGE_HASH_TABLE_SIZE];

static void OutputEdge(EdgeStruct *PEtemp);
static int EmitOutEdgeHashTable(void);
static struct EdgeStruct *FoundMatchEdge(int *Colinear, EdgeStruct *PEdge);
static int FuzzSamePoints(float Pt1[3], float Pt2[3]);
static int ColinearPoints(float Pt1[3], float Pt2[3], float Pt3[3]);
static int VisibleEdge(int EdgeMinY, EdgeStruct *PEdge,
					PolygonStruct *PolyHashTable[]);
static int VisibleEdgeOnePoly(EdgeStruct *PEdge, PolygonStruct *PPoly);
static int InverseMatrix(MatrixType A);
static void SwapDouble(double *D1, double *D2);

/*****************************************************************************
* Routine to test for visibility all edges in EdgeHashTable and	display	or   *
* output the visible ones only.	It is assumed that only	totally	visible	or   *
* invisible edges are in table (Pass 3 broke all other kind of edges).	     *
*****************************************************************************/
void OutVisibleEdges(void)
{
    int	i, PolyLineCount = 0;
    long SaveTime = time(NULL);
    struct EdgeStruct *PEtemp, *PEtail;

    fprintf(stderr, "\nPass 4, Edges [%5d] =      ", EdgeCount);
    for (i=0; i<EDGE_HASH_TABLE_SIZE; i++) OutEdgeHashTable[i] =
							(EdgeStruct *) NULL;
    if (!InverseMatrix(GlblViewMat)) {
	fprintf(stderr,
	        "\nNo inverse for matrix transformation, output is in screen space\n");
	GenUnitMat(GlblViewMat);	     /* No inverse transform at all. */
    }
    printf("Automatic Hidden Line generator\t\tGershon Elber\t\tJan. 1989\n");

    EdgeCount =	0;

    for (i=0; i<EDGE_HASH_TABLE_SIZE; i++)
	if ((PEtemp = EdgeHashTable[i]) != NULL)/* If any edge in that list. */
	    while (PEtemp) {
		EdgeCount++;
		fprintf(stderr, "\b\b\b\b\b%5d", EdgeCount);
		PEtail = PEtemp	-> Pnext;	   /* OutputEdge destroy it. */
		if ((!PEtemp -> Internal || InternalFlag) &&
		    VisibleEdge(i, PEtemp, PolyHashTable))
		    OutputEdge(PEtemp);
		PEtemp = PEtail;
	    }

    PolyLineCount =	EmitOutEdgeHashTable();     /* Output all polylines. */
    printf("[OBJECT Hidden");	       /* Output the OBJECT keyword: */
    for (i=0; i<PolyLineCount; i++)	{
	if (i % 10 == 0) printf("\n\t");
	printf("P%-3d ", i);
    }
    printf("]\n");
    fprintf(stderr, ",  %d seconds.", time(NULL) - SaveTime);
}

/*****************************************************************************
* Routine to output one	edge to	the OutEdgeHashTable.			     *
* The edge is inserted into OutEdgeHashTable to	be able	to match it into     *
* the longest path possible - connecting edges into edge sequence.	     *
* Each edge is inserted	by its Ymin, which will	cause the paths	be in	     *
* increased Y order...							     *
*****************************************************************************/
static void OutputEdge(EdgeStruct *PEtemp)
{
    int	Level;

    Level = (int) ((PEtemp -> Vertex[0]	-> Coord[1] + 1.0) *
							EDGE_HASH_TABLE_SIZE2);
    PEtemp -> Pnext = OutEdgeHashTable[Level];
    OutEdgeHashTable[Level] = PEtemp;
}

/*****************************************************************************
* Routine to scan the OutEdgeHashTable,	and Emit the longest paths found in  *
* it by	searching for consecutive edges	and combining colinear edges.	     *
*****************************************************************************/
static int EmitOutEdgeHashTable(void)
{
    int	i, j, Colinear, VertexCount = 0, VertexBaseCount, PolylineCount = 0;
    double Vec[3];
    struct EdgeStruct *PEtemp, *PEnext;

    for	(i=0; i<EDGE_HASH_TABLE_SIZE; i++)	 /* Scan all the hash table. */
	while (OutEdgeHashTable[i]) {		 /* Scan all edges in entry. */
	    PEtemp = OutEdgeHashTable[i];	  /* Take edge out of table. */
	    OutEdgeHashTable[i]	= OutEdgeHashTable[i] -> Pnext;
	    VertexBaseCount = VertexCount;

	    /* Print first vertex of polyline: */
	    MultVecby4by4(Vec, PEtemp -> Vertex[0] -> Coord, GlblViewMat);
	    printf("[VERTEX V%-4d %11lf %11lf %11lf]\n", VertexCount++,
		Vec[0], Vec[1], Vec[2]);

	    while ((PEnext = FoundMatchEdge(&Colinear, PEtemp)) != NULL) {
		if (!Colinear) {       /* Print PEtemp only if not colinear! */
		    /* Print next vertices of polyline:	*/
		    for	(j=0; j<3; j++)	Vec[j] =
			PEtemp -> Vertex[1] -> Coord[j];
		    MultVecby4by4(Vec, PEtemp -> Vertex[1] -> Coord,
								GlblViewMat);
		    printf("[VERTEX V%-4d %11lf %11lf %11lf]\n", VertexCount++,
			Vec[0], Vec[1], Vec[2]);
		}
		PEtemp = PEnext;
	    }

	    /* Print last vertex of polyline: */
	    MultVecby4by4(Vec, PEtemp -> Vertex[1] -> Coord, GlblViewMat);
	    printf("[VERTEX V%-4d %11lf %11lf %11lf]\n", VertexCount++,
		Vec[0], Vec[1], Vec[2]);

	    printf("[POLYLINE P%-4d", PolylineCount++);
	    for	(j=VertexBaseCount; j<VertexCount; j++)	{
		if ((j-VertexBaseCount)	% 10 ==	0) printf("\n\t");
		printf("V%-4d ", j);
	    }
	    printf("]\n");
	}
    return PolylineCount;
}

/*****************************************************************************
* Routine to scan the OutEdgeHashTable for a match edge	if any,	delete it    *
* from HashTable and return it.	if colinear with PEdge set Colinear to TRUE. *
* return FALSE if no match found (end of polyline).			     *
*****************************************************************************/
static struct EdgeStruct *FoundMatchEdge(int *Colinear, EdgeStruct *PEdge)
{
    int	Level;
    struct EdgeStruct *PEtemp, *PElast;

    Level = (int) ((PEdge -> Vertex[1] -> Coord[1] + 1.0) *
							EDGE_HASH_TABLE_SIZE2);
    PEtemp = PElast = OutEdgeHashTable[Level];
    while (PEtemp) {			    /* First scan for colinear edge. */
	if (FuzzSamePoints(PEdge -> Vertex[1] -> Coord,
			   PEtemp -> Vertex[0] -> Coord) &&
	    ColinearPoints(PEdge -> Vertex[0] -> Coord,
			   PEdge -> Vertex[1] -> Coord,
			   PEtemp -> Vertex[1] -> Coord)) {
	    *Colinear =	TRUE;
	    /* Delete that edge	from hash table	structure: */
	    if (PEtemp == PElast)
		 OutEdgeHashTable[Level] = OutEdgeHashTable[Level] -> Pnext;
	    else PElast	-> Pnext = PEtemp -> Pnext;

	    if (FuzzSamePoints(PEtemp -> Vertex[0] -> Coord,
			       PEtemp -> Vertex[1] -> Coord))
		 return	FoundMatchEdge(Colinear, PEdge);/* Call recursively! */
	    else return	PEtemp;
	}
	PElast = PEtemp;
	PEtemp = PEtemp	-> Pnext;
    }

    PEtemp = PElast = OutEdgeHashTable[Level];
    while (PEtemp) {			    /* Next scan for any match edge. */
	if (FuzzSamePoints(PEdge -> Vertex[1] -> Coord,
			   PEtemp -> Vertex[0] -> Coord)) {
	    *Colinear =	FALSE;
	    /* Delete that edge	from hash table	structure: */
	    if (PEtemp == PElast)
		 OutEdgeHashTable[Level] = OutEdgeHashTable[Level] -> Pnext;
	    else PElast	-> Pnext = PEtemp -> Pnext;

	    if (FuzzSamePoints(PEtemp -> Vertex[0] -> Coord,
			       PEtemp -> Vertex[1] -> Coord))
		 return	FoundMatchEdge(Colinear, PEdge);/* Call	recursively! */
	    else return	PEtemp;
	}
	PElast = PEtemp;
	PEtemp = PEtemp	-> Pnext;
    }

    return (EdgeStruct *) NULL;		      /* No match - end of polyline. */
}

/*****************************************************************************
* Routine to test if two points	are almost the same, and returns TRUE if so. *
*****************************************************************************/
static int FuzzSamePoints(float Pt1[3], float Pt2[3])
{
    return (APX_EQ(Pt1[0], Pt2[0]) &&
	    APX_EQ(Pt1[1], Pt2[1]) &&
	    APX_EQ(Pt1[2], Pt2[2]));
}

/*****************************************************************************
* Routine to test if three points are colinear,	and returns TRUE if so...    *
* Algorithm: cross product should be zero if colinear...		     *
* Note this routine does not return TRUE if Pt2	is not between Pt1..Pt3.     *
*****************************************************************************/
static int ColinearPoints(float Pt1[3], float Pt2[3], float Pt3[3])
{
    int	i;
    float Xout, Yout, Zout, U[3], V[3], temp;

    for	(i=0; i<3; i++)	{
	U[i] = Pt2[i] -	Pt1[i];
	V[i] = Pt3[i] -	Pt2[i];
    }
    temp = sqrt(SQR(U[0]) + SQR(U[1]) + SQR(U[2]));
    if (APX_EQ(temp, 0.0)) return TRUE;
    for	(i=0; i<3; i++)	U[i] = U[i] / temp;

    temp = sqrt(SQR(V[0]) + SQR(V[1]) + SQR(V[2]));
    if (APX_EQ(temp, 0.0)) return TRUE;
    for	(i=0; i<3; i++)	V[i] = V[i] / temp;

    /* Xoutput = Uy * Vz - Uz *	Vy. */
    Xout = U[1]	/* Uy */  * V[2] /* Vz */  -
	   U[2]	/* Uz */  * V[1] /* Vy */;

    /* Youtput = Uz * Vx - Ux *	Vz. */
    Yout = U[2]	/* Uz */  * V[0] /* Vx */  -
	   U[0]	/* Ux */  * V[2] /* Vz */;

    /* Zoutput = Ux * Vy - Uy *	Vx. */
    Zout = U[0]	/* Ux */  * V[1] /* Vy */  -
	   U[1]	/* Uy */  * V[0] /* Vx */;

    if (APX_EQ(Xout, 0.0) && APX_EQ(Yout, 0.0) && APX_EQ(Zout, 0.0) &&
	((MIN(Pt1[0], Pt3[0]) < Pt2[0]) && (MAX(Pt1[0], Pt3[0]) > Pt2[0]) &&
	 (MIN(Pt1[1], Pt3[1]) < Pt2[1]) && (MAX(Pt1[1], Pt3[1]) > Pt2[1])))
	 return	TRUE;
    else return	FALSE;
}

/*****************************************************************************
* Routine to test the visibility of the	given edge relative to all polygons  *
* in polygon list. return TRUE if the edge is visible. It is assumed that    *
* the edge is whole visible or whole invisible (Pass 3 broke the edge if     *
* that whas not	true). Also it is assumed the polygons are all convex.	     *
* A short cut is made to test the edge only against the	needed polygons	     *
* only, using the polygon hash table and Y level sorting.		     *
*****************************************************************************/
static int VisibleEdge(int EdgeMinY, EdgeStruct * PEdge,
						PolygonStruct *PolyHashTable[])
{
    int	EdgeMaxY, i;
    PolygonStruct *PPolyList, *PPolyLast;

    EdgeMaxY = MIN((MAX(PEdge -> Vertex[0] -> Coord[1],
			PEdge -> Vertex[1] -> Coord[1]) + 1.0) *
						EDGE_HASH_TABLE_SIZE2,
		   EDGE_HASH_TABLE_SIZE1);

    /* Test the	edge only in the interval 0..EdgeMaxY as these are the only  */
    /* which might hide	that edge. Also	polygons with MaxY less	than	     */
    /* EdgeMinY	are deleted from PolyHashTable as it is	assumed	the edges    */
    /* are scanned in increasing order of the EdgeHashTable.		     */
    for	(i=0; i<=EdgeMaxY; i++) {
	PPolyList = PPolyLast =	PolyHashTable[i];
	while (PPolyList) {
	    if ((PPolyList -> Ymax + 1.0) * EDGE_HASH_TABLE_SIZE2 < EdgeMinY)
		/* Delete this polygon!	*/
		if (PPolyList == PPolyLast) {
		    PolyHashTable[i] = PPolyLast = PPolyList -> Pnext;
		    PPolyList =	PPolyList -> Pnext;
		}
		else {
		    PPolyList =	PPolyList -> Pnext;
		    PPolyLast -> Pnext = PPolyList;
		}
	    else {	     /* Polygon still active - test edge against it. */
		/* If found one	polygon	that hides this	edge return FALSE... */

		if (!VisibleEdgeOnePoly(PEdge, PPolyList)) return FALSE;
		PPolyLast = PPolyList;
		PPolyList = PPolyList -> Pnext;
	    }
	}
    }

    return TRUE;
}

/*****************************************************************************
* Routine to test the visibility of the	given edge relative to one polygon   *
* return TRUE if the edge is visible. It is assumed that the edge is whole   *
* visible or whole invisible (Pass 3 broke the edge if that whas not true).  *
* Also it is assumed the polygons are all convex.			     *
*****************************************************************************/
static int VisibleEdgeOnePoly(EdgeStruct *PEdge, PolygonStruct *PPoly)
{
    int	i;
    float MidPt[3], Temp;
    struct VertexStruct *PList;

    for	(i=0; i<3; i++)	MidPt[i] =	    /* Calc a mid point on the edge: */
		(PEdge -> Vertex[0] -> Coord[i] +
		 PEdge -> Vertex[1] -> Coord[i]) / 2.0;

    /* Test if edges is	out of polygon boundary	box: */
    if (MidPt[0] > PPoly -> Xmax || MidPt[0] < PPoly ->	Xmin ||
	MidPt[1] > PPoly -> Ymax || MidPt[1] < PPoly ->	Ymin) return TRUE;

    Temp = PPoly -> Plane[0] * MidPt[0]	+
	   PPoly -> Plane[1] * MidPt[1]	+
	   PPoly -> Plane[2] * MidPt[2]	+
	   PPoly -> Plane[3];
    if ((Temp >= 0) || (APX_EQ(Temp, 0.0)))
	return TRUE;				/* Edge in front of polygon. */

    /* We cannt	escape it - we now must	find if	the point is included in */
    /* polygon area. We	assume the polygon is convex and use the fact	 */
    /* that the	polygon	list is	ordered	such that the cross product	 */
    /* of two consecutive vertices and the point is positive for all	 */
    /* consecutive vertices pairs iff the point	is in the polygon.	 */
    PList = PPoly -> PVertex;
    while (PList -> Pnext) {
	if (CrossProd(PList -> Coord, PList -> Pnext -> Coord, MidPt) < 0)
	    return TRUE;				  /* Out of polygon! */
	PList =	PList -> Pnext;
    }
    /* Now test	last polygon edge (last	point, first point in poly list) */
    if (CrossProd(PList -> Coord, PPoly -> PVertex -> Coord, MidPt) < 0)
	return TRUE;					  /* Out of polygon! */

    return FALSE;
}

/*****************************************************************************
* Procedure to calculate the INVERSE of	a given	MatrixType matrix A which is *
* modified to the inverted matrix. Return TRUE if inverse exists.	     *
*****************************************************************************/
static int InverseMatrix(MatrixType A)
{
    MatrixType InvA;
    int	i, j, k;
    double V;

    GenUnitMat(InvA);

    for	(i=0; i<4; i++)	{
	V = A[i][i];				      /* Find the new pivot. */
	k = i;
	for (j=i+1; j<4; j++) if (A[j][i] > V) {
	    /* Find maximum on col i, row i+1..4. */
	    V =	A[j][i];
	    k =	j;
	}
	j = k;

	if (i != j) for	(k=0; k<4; k++)	{
	    SwapDouble(&A[i][k], &A[j][k]);
	    SwapDouble(&InvA[i][k], &InvA[j][k]);
	}

    	for (j=i+1; j<4; j++) {		 /* Eliminate col i from row i+1..4. */
	    V =	A[j][i]	/ A[i][i];
	    for	(k=0; k<4; k++)	{
		A[j][k]	   -= V	* A[i][k];
		InvA[j][k] -= V	* InvA[i][k];
	    }
	}
    }

    for	(i=4-1;	i>=0; i--) {			       /* Back Substitution. */
	if (A[i][i] == 0) return FALSE;	   /* Error */

	for (j=0; j<i; j++) {		 /* Eliminate col i from row 1..i-1. */
	    V =	A[j][i]	/ A[i][i];
	    for	(k=0; k<4; k++)	{
		/* A ->	[j][k] -= V * A[i][k]; */
		InvA[j][k] -= V	* InvA[i][k];
	    }
	}
    }

    for	(i=0; i<4; i++)			    /* Normalize the inverse Matrix. */
	for (j=0; j<4; j++)
	    InvA[i][j] /= A[i][i];
    for	(i=0; i<4; i++)			       /* Copy inverted	matrix to A. */
	for (j=0; j<4; j++)
	    A[i][j] = InvA[i][j];

    return TRUE;
}

/*****************************************************************************
* Procedure to swap two	double parameters:				     *
*****************************************************************************/
static void SwapDouble(double *D1, double *D2)
{
    double D;

    D =	(*D1);
    *D1	= (*D2);
    *D2	= D;
}