/*****************************************************************************
*   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 "iritprsr.h"

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

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

static RealType MinNormalDotProd; /* Min. normal dot prod 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(IPObjectStruct *PObject);
static void PrepareOneObject(IPObjectStruct *PObject);
static void AverageVrtxNormals(IPObjectStruct *PObject);
static void InsertVerticesToHashTable(VertexHashStruct *VerticesHashTable,
						      IPObjectStruct *PObject);
static void UpdateVerticesFromHashTable(VertexHashStruct *VerticesHashTable,
						      IPObjectStruct *PObject);
static int SameVertexAndNormals(IPVertexStruct *PVertex1, RealType *Normal1,
						  VertexHashStruct *PHVertex2);
static RealType EvalNormalsAngle(IPVertexStruct *PVertex1, RealType *Normal1,
						  VertexHashStruct *PHVertex2);
static int EvalVertexKey(IPVertexStruct *PVertex);
static void EvalOnePolygon(IPPolygonStruct *PPolygon, int ColorIndex);
static int CosineShading(RealType *Normal, int ColorIndex);

/*****************************************************************************
* Routine to evaluate the vertices colors for both Flat/Gouraud shading, as  *
* indices for the color table.						     *
*****************************************************************************/
void EvalVrtxColors(IPObjectStruct *PObjects)
{
    int	i;
    long
	SaveTime = time(NULL);

    PolyHashTable = (IPPolygonStruct **)
	MyMalloc(sizeof(IPPolygonStruct *) * GlblShadeInfo.ScrnYSize *
						GlblShadeInfo.SubSamplePixel);

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

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

    PrepareAllObjects(PObjects);

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

/*****************************************************************************
* Scan all objects.							     *
*****************************************************************************/
static void PrepareAllObjects(IPObjectStruct *PObjects)
{
    while (PObjects) {
	PrepareOneObject(PObjects);
	PObjects = PObjects -> Pnext;
    }
}

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

    if (GlblShadeInfo.Gouraud) {
	/* 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: */
	switch (PList -> Type) {
	    case IP_POLYGON:
		/* If BackFacing is TRUE, and this polygon has plane         */
		/* suggesting it is back facing - dont insert to hash table. */
		if (!(GlblShadeInfo.BackFacing && PList -> Plane[2] > 0.0)) {
		    Level = PList -> Ymin;
		    Level = BOUND(Level, 0, GlblShadeInfo.ScrnYSize *
			     GlblShadeInfo.SubSamplePixel); /* Be 100% safe. */
		    if (Level < GlblShadeInfo.ScrnYSize *
						GlblShadeInfo.SubSamplePixel) {
			/* Concat to hash table at level Level: */
			PList -> Pnext = PolyHashTable[Level];
			PolyHashTable[Level] = PList;
		    }
		}
		else {
		    if (!GlblShadeInfo.Gouraud) PolyCount--;/* One less poly.*/
		}
		break;
	    case IP_POLYLINE:
		if (!WasPolyline) {
		    WasPolyline = TRUE;
		    if (GlblMore)
			fprintf(stderr,
			    "\nPolylines are not supported, ignored.\n");
		}
		break;
	}

	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 GlblShadeInfo.NrmlAvgDegree degrees.			     *
* 2. For each vertex find its closest averaged normal as evaluated at 1, and *
*    calculate color for it.						     *
*****************************************************************************/
static void AverageVrtxNormals(IPObjectStruct *PObject)
{
    int i;
    VertexHashStruct *VerticesHashTable;

    /* Prepare maximum degree allowed to merge two normals in cosine form:   */
    MinNormalDotProd = cos(GlblShadeInfo.NrmlAvgDegree * 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(GlblNumOfVerts * 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,
						      IPObjectStruct *PObject)
{
    int i, Key;
    RealType EntryRatio, *Normal, *OldNormal;
    IPPolygonStruct *PPolygon;
    IPVertexStruct *PVertex;

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

	/* If BackFacing is TRUE, and this polygon has plane        */
	/* suggesting it is back facing - dont count it.	    */
	if (!(GlblShadeInfo.BackFacing && 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(RealType));
		++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,
						      IPObjectStruct *PObject)
{
    static int
	Flip = FALSE;
    int Key;
    RealType *Normal, MaxCosAngle, CosAngle;
    IPPolygonStruct *PPolygon;
    IPVertexStruct *PVertex;
    VertexHashStruct *PHVertex;

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

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

	PVertex = PPolygon -> PVertex;
	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 (IP_HAS_VRTX_NORMAL(PVertex)) {
		/* 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 GlblShadeInfo.NrmlAvgDegree.		     *
*****************************************************************************/
static int SameVertexAndNormals(IPVertexStruct *PVertex1, RealType *Normal1,
						VertexHashStruct *PHVertex2)
{
    RealType *Normal2,
	*Coord1 = PVertex1 -> Coord,
	*Coord2 = PHVertex2 -> PVertex -> Coord;

    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 RealType EvalNormalsAngle(IPVertexStruct *PVertex1, RealType *Normal1,
						VertexHashStruct *PHVertex2)
{
    RealType *Normal2,
	*Coord1 = PVertex1 -> Coord,
	*Coord2 = PHVertex2 -> PVertex -> Coord;

    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(IPVertexStruct *PVertex)
{
    int Key;
    RealType
	*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 global GlblShadeInfo structure. Given Object color (as index    *
* into table), the exact level is evaluationed using each polygon normal.    *
*****************************************************************************/
static void EvalOnePolygon(IPPolygonStruct *PPolygon, int ColorIndex)
{
    int Color;
    IPVertexStruct
	*V = PPolygon -> PVertex;

    if (GlblShadeInfo.Gouraud) {
	/* 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);
	for (; V != NULL; V = V -> Pnext) V -> Color = Color;
    }
}

/*****************************************************************************
* Routine to evaluate cosine shading given a normal vector. Uses the global  *
* shading structure GlblShadeInfo 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 + GlblShadeInfo.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(RealType *Normal, int ColorIndex)
{
    int NewColorIndex;
    RealType Intensity;

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

    NewColorIndex = ColorIndex +
	      ((int) ((GlblShadeInfo.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 + GlblShadeInfo.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.					     *
*   Returns FALSE if failed to evaluate the PLANE equation.		     *
*****************************************************************************/
int UpdateEqnPolygon(IPPolygonStruct *PPolygon, int SetFlipDir)
{
    int	i, FlipPlaneDir;
    RealType Len, V1[3], V2[3], *Coord, *CoordNext, *CoordNextNext, Plane[3],
	MaxPlane[3],
	MaxLen = 0.0;
    IPVertexStruct
	*VList = PPolygon -> PVertex;

    /* Search for 3 consequtive non-colinear point from polygon: */
    for (; VList -> Pnext -> Pnext != NULL; VList = VList -> Pnext) {
	Coord = VList -> Coord;
	CoordNext = VList -> Pnext -> Coord;
	CoordNextNext = VList -> Pnext -> Pnext -> Coord;

	if (!PT_EQ(Coord, CoordNext) &&
	    !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;
	    }
	}
    }

    if (ABS(MaxLen) < SQR(EPSILON)) { /* Fail to find 3 non-colinear points. */
	if (GlblMore) {
	    fprintf(stderr,
		"\nError: Invalid polygon (%d) found in file (zero edge length/colinear vertices):\n",
		PolyCount);
	    PrintPolyContent(PPolygon);
	}
	PPolygon -> Plane[0] = PPolygon -> Plane[1] = 0.0;  /* Pick Z = 1.0. */
	PPolygon -> Plane[2] = PPolygon -> Plane[2] = 1.0;
	return FALSE;
    }

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

    for (i = 0; i < 3; i++)
	PPolygon -> Plane[i] = (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(RealType Pt1[3], RealType Pt2[3], RealType Pt3[3])
{
    RealType 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(IPPolygonStruct *PPoly)
{
    IPVertexStruct
	*VList = PPoly -> PVertex;

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