/*****************************************************************************
*   Routines to	evaluate vertices colors for all polygons in data. It is     *
* assumed the parser has updated plane equations for all polygons at this    *
* stage (if from reading it in, or by regenerating it), and they were mapped *
* correctly using the global transformation applied to them.		     *
*   The vertices are updated for both Flat or Gouraud shading, and color is  *
* evaluated for them, as a color index into the color table. Because of the  *
* way the table is constructed and because two vertices in same polygon can  *
* have two levels of same color, it is o.k. to interplated the color	     *
* indices (We assume no specular affects)!				     *
*   As side affect, we use this pass to hash the polygon into global hash    *
* table PolyHasTable, sorted by their Ymin values.			     *
*									     *
* Written by:  Gershon Elber				Ver 2.0, Mar. 1990   *
*****************************************************************************/

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

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

#define	MIN_VERTEX_HASH_SIZE	1000  /* Minimum number of entries in table. */

typedef struct VertexHashStruct {	 /* Used to hash an object vertices. */
    int Entry;
    float Normal[3];
    VertexStruct *PVertex;
} VertexHashStruct;

static float MinNormalDotProd; /* Min. normal dot product to accept as same. */
static int VertexHashSize = 0;	/* Size of vertex hash table size allocated. */
static int PolyCount = 0;    /* Number of polygons processed in this module: */

static void PrepareAllObjects(FileDescription **FD);
static void VisitObjectTree(BinTree *PBinTree);
static void PrepareOneObject(ObjectStruct *PObject);
static void AverageVrtxNormals(ObjectStruct *PObject);
static void InsertVerticesToHashTable(VertexHashStruct *VerticesHashTable,
							ObjectStruct *PObject);
static void UpdateVerticesFromHashTable(VertexHashStruct *VerticesHashTable,
							ObjectStruct *PObject);
static int SameVertexAndNormals(VertexStruct *PVertex1, float *Normal1,
						VertexHashStruct *PHVertex2);
static float EvalNormalsAngle(VertexStruct *PVertex1, float *Normal1,
						VertexHashStruct *PHVertex2);
static int EvalVertexKey(VertexStruct *PVertex);
static void EvalOnePolygon(PolygonStruct *PPolygon, int ColorIndex);
static int CosineShading(float *Normal, int ColorIndex);

/*****************************************************************************
* Routine to evaluate the vertices colors for both Flat/Gouraud shading, as  *
* indices for the color table.						     *
*****************************************************************************/
void EvalVrtxColors(FileDescription **FD, int NumOfObjects, char **Objects)
{
    int	i;
    long SaveTime = time(NULL);
    struct ObjectStruct	*PObject;

    PolyHashTable = (PolygonStruct **)
	MyMalloc(sizeof(PolygonStruct *) * ScreenYSize *
					ShadingInfo.SubSamplePixel);

    for	(i=0; i<ScreenYSize * ShadingInfo.SubSamplePixel; i++)
	PolyHashTable[i] = NULL;

    fprintf(stderr, "\nPass 3, Polys [%4d] =      ", NumOfPolygons);

    if (NumOfObjects > 0)	   /* There was something on command line... */
	for (i=0; i<NumOfObjects; i++) {
	    if ((PObject=SearchObject(FD, *Objects)) != (ObjectStruct *) NULL)
		PrepareOneObject(PObject);
	}
    else {		    /* Prepare all objects by scanning object trees. */
	PrepareAllObjects(FD);
    }

    fprintf(stderr, ",  %ld seconds.", time(NULL) - SaveTime);
}

/*****************************************************************************
* Scan all objects.							     *
*****************************************************************************/
static void PrepareAllObjects(FileDescription **FD)
{
    while (*FD)	VisitObjectTree((*FD++) -> ObjectPointer);
}

/*****************************************************************************
* Scanning all the object in tree PBinTree and prepare them.		     *
*****************************************************************************/
static void VisitObjectTree(BinTree *PBinTree)
{
    if (PBinTree == (BinTree *)	NULL) return;

    VisitObjectTree(PBinTree -> right);

    PrepareOneObject(PBinTree -> Data.PObject);

    VisitObjectTree(PBinTree -> left);
}

/*****************************************************************************
* Routine to prepare one object Object.					     *
*****************************************************************************/
static void PrepareOneObject(ObjectStruct *PObject)
{
    static int WasPolyline = FALSE;
    int Level, ColorIndex = PObject -> Color;
    struct PolygonStruct *Ptemp, *PList = PObject -> PPolygon;

    if (GouraudFlag) {
	/* Need to average all vertices normals from adjacent polygons,      */
	/* so we can interpolate them, and select color for them directly.   */
	AverageVrtxNormals(PObject);
    }

    while (PList) {
	Ptemp = PList -> Pnext;

	EvalOnePolygon(PList, ColorIndex);

	/* And add polygon into polygon hash table sorted by Ymin: */
	if (!PList -> Polyline) {
	    /* If BackFacingFlag is TRUE, and this polygon has plane     */
	    /* suggesting it is back facing - dont insert to hash table. */
	    if (!(BackFacingFlag && PList -> Plane[2] < 0.0)) {
		Level = PList -> Ymin;
		Level = BOUND(Level, 0, ScreenYSize *
				ShadingInfo.SubSamplePixel);/* Be 100% safe. */
		if (Level < ScreenYSize * ShadingInfo.SubSamplePixel) {
		    /* Concat to hash table at level Level: */
		    PList -> Pnext = PolyHashTable[Level];
		    PolyHashTable[Level] = PList;
		}
	    }
	    else {
		if (!GouraudFlag) PolyCount--;	      /* One less polygon... */
	    }
	}
	else {
	    if (!WasPolyline) {
		WasPolyline = TRUE;
		if (MoreFlag)
		    fprintf(stderr, "\nPolylines are not supported, ignored\n");
	    }
	}

	PList =	Ptemp;
    }
}

/*****************************************************************************
* Routine to update average normals of adjacent polygons at common vertices  *
* so gouraud shading can be applied. Does 2 passes on all object vertices:   *
* 1. Combine (and average) common vertices - vertices with normal no differ  *
*    more than NormalAvgDegree degrees.					     *
* 2. For each vertex find its closest averaged normal as evaluated at 1, and *
*    calculate color for it.						     *
*****************************************************************************/
static void AverageVrtxNormals(ObjectStruct *PObject)
{
    int i;
    VertexHashStruct *VerticesHashTable;

    /* Prepare maximum degree allowed to merge two normals in cosine form:   */
    MinNormalDotProd = cos(NormalAvgDegree * M_PI / 180.0);

    /* Allocate hash table twice as big as number of possible entries to     */
    /* reduce the average hit ratio, with minimum of MIN_VERTEX_HASH_SIZE.   */
    VertexHashSize = MAX(NumOfVertices * 2, MIN_VERTEX_HASH_SIZE);
    VerticesHashTable = (VertexHashStruct *)
	MyMalloc(VertexHashSize * sizeof(VertexHashStruct));
    for (i=0; i<VertexHashSize; i++) VerticesHashTable[i].Entry = 0;

    InsertVerticesToHashTable(VerticesHashTable, PObject);

    UpdateVerticesFromHashTable(VerticesHashTable, PObject);

    free((char *) VerticesHashTable);
}

/*****************************************************************************
* Routine to insert all vertices in object into the hash table. Each vertex  *
* is entered at place (key) as select via EvalVertexKey routine. If no place *
* the next ones are scaned until free (NULL) is found (no double key fansy   *
* hashing techniques...). Note however that while scanning the non NULL      *
* entries the vertex position is compared for equality, and its normal for   *
* equality upto AvgNormalDegree, and if both hold, the two are merged into   *
* one vertex in that position but with averaged normal.			     *
*****************************************************************************/
static void InsertVerticesToHashTable(VertexHashStruct *VerticesHashTable,
							ObjectStruct *PObject)
{
    int i, Key;
    float EntryRatio, *Normal, *OldNormal;
    PolygonStruct *PPolygon;
    VertexStruct *PVertex;

    for (PPolygon = PObject -> PPolygon; PPolygon != NULL;
	 PPolygon = PPolygon -> Pnext) {
	Normal = PPolygon -> Plane;

	/* If BackFacingFlag is TRUE, and this polygon has plane    */
	/* suggesting it is back facing - dont count it.	    */
	if (!(BackFacingFlag && PPolygon -> Plane[2] < 0.0))
	    PolyCount++;
	fprintf(stderr, "\b\b\b\b\b%5d", PolyCount);

	for (PVertex = PPolygon -> PVertex; PVertex != NULL;
	     PVertex = PVertex -> Pnext) {
	    Key = EvalVertexKey(PVertex);
	    while (VerticesHashTable[Key].Entry != 0) {
		/* Test if should be combined with old vertex: */
		if (SameVertexAndNormals(PVertex, Normal,
						    &VerticesHashTable[Key])) {
		    /* Megre the normals: */
		    EntryRatio = 1.0 / ++VerticesHashTable[Key].Entry;
		    OldNormal = VerticesHashTable[Key].Normal;
		    for (i=0; i<3; i++) OldNormal[i] =
			OldNormal[i] * (1.0 - EntryRatio) +
			Normal[i] * EntryRatio;
		    break;
		}
		Key = ++Key % VertexHashSize;
	    }
	    if (VerticesHashTable[Key].Entry == 0) {
		/* Could not merge the vertex with old one - do it now: */
		VerticesHashTable[Key].PVertex = PVertex;
		GEN_COPY(VerticesHashTable[Key].Normal, Normal,
							3 * sizeof(float));
		++VerticesHashTable[Key].Entry;
	    }
	}
    }
}

/*****************************************************************************
* Routine to scan all vertices in object and update their normals to the     *
* close one at that point as found in the hash table.			     *
*****************************************************************************/
static void UpdateVerticesFromHashTable(VertexHashStruct *VerticesHashTable,
							ObjectStruct *PObject)
{
    static Flip = FALSE;
    int Key;
    float *Normal, MaxCosAngle, CosAngle;
    PolygonStruct *PPolygon;
    VertexStruct *PVertex;
    VertexHashStruct *PHVertex;

    for (PPolygon = PObject -> PPolygon; PPolygon != NULL;
	 PPolygon = PPolygon -> Pnext) {
	Normal = PPolygon -> Plane;

	if (Flip) {
	    Flip = FALSE;
	    fprintf(stderr, "a\b");
	}
	else {
	    Flip = TRUE;
	    fprintf(stderr, "A\b");
	}

	for (PVertex = PPolygon -> PVertex; PVertex != NULL;
	     PVertex = PVertex -> Pnext) {
	    PHVertex = NULL;
	    MaxCosAngle = -INFINITY;
	    Key = EvalVertexKey(PVertex);
	    while (VerticesHashTable[Key].Entry != 0) {
		/* Search for the closest Normal at this vertex point: */
		if ((CosAngle = (EvalNormalsAngle(PVertex, Normal,
				    &VerticesHashTable[Key]))) > MaxCosAngle) {
		    PHVertex = &VerticesHashTable[Key];
		    MaxCosAngle = CosAngle;
		}
		Key = ++Key % VertexHashSize;
	    }
	    if (PVertex -> HasNormal) {
		/* Use defined normal for this vertex, if has one. */
		PVertex -> Color = CosineShading(PVertex -> Normal,
							PObject -> Color);
	    }
	    else
	    if (PHVertex) {    /* Should always be non NULL but who knows... */
		PVertex -> Color = CosineShading(PHVertex -> Normal,
							PObject -> Color);
	    }
	    else {
		PVertex -> Color = CosineShading(PPolygon -> Plane,
							PObject -> Color);
	    }
	}
    }
}

/*****************************************************************************
* Routine to compare two vertices to be exactly equal in their position and  *
* to have same normal up to NormalAvgDegree.				     *
*****************************************************************************/
static int SameVertexAndNormals(VertexStruct *PVertex1, float *Normal1,
						VertexHashStruct *PHVertex2)
{
    float *Coord1 = PVertex1 -> Coord,
	  *Coord2 = PHVertex2 -> PVertex -> Coord,
	  *Normal2;

    if (!APX_EQ(Coord1[0], Coord2[0]) ||
	!APX_EQ(Coord1[1], Coord2[1]) ||
	!APX_EQ(Coord1[2], Coord2[2])) return FALSE;

    Normal2 = PHVertex2 -> Normal;

    return (Normal1[0] * Normal2[0] +
	    Normal1[1] * Normal2[1] +
	    Normal1[2] * Normal2[2] > MinNormalDotProd);
}

/*****************************************************************************
* Routine to evaluate the angle between the given two vertices if they are   *
* the same, and return the angle cosine value. (-INFINITY is returned if     *
* not same point...).							     *
*****************************************************************************/
static float EvalNormalsAngle(VertexStruct *PVertex1, float *Normal1,
						VertexHashStruct *PHVertex2)
{
    float *Coord1 = PVertex1 -> Coord,
	  *Coord2 = PHVertex2 -> PVertex -> Coord,
	  *Normal2;

    if (!APX_EQ(Coord1[0], Coord2[0]) ||
	!APX_EQ(Coord1[1], Coord2[1]) ||
	!APX_EQ(Coord1[2], Coord2[2])) return -INFINITY;

    Normal2 = PHVertex2 -> Normal;

    return Normal1[0] * Normal2[0] +
	   Normal1[1] * Normal2[1] +
	   Normal1[2] * Normal2[2];
}

/*****************************************************************************
* Routine to evaluate integer key in the range 0 .. VertexHashSize - 1	     *
* for the given vertex PVertex.						     *
*****************************************************************************/
static int EvalVertexKey(VertexStruct *PVertex)
{
    int Key;
    float *Coord = PVertex->Coord;

    Key = ((int) (((long) (Coord[0] * 239 + Coord[1] * 677 + Coord[2] * 109)) %
							VertexHashSize));

    return Key;
}

/*****************************************************************************
* Routine to update the vertices colors, which are indices tp the ColorMap   *
* as save in the global ShadingInfo structure. Given Object color (as index  *
* into table), the exact level is evaluationed using each polygon normal.    *
*****************************************************************************/
static void EvalOnePolygon(PolygonStruct *PPolygon, int ColorIndex)
{
    int Color;
    struct VertexStruct *V = PPolygon -> PVertex;

    if (GouraudFlag) {
	/* Vertices were already updated by the Object averaging of normals  */
	/* routine (see AverageVrtxNormals routine).			     */
    }
    else {
	fprintf(stderr, "\b\b\b\b\b%5d", ++PolyCount);

	/* This is much easier - set all vertices color to same Color.       */
	Color = CosineShading(PPolygon -> Plane, ColorIndex);
	do {
	    V -> Color = Color;
	    V = V -> Pnext;
	}
	while (V != NULL);
    }
}

/*****************************************************************************
* Routine to evaluate cosine shading given a normal vector. Uses the global  *
* shading structure ShadingInfo to evaluate it. Returns a color from the     *
* color map defined for the current scene.				     *
* It is assumed the given normal is normalized to size 1.0, and that the     *
* intensity levels of the Object Color are saved in indices PObject -> Color *
* to PObject -> Color + ShadingInfo.LevelsPerColor - 1.	These colors are     *
* interpolated to begin from the ambient level, so we dont have to consider  *
* ambient light here - just the cosine factor.				     *
*****************************************************************************/
static int CosineShading(float *Normal, int ColorIndex)
{
    int NewColorIndex;
    double Intensity;

    if (ShadingInfo.TwoSources) {
	/* Take the absolute value of the dot product as intensity: */
	Intensity = Normal[0] * ShadingInfo.LightSource[0] +
		    Normal[1] * ShadingInfo.LightSource[1] +
		    Normal[2] * ShadingInfo.LightSource[2];
	Intensity = ABS(Intensity);
    }
    else {
	/* Dot product of two normals is in [-1..1] range. Make it [0..1]: */
	Intensity = (Normal[0] * ShadingInfo.LightSource[0] +
		     Normal[1] * ShadingInfo.LightSource[1] +
		     Normal[2] * ShadingInfo.LightSource[2] + 1.0) * 0.5;
    }

    NewColorIndex = ColorIndex +
		((int) ((ShadingInfo.LevelsPerColor - 1) * (1.0 - Intensity)));
    /* This should never happen, but if it does, the error is so fatal       */
    /* (generate different color tone instead...) we double check this one.  */
    return BOUND(NewColorIndex, ColorIndex,
				ColorIndex + ShadingInfo.LevelsPerColor - 1);
}

/*****************************************************************************
* Routine to update plane equation of the given	polygon:		     *
*   It is assumed that at list 3 points in polygon do exists, and pick the   *
* tuple that has biggest length for maximum accuracy.			     *
*   Note the polygon must be convex.					     *
*   if SetFlipDir, then it is assumed Plane equation is given and the	     *
* FlipPlaneDir entry in PPolygon is set to know if direction should be	     *
* flipped.								     *
*   Returns FALSE if failed to evaluate the PLANE equation.		     *
*****************************************************************************/
int UpdateEqnPolygon(PolygonStruct *PPolygon, int SetFlipDir)
{
    int	i;
    float MaxLen = 0.0, Len, V1[3], V2[3], *Coord, *CoordNext, *CoordNextNext,
	Plane[3], MaxPlane[3];
    struct VertexStruct *PList = PPolygon -> PVertex;

    do {	/* Search for 3 consequtive non-colinear point from polygon: */
	Coord = PList -> Coord;
	CoordNext = PList -> Pnext -> Coord;
	CoordNextNext = PList -> Pnext -> Pnext -> Coord;

	if (!APX_PT_EQ(Coord, CoordNext) &&
	    !APX_PT_EQ(CoordNext, CoordNextNext)) {

	    for (i=0; i<3; i++) {    /* Prepare two vectors on polygon plane */
		V1[i] = Coord[i] - CoordNext[i];
		V2[i] = CoordNext[i] - CoordNextNext[i];
	    }

	    /* Find plane normal by a cross product of two vectors on plane: */
	    Plane[0] = V1[1] * V2[2] - V1[2] * V2[1];
	    Plane[1] = V1[2] * V2[0] - V1[0] * V2[2];
	    Plane[2] = V1[0] * V2[1] - V1[1] * V2[0];

	    /* Find vector Len. - we are looking for the biggest: */
	    Len = sqrt(SQR(Plane[0]) + SQR(Plane[1]) + SQR(Plane[2]));
	    if (Len > MaxLen) {
		for (i=0; i<3; i++) MaxPlane[i] = Plane[i];
		MaxLen = Len;
	    }
	}

	PList = PList -> Pnext;				  /* Try next tuple. */
    } while (PList -> Pnext -> Pnext != NULL);

    if (ABS(MaxLen) < SQR(EPSILON)) { /* Fail to find 3 non-colinear points. */
	if (MoreFlag) {
	    fprintf(stderr,
	        "\nError: Invalid polygon (%d) found in file (zero edge length/colinear vertices):\n",
		PolyCount);
	    PrintPolyContent(PPolygon);
	}
	return FALSE;
    }

    if (SetFlipDir) {
	PPolygon -> FlipPlaneDir =
	     PPolygon -> Plane[0] * MaxPlane[0] +
	     PPolygon -> Plane[1] * MaxPlane[1] +
	     PPolygon -> Plane[2] * MaxPlane[2] < 0.0;
    }

    for (i=0; i<3; i++)
	PPolygon -> Plane[i] = (PPolygon -> FlipPlaneDir ? -1 : 1)
						* MaxPlane[i] / MaxLen;

    PPolygon ->	Plane[3] =
	(- Coord[0] * PPolygon -> Plane[0]
	 - Coord[1] * PPolygon -> Plane[1]
	 - Coord[2] * PPolygon -> Plane[2]);

    return TRUE;
}

/*****************************************************************************
* Routine to evaluate the cross	product	of 3 points projected to Z = 0 plane *
* and return the sign of the result (Only Z component).			     *
*****************************************************************************/
int CrossProd(float Pt1[3], float Pt2[3], float Pt3[3])
{
    float Zout;

    /* U = Pt2 - Pt1,  V = Pt3 - Pt2,		Zoutput	= Ux * Vy - Uy * Vx. */
    Zout = (Pt2[0] - Pt1[0]) /*	Ux */  * (Pt3[1] - Pt2[1]) /* Vy */  -
	   (Pt2[1] - Pt1[1]) /*	Uy */  * (Pt3[0] - Pt2[0]) /* Vx */;
    if (APX_EQ(Zout, 0.0)) return 0;
    if (Zout < 0.0)
	 return	-1;
    else return	1;
}

/*****************************************************************************
* Routine to print the content of a given edge:				     *
*****************************************************************************/
void PrintPolyContent(PolygonStruct *PPoly)
{
    struct VertexStruct *PList = PPoly -> PVertex;

    while (PList) {
	fprintf(stderr, "   %12f %12f %12f\n",
	    PList -> Coord[0],
	    PList -> Coord[1],
	    PList -> Coord[2]);
	PList =	PList -> Pnext;
    }
}