/******************************************************************************
* GeomVals.c - Area, Volume, and counts on polygonal objects.		      *
*******************************************************************************
* Written by Gershon Elber, March 1990.					      *
******************************************************************************/

#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include <string.h>
#include "allocate.h"
#include "convex.h"
#include "geomat3d.h"
#include "geomvals.h"

static RealType PolygonXYArea(IPVertexStruct *VHead);
static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3);

/*****************************************************************************
*  Routine to evaluate the Area of the given geom. object, in object unit.   *
* Algorithm (for each polygon):					V3	     *
* 1. Set Polygon Area to be zero.			       /\	     *
*    Make a copy of the original polygon		     /	  \ 	     *
*    and transform it to a XY parallel plane.		   /	    \V2	     *
*    Find the minimum Y value of the polygon	       V4/	     |	     *
*    in the XY plane					 \	     |	     *
* 2. Let V(0) be the first vertex, V(n) the last one       \         |	     *
*    for i goes from 0 to n-1 add to Area the area of	     \_______|       *
*    below edge V(i), V(i+1):				     V0      V1      *
*    PolygonArea += (V(i+1)x - V(i)x) * (V(i+1)y' - V(i)y') / 2		     *
*    where V(i)y' is V(i)y - MinY, where MinY is the polygon minimum Y value *
* 3. Note that the result of step 2 is the area of the polygon itself.       *
*    However it might be negative, so take the absolute result of step 2 and *
*    add it to the global ObjectArea.					     *
* Note step 2 is performed by another aux. routine below: PolygonXYArea.     *
*****************************************************************************/
double PolyObjectArea(IPObjectStruct *PObj)
{
    RealType
	ObjectArea = 0.0;
    MatrixType RotMat;
    IPPolygonStruct *Pl;
    IPVertexStruct *V, *VHead;

    if (!IP_IS_POLY_OBJ(PObj))
	IritFatalError("Geometric property requested on non polygonal object");

    if (IP_IS_POLYLINE_OBJ(PObj)) {
	return 0.0;
    }

    Pl = PObj -> U.Pl;
    while (Pl != NULL) {
	/* Dont trans. original object. */
	V = VHead = CopyVertexList(Pl -> PVertex);
	/* Create the trans. matrix to transform the polygon to XY parallel  */
	GenRotateMatrix(RotMat, Pl -> Plane);			   /* plane. */
	do {
	    MatMultVecby4by4(V -> Coord, V -> Coord, RotMat);

	    V = V -> Pnext;
	}
	while (V != NULL && V != VHead);

	ObjectArea += PolygonXYArea(VHead);

	IPFreeVertexList(VHead);		  /* Free the vertices list. */

	Pl = Pl -> Pnext;
    }

    return ObjectArea;
}

/*****************************************************************************
*  Routine to evaluate the area of the given polygon projection on the XY    *
* plane. Note the polygon does not have to be on a XY parallel plane, as     *
* only its XY projection is considered (Z is ignored). Returns the area of   *
* its XY parallel projection.						     *
*  See GeomObjectArea above for algorithm:				     *
*****************************************************************************/
static RealType PolygonXYArea(IPVertexStruct *VHead)
{
    RealType PolygonArea = 0.0, MinY;
    IPVertexStruct *V = VHead, *Vnext;

    MinY = V -> Coord[1];
    V = V -> Pnext;
    while (V != VHead && V != NULL /* Should not happen! */) {
	if (MinY > V -> Coord[1])
	    MinY = V -> Coord[1];

	V = V -> Pnext;
    }

    Vnext = V -> Pnext;
    MinY *= 2.0;		  /* Instead of subtracting twice each time. */
    do {
	/* Evaluate area below edge V-Vnext relative to Y level MinY. Note   */
	/* it can come out negative, but thats o.k. as the sum of all these  */
	/* quadraliterals should be exactly the area (up to correct sign).   */
	PolygonArea += (Vnext -> Coord[0] - V -> Coord[0]) *
			(Vnext -> Coord[1] + V -> Coord[1] - MinY) / 2.0;
	V = Vnext;
	Vnext = V -> Pnext;
    }
    while (V != VHead && V != NULL /* Should not happen! */);

    return ABS(PolygonArea);
}

/*****************************************************************************
*   Routine to evaluate the Volume of the given geom object, in object unit. *
*   This routine has a side effect that all non-convex polygons will be      *
* splitted to convex ones.						     *
* Algorithm (for each polygon, and let ObjMinY be the minimum OBJECT Y):     *
*								V3	     *
* 1. Set Polygon Area to be zero.			       /\	     *
*    Let V(0) be the first vertex, V(n) the last.	     /	  \ 	     *
*    For i goes from 1 to n-1 form triangles		   /	    \V2	     *
*    by V(0), V(i), V(i+1). For each such	       V4/	     |	     *
*    triangle do:					 \	     |	     *
*    1.1. Find the vertex (out of V(0), V(i), V(i+1))      \         |	     *
*         with the minimum Z - TriMinY.			     \_______|       *
*    1.2. The volume below V(0), V(i), V(i+1) triangle,	     V0      V1      *
*	  relative to ObjMinZ level, is the sum of:			     *
*	  1.2.1. volume of V'(0), V'(i), V'(i+1) - the			     *
*		 area of projection of V(0), V(i), V(i+1) on XY parallel     *
*		 plane, times (TriMinZ - ObjMinZ).			     *
*	  1.2.2. Assume V(0) is the one with the PolyMinZ. Let V"(i) and     *
*		 V"(i+1) be the projections of V(i) and V(i+1) on the plane  *
*		 Z = PolyZMin. The volume above 1.2.1. and below the polygon *
*		 (triangle!) will be: the area of quadraliteral V(i), V(i+1),*
*		 V"(i+1), V"(i), times distance of V(0) for quadraliteral    *
*		 plane divided by 3.					     *
*    1.3. If Z component of polygon normal is negative add 1.2. result to    *
*	  ObjectVolume, else subtract it.				     *
*****************************************************************************/
double PolyObjectVolume(IPObjectStruct *PObj)
{
    int PlaneExists;
    RealType ObjVolume = 0.0, ObjMinZ, TriMinZ, Area, PolygonVolume, Dist;
    PointType Pt1;
    PlaneType Plane;
    IPPolygonStruct *Pl;
    IPVertexStruct *V, *VHead, *Vnext, *Vtemp;

    if (!IP_IS_POLY_OBJ(PObj))
	IritFatalError("Geometric property requested on non polygonal object");

    if (IP_IS_POLYLINE_OBJ(PObj)) {
	return 0.0;
    }

    ObjMinZ = INFINITY;	 /* Find Object minimum Z value (used as min level). */
    Pl = PObj -> U.Pl;
    while (Pl != NULL) {
	V = VHead = Pl -> PVertex;
	do {
	    if (V -> Coord[2] < ObjMinZ)
		ObjMinZ = V -> Coord[2];
	    V = V -> Pnext;
	}
	while (V != VHead && V != NULL);
	Pl = Pl -> Pnext;
    }

    ConvexPolyObject(PObj);	       /* Make sure all polygons are convex. */
    Pl = PObj -> U.Pl;
    while (Pl != NULL) {
	PolygonVolume = 0.0; /* Volume below poly relative to ObjMinZ level. */
	/* We set VHead to be vertex with min Z: */
	V = Vtemp = VHead = Pl -> PVertex;
	do {
	    if (V -> Coord[2] < Vtemp -> Coord[2])
		Vtemp = V;
	    V = V -> Pnext;
	}
	while (V != VHead && V != NULL);
	VHead = Vtemp;	   /* Now VHead is the one with lowest Z in polygon! */
	TriMinZ = VHead -> Coord[2];	     /* Save this Z for fast access. */

	V = VHead -> Pnext;
	Vnext = V -> Pnext;
	do {
	    /* VHead, V, Vnext form the triangle - find volume 1.2.1. above: */
	    Area = Polygon3VrtxXYArea(VHead -> Coord, V -> Coord,
							       Vnext -> Coord);
	    PolygonVolume += Area * (TriMinZ - ObjMinZ);

	    /* VHead, V, Vnext form the triangle - find volume 1.2.2. above: */
	    Area = sqrt(SQR(V -> Coord[0] - Vnext -> Coord[0]) + /* XY dist. */
			SQR(V -> Coord[1] - Vnext -> Coord[1])) *
		   ((V -> Coord[2] + Vnext -> Coord[2]) / 2.0 - TriMinZ);
	    PT_COPY(Pt1, V -> Coord);
	    Pt1[2] = TriMinZ;
	    if ((PlaneExists =
		 CGPlaneFrom3Points(Plane, V -> Coord,
				    Vnext -> Coord, Pt1)) == 0) {
		/* Try second pt projected to Z = TriMinZ plane as third pt. */
		PT_COPY(Pt1, Vnext -> Coord);
		Pt1[2] = TriMinZ;
		PlaneExists =
			CGPlaneFrom3Points(Plane, V -> Coord, Vnext -> Coord,
									Pt1);
	    }
	    if (PlaneExists) {
		Dist = CGDistPointPlane(VHead -> Coord, Plane);
		PolygonVolume += Area * ABS(Dist) / 3.0;
	    }

	    V = Vnext;
	    Vnext = V -> Pnext;
	}
	while (Vnext != VHead);

	if (Pl -> Plane[2] < 0.0)
	    ObjVolume += PolygonVolume;
	else
	    ObjVolume -= PolygonVolume;

	Pl = Pl -> Pnext;
    }

    return ObjVolume;
}

/*****************************************************************************
*  Routine to evaluate the area of the given triangle projected to the XY    *
* plane, given as 3 Points. Only the X & Y components are considered.	     *
*  See PolyObjectArea above for algorithm:				     *
*****************************************************************************/
static RealType Polygon3VrtxXYArea(PointType Pt1, PointType Pt2, PointType Pt3)
{
    RealType PolygonArea = 0.0, MinY;

    MinY = MIN(Pt1[1], MIN(Pt2[1], Pt3[1])) * 2.0;

    PolygonArea += (Pt2[0] - Pt1[0]) * (Pt2[1] + Pt1[1] - MinY) / 2.0;
    PolygonArea += (Pt3[0] - Pt2[0]) * (Pt3[1] + Pt2[1] - MinY) / 2.0;
    PolygonArea += (Pt1[0] - Pt3[0]) * (Pt1[1] + Pt3[1] - MinY) / 2.0;

    return ABS(PolygonArea);
}

/*****************************************************************************
*  Routine to count number of polygons in given geometric object.	     *
*****************************************************************************/
double PolyCountPolys(IPObjectStruct *PObj)
{
    int Count = 0;
    IPPolygonStruct *Pl;

    if (!IP_IS_POLY_OBJ(PObj))
	IritFatalError("Geometric property requested on non polygonal object");

    Pl = PObj -> U.Pl;
    while (Pl != NULL) {
	Count++;
	Pl = Pl -> Pnext;
    }

    return ((double) Count);
}