/******************************************************************************
* 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);

/*****************************************************************************
* DESCRIPTION:                                                               M
*  Routine to evaluate the Area of the given geom. object, in object unit.   M
* Algorithm (for each polygon):					V3	     V
* 1. Set Polygon Area to be zero.			       /\	     V
*    Make a copy of the original polygon		     /	  \ 	     V
*    and transform it to a XY parallel plane.		   /	    \V2	     V
*    Find the minimum Y value of the polygon	       V4/	     |	     V
*    in the XY plane.					 \	     |	     V
* 2. Let V(0) be the first vertex, V(n) the last one.      \         |	     V
*    For i goes from 0 to n-1 add to Area the area   	     \_______|       V
*    below edge V(i), V(i+1):				     V0      V1      V
*    PolygonArea += (V(i+1)x - V(i)x) * (V(i+1)y' - V(i)y') / 2		     V
*    where V(i)y' is V(i)y - MinY, where MinY is polygon minimum Y value.    V
* 3. The result of step 2 is the area of the polygon itself. 		     V
*    However, it might be negative, so take the absolute result of step 2    V
*    and add it to the global ObjectArea.				     V
* Note step 2 is performed by another auxiliary routine: PolygonXYArea.      M
*                                                                            *
* PARAMETERS:                                                                M
*   PObj:     A polyhedra object to compute its surface area.                M
*                                                                            *
* RETURN VALUE:                                                              M
*   double:   The area of object PObj.                                       M
*                                                                            *
* KEYWORDS:                                                                  M
*   PolyObjectArea, area                                                     M
*****************************************************************************/
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;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Routine to evaluate the area of the given polygon projected on the XY    *
* plane.	                                                             *
*   The polygon does not have to be on a XY parallel plane, as only its XY   *
* projection is considered (Z is ignored).				     *
*   Returned the area of its XY parallel projection.			     *
*   See GeomObjectArea above for algorithm.				     *
*                                                                            *
* PARAMETERS:                                                                *
*   VHead:     Of vertex list to compute area for.                           *
*                                                                            *
* RETURN VALUE:                                                              *
*   RealType:  Computed area.                                                *
*****************************************************************************/
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) {
	if (MinY > V -> Coord[1])
	    MinY = V -> Coord[1];

	V = V -> Pnext;
    }

    if (V == NULL)
	V = VHead;
    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 == NULL ? VHead : V -> Pnext;
    }
    while (V != VHead && V != NULL);

    return ABS(PolygonArea);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to evaluate the Volume of the given geom object, in object unit. M
*   This routine has a side effect that all non-convex polygons will be      M
* splitted to convex ones.						     M
* Algorithm (for each polygon, and let ObjMinY be the minimum OBJECT Y):     M
*								V3	     V
* 1. Set Polygon Area to be zero.			       /\	     V
*    Let V(0) be the first vertex, V(n) the last.	     /	  \ 	     V
*    For i goes from 1 to n-1 form triangles		   /	    \V2	     V
*                   by V(0), V(i), V(i+1).             V4/	     |	     V
*    For each such triangle di:				 \	     |	     V
*    1.1. Find the vertex (out of V(0), V(i), V(i+1))      \         |	     V
*         with the minimum Z - TriMinY.			     \_______|       V
*    1.2. The volume below V(0), V(i), V(i+1) triangle,	     V0      V1      V
*	  relative to ObjMinZ level, is the sum of:			     V
*	  1.2.1. volume of V'(0), V'(i), V'(i+1) - the			     V
*		 area of projection of V(0), V(i), V(i+1) on XY parallel     V
*		 plane, times (TriMinZ - ObjMinZ).			     V
*	  1.2.2. Assume V(0) is the one with the PolyMinZ. Let V"(i) and     V
*		 V"(i+1) be the projections of V(i) and V(i+1) on the plane  V
*		 Z = PolyZMin. The volume above 1.2.1. and below the polygon V
*		 (triangle!) will be: the area of quadraliteral V(i), V(i+1),V
*		 V"(i+1), V"(i), times distance of V(0) for quadraliteral    V
*		 plane divided by 3.					     V
*    1.3. If Z component of polygon normal is negative add 1.2. result to    V
*	  ObjectVolume, else subtract it.				     V
*                                                                            *
* PARAMETERS:                                                                M
*   PObj:      To compute volume for.                                        M
*                                                                            *
* RETURN VALUE:                                                              M
*   double:    Computed volume.                                              M
*                                                                            *
* KEYWORDS:                                                                  M
*   PolyObjectVolume, volume                                                 M
*****************************************************************************/
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);
	if (V == NULL)
	    IritFatalError("For VOLUME polygons must be closed loops. Open loop detected.\n");

	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;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*  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.				     *
*                                                                            *
* PARAMETERS:                                                                *
*   Pt1, Pt2, Pt3: Vertices of triangle to compute area for.                 *
*                                                                            *
* RETURN VALUE:                                                              *
*   RealType:      Area of triangle in the XY plane.                         *
*****************************************************************************/
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);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*  Routine to count the number of polygons in the given geometric object.    M
*                                                                            *
* PARAMETERS:                                                                M
*   PObj:       To count number of polygons in.                              M
*                                                                            *
* RETURN VALUE:                                                              M
*   double:     Number of polygons, an integer value returned as double.     M
*                                                                            *
* KEYWORDS:                                                                  M
*   PolyCountPolys                                                           M
*****************************************************************************/
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);
}
