/*****************************************************************************
*   "Irit" - the 3d (not only polygonal) solid modeller.		     *
*									     *
* Written by:  Gershon Elber				Ver 0.2, Mar. 1990   *
******************************************************************************
*   Module to generate the geometric primitives defined in the system. The   *
* primitives currently defined are:					     *
* 1. BOX - main planes parallel box.					     *
* 2. GBOX - generalized box - 6 arbitrary planes.			     *
* 3. CYLIN - cylinder with any main direction.				     *
* 4. CONE, CONE2 - cone with any main direction (two bases).		     *
* 5. SPHERE								     *
* 6. TORUS - with any main direction.					     *
* 7. PLANE - non closed, single polygon object: circle with resolution edges *
* 8. POLY - directly define single polygon object by specifing its vertices. *
*   In addition, the following lower level operations are defined to create  *
* objects - EXTRUDE, and SURFREV, both require a polygon and a vector to     *
* extrude/rotate the polygon along.					     *
*****************************************************************************/

#include <stdio.h>
#include <math.h>
#include "irit_sm.h"
#include "geomat3d.h"
#include "allocate.h"
#include "attribut.h"
#include "convex.h"
#include "primitiv.h"

#define	MIN_RESOLUTION 4

static int GlblResolution = 16;

static IPPolygonStruct *GenInsidePoly(IPPolygonStruct *Pl);
static IPPolygonStruct *GenPolygon4Vrtx(VectorType V1,
					VectorType V2,
					VectorType V3,
					VectorType V4,
					VectorType Vin,
					IPPolygonStruct *Pnext);
static IPPolygonStruct *GenPolygon3Vrtx(VectorType V1,
					VectorType V2,
					VectorType V3,
					VectorType Vin,
					IPPolygonStruct *Pnext);
static void UpdateVertexNormal(NormalType Normal,
			       PointType Pt,
			       PointType InPt,
			       int Perpendicular,
			       PointType PerpPt);

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a BOX geometric object defined by Pt - the minimum     V
* 3d point, and Width - Dx Dy & Dz vector.		4		     V
*   Order of vertices is as                         5       7		     V
* follows in the picture:                           |   6   |		     V
*						    |   |   |		     V
* (Note vertex 0 is hidden behind edge 2-6)	    |	|   |		     V
*						    1   |   3                V
*							2		     V
*   All dimensions can be negative, denoting the reversed direction.         M
*                                                                            *
* PARAMETERS:                                                                M
*   Pt:          Low end corner of BOX.                                      M
*   WidthX:      Width of BOX (X axis).                                      M
*   WidthY:      Depth of BOX( Y axis).                                      M
*   WidthZ:      Height of BOX( Z axis).                                     M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:   A BOX primitive                                      M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenBOXObject, box, primitives                                        M
*****************************************************************************/
IPObjectStruct *PrimGenBOXObject(VectorType Pt,
				 RealType WidthX,
				 RealType WidthY,
				 RealType WidthZ)
{
    VectorType Dir1, Dir2, Dir3;

    PT_CLEAR(Dir1);	Dir1[0] = WidthX;      /* Prepare direction vectors. */
    PT_CLEAR(Dir2);	Dir2[1] = WidthY;	   /* Parallel to main axes. */
    PT_CLEAR(Dir3);	Dir3[2] = WidthZ;		   /* For GBOX call. */

    return PrimGenGBOXObject(Pt, Dir1, Dir2, Dir3);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a GBOX geometric object defined by Pt - the minimum    M
* 3d point, and 3 direction Vectors Dir1, Dir2, Dir3.                        M
*   If two of the direction vectors are parallel the GBOX degenerates to     M
* a zero volume object. A NULL pointer is returned in that case.   	     M
* 							4		     V
* Order of vertices is as                           5       7		     V
* follows in the picture:                           |   6   |		     V
*						    |   |   |		     V
* (Note vertex 0 is hidden behind edge 2-6)	    |	|   |		     V
*						    1   |   3                V
*							2		     V
*                                                                            *
* PARAMETERS:                                                                M
*   Pt:                Low end corner of GBOX.				     M
*   Dir1, Dir2, Dir3:  Three independent directional vectors to define GBOX. M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:  A GBOX primitive.                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenGBOXObject, general box, box, primitives                          M
*****************************************************************************/
IPObjectStruct *PrimGenGBOXObject(VectorType Pt,
				  VectorType Dir1,
				  VectorType Dir2,
				  VectorType Dir3)
{
    int i;
    VectorType Temp;
    VectorType V[8];				  /* Hold 8 vertices of BOX. */
    IPVertexStruct *PVertex;
    IPPolygonStruct *PPolygon;
    IPObjectStruct *PBox;

    GMVecCrossProd(Temp, Dir1, Dir2);
    if (APX_EQ(PT_LENGTH(Temp), 0.0))
	return NULL;
    GMVecCrossProd(Temp, Dir2, Dir3);
    if (APX_EQ(PT_LENGTH(Temp), 0.0))
	return NULL;
    GMVecCrossProd(Temp, Dir3, Dir1);
    if (APX_EQ(PT_LENGTH(Temp), 0.0))
	return NULL;

    /* Also the 0..7 sequence is binary decoded such that bit 0 is Dir1, */
    /* bit 1 Dir2, and bit 2 is Dir3 increment:				 */
    for (i = 0; i < 8; i++) {
	PT_COPY(V[i], Pt);

	if (i & 1)
	    PT_ADD(V[i], V[i], Dir1);
	if (i & 2)
	    PT_ADD(V[i], V[i], Dir2);
	if (i & 4)
	    PT_ADD(V[i], V[i], Dir3);
    }

    PBox = GenPolyObject("", NULL, NULL); /* Generate the BOX object itself: */

    /* And generate the 6 polygons (Bottom, top and 4 sides in this order):  */
    PBox -> U.Pl = GenPolygon4Vrtx(V[0], V[1], V[3], V[2], V[4], PBox -> U.Pl);
    PBox -> U.Pl = GenPolygon4Vrtx(V[6], V[7], V[5], V[4], V[0], PBox -> U.Pl);
    PBox -> U.Pl = GenPolygon4Vrtx(V[4], V[5], V[1], V[0], V[2], PBox -> U.Pl);
    PBox -> U.Pl = GenPolygon4Vrtx(V[5], V[7], V[3], V[1], V[0], PBox -> U.Pl);
    PBox -> U.Pl = GenPolygon4Vrtx(V[7], V[6], V[2], V[3], V[1], PBox -> U.Pl);
    PBox -> U.Pl = GenPolygon4Vrtx(V[6], V[4], V[0], V[2], V[3], PBox -> U.Pl);

    /* Update the vertices normals using the polygon plane equation: */
    for (PPolygon = PBox -> U.Pl;
	 PPolygon != NULL;
	 PPolygon = PPolygon -> Pnext) {
	PVertex = PPolygon -> PVertex;
	do {
	    PT_COPY(PVertex -> Normal, PPolygon -> Plane);
	    PVertex = PVertex -> Pnext;
	}
	while (PVertex != PPolygon -> PVertex);
    }

    return PBox;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a CONE geometric object defined by Pt - the base       M
* 3d center point, Dir - the cone direction and height, and base radius R.   M
*   See also PrimSetResolution on fineness control of approximation of the   M
* primitive using flat faces.                                                M
*                                                                            *
* PARAMETERS:                                                                M
*   Pt:         Center location of Base of CONE.                             M
*   Dir:        Direction and distance from Pt to apex of CONE.              M
*   R:          Radius of Base of the cone.                                  M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:    A CONE Primitive.                                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenCONEObject, cone, primitives                                      M
*****************************************************************************/
IPObjectStruct *PrimGenCONEObject(VectorType Pt, VectorType Dir, RealType R)
{
    int i;
    RealType Angle, AngleStep;
    PointType LastCirclePt, CirclePt, ApexPt;
    NormalType LastCircleNrml, CircleNrml, ApexNrml;
    MatrixType Mat;
    IPVertexStruct *VBase, *PVertex;
    IPPolygonStruct *PBase;
    IPObjectStruct *PCone;

    CGGenTransMatrixZ2Dir(Mat, Pt, Dir, R);   /* Transform from unit circle. */

    PT_COPY(ApexPt, Pt);		   /* Find the apex point: Pt + Dir. */
    PT_ADD(ApexPt, ApexPt, Dir);
    PT_NORMALIZE(Dir);

    PCone = GenPolyObject("", NULL, NULL);   /* Gen. the CONE object itself: */
    /* Also allocate the base polygon header with first vertex on it: */
    PBase = IPAllocPolygon(0, 0, VBase = IPAllocVertex(0, 0, NULL, NULL), NULL);

    LastCirclePt[0] = 1.0;		/* First point is allways Angle = 0. */
    LastCirclePt[1] = 0.0;
    LastCirclePt[2] = 0.0;
    MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);

    UpdateVertexNormal(LastCircleNrml, LastCirclePt, Pt, TRUE, ApexPt);

    PT_COPY(VBase -> Coord, LastCirclePt);/* Update first pt in base polygon.*/
    PT_COPY(VBase -> Normal, Dir);

    AngleStep = M_PI * 2 / GlblResolution;

    for (i = 1; i <= GlblResolution; i++) {   /* Pass the whole base circle. */
	Angle = AngleStep * i;		     /* Prevent from additive error. */

	CirclePt[0] = cos(Angle);
	CirclePt[1] = sin(Angle);
	CirclePt[2] = 0.0;
	MatMultVecby4by4(CirclePt, CirclePt, Mat);

 	UpdateVertexNormal(CircleNrml, CirclePt, Pt, TRUE, ApexPt);

	PCone -> U.Pl = GenPolygon3Vrtx(LastCirclePt, ApexPt,
					      CirclePt, Pt, PCone -> U.Pl);

	/* Update the normals for this cone side polygon vertices: */
	PVertex = PCone -> U.Pl -> PVertex;
	PT_COPY(PVertex -> Normal, LastCircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);
	PVertex = PVertex -> Pnext;
	/* The apex normal is the average of the two base vertices: */
	PT_ADD(ApexNrml, CircleNrml, LastCircleNrml);
	PT_NORMALIZE(ApexNrml);
	PT_COPY(PVertex -> Normal, ApexNrml);
	IP_SET_NORMAL_VRTX(PVertex);
	PVertex = PVertex -> Pnext;
	PT_COPY(PVertex -> Normal, CircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);

	/* And add this vertex to base polygon: */
	if (i == GlblResolution)       /* Its last point - make it circular. */
	    VBase -> Pnext = PBase -> PVertex;
	else {
	    VBase -> Pnext = IPAllocVertex(0, 0, NULL, NULL);
	    VBase = VBase -> Pnext;
	    PT_COPY(VBase -> Normal, Dir);
	    PT_COPY(VBase -> Coord, CirclePt);
	}

	PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
	PT_COPY(LastCircleNrml, CircleNrml);
    }

    /* Update base polygon plane equation. */
    IritPrsrUpdatePolyPlane2(PBase, ApexPt);
    IP_SET_CONVEX_POLY(PBase);		       /* Mark it as convex polygon. */
    PBase -> Pnext = PCone -> U.Pl;  /* And stick it into the cone polygons. */
    PCone -> U.Pl = PBase;

    return PCone;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a truncated CONE, CON2, geometric object defined by    M
* Pt - the base 3d center point, Dir - the cone direction and height, and    M
* two base radii R1 and R2.						     M
*   See also PrimSetResolution on fineness control of approximation of the   M
* primitive using flat faces.                                                M
*                                                                            *
* PARAMETERS:                                                                M
*   Pt:      Center location of Base of CON2.                                M
*   Dir:     Direction and distance from Pt to center of other base of CON2. M
*   R1, R2:  Two base radii of the truncated CON2                            M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:    A CON2 Primitive.                                   M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenCONE2Object, cone, primitives                                     M
*****************************************************************************/
IPObjectStruct *PrimGenCONE2Object(VectorType Pt,
				   VectorType Dir,
				   RealType R1,
				   RealType R2)
{
    int i;
    RealType Angle, AngleStep;
    PointType LastCirclePt, CirclePt, ApexPt, LastApexPt1, ApexPt1;
    NormalType LastCircleNrml, CircleNrml;
    VectorType InvDir;
    MatrixType Mat1, Mat2;
    IPVertexStruct *VBase1, *VBase2, *PVertex;
    IPPolygonStruct *PBase1, *PBase2;
    IPObjectStruct *PCone;

    PT_COPY(ApexPt, Pt);		   /* Find the apex point: Pt + Dir. */
    PT_ADD(ApexPt, ApexPt, Dir);
    PT_NORMALIZE(Dir);
    PT_COPY(InvDir, Dir);
    PT_SCALE(InvDir, -1.0);

    CGGenTransMatrixZ2Dir(Mat1, Pt, Dir, R1);    /* Trans. from unit circle. */
    CGGenTransMatrixZ2Dir(Mat2, ApexPt, Dir, R2);

    PCone = GenPolyObject("", NULL, NULL);   /* Gen. the CONE object itself: */
    /* Also allocate the base polygon header with first vertex on it: */
    PBase1 = IPAllocPolygon(0, 0,
			    VBase1 = IPAllocVertex(0, 0, NULL, NULL), NULL);
    PBase2 = IPAllocPolygon(0, 0,
			    VBase2 = IPAllocVertex(0, 0, NULL, NULL), NULL);

    /* First point is allways at Angle = 0. */
    LastCirclePt[0] = LastApexPt1[0] = 1.0;
    LastCirclePt[1] = LastApexPt1[1] = 0.0;
    LastCirclePt[2] = LastApexPt1[2] = 0.0;
    MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat1);
    MatMultVecby4by4(LastApexPt1, LastApexPt1, Mat2);

    UpdateVertexNormal(LastCircleNrml, LastCirclePt, Pt, TRUE, ApexPt);

    PT_COPY(VBase1 -> Coord, LastCirclePt);/* Update 1st pt in base1 polygon.*/
    PT_COPY(VBase1 -> Normal, Dir);
    PT_COPY(VBase2 -> Coord, LastApexPt1);/* Update 1st pt in base2 polygon. */
    PT_COPY(VBase2 -> Normal, InvDir);

    AngleStep = M_PI * 2 / GlblResolution;

    for (i = 1; i <= GlblResolution; i++) {   /* Pass the whole base circle. */
	Angle = AngleStep * i;		     /* Prevent from additive error. */

	CirclePt[0] = ApexPt1[0] = cos(Angle);
	CirclePt[1] = ApexPt1[1] = sin(Angle);
	CirclePt[2] = ApexPt1[2] = 0.0;
	MatMultVecby4by4(CirclePt, CirclePt, Mat1);
	MatMultVecby4by4(ApexPt1, ApexPt1, Mat2);

 	UpdateVertexNormal(CircleNrml, CirclePt, Pt, TRUE, ApexPt);

	PCone -> U.Pl = GenPolygon4Vrtx(LastCirclePt, LastApexPt1, ApexPt1,
					  CirclePt, Pt, PCone -> U.Pl);

	/* Update the normals for this cone side polygon vertices: */
	PVertex = PCone -> U.Pl -> PVertex;
	PT_COPY(PVertex -> Normal, LastCircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);
	PVertex = PVertex -> Pnext;
	PT_COPY(PVertex -> Normal, LastCircleNrml );
	IP_SET_NORMAL_VRTX(PVertex);
	PVertex = PVertex -> Pnext;
	PT_COPY(PVertex -> Normal, CircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);
	PVertex = PVertex -> Pnext;
	PT_COPY(PVertex -> Normal, CircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);

	/* And add these vertices to base polygons: */
	if (i == GlblResolution) {     /* Its last point - make it circular. */
	    VBase1 -> Pnext = PBase1 -> PVertex;
	    VBase2 -> Pnext = PBase2 -> PVertex;
	}
	else {
	    VBase1 -> Pnext = IPAllocVertex(0, 0, NULL, NULL);
	    VBase1 = VBase1 -> Pnext;
	    PT_COPY(VBase1 -> Coord, CirclePt);
	    PT_COPY(VBase1 -> Normal, Dir);
	    VBase2 -> Pnext = IPAllocVertex(0, 0, NULL, NULL);
	    VBase2 = VBase2 -> Pnext;
	    PT_COPY(VBase2 -> Coord, ApexPt1);
	    PT_COPY(VBase2 -> Normal, InvDir);
	}

	PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
	PT_COPY(LastApexPt1, ApexPt1);
	PT_COPY(LastCircleNrml, CircleNrml);
    }

    /* Update base polygon plane equation. */
    IritPrsrUpdatePolyPlane2(PBase1, ApexPt);
    IP_SET_CONVEX_POLY(PBase1);		       /* Mark it as convex polygon. */
    /* Update base polygon plane equation. */
    IritPrsrUpdatePolyPlane2(PBase2, Pt);
    IP_SET_CONVEX_POLY(PBase2);		       /* Mark it as convex polygon. */

    PBase1 -> Pnext = PCone -> U.Pl;    /* And stick into the cone polygons. */
    PCone -> U.Pl = PBase1;
    PBase2 -> Pnext = PCone -> U.Pl;
    PCone -> U.Pl = PBase2;

    return PCone;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a CYLINder geometric object defined by Pt - the base   M
* 3d center point, Dir - the cylinder direction and height, and radius R.    M
*   See also PrimSetResolution on fineness control of approximation of the   M
* primitive using flat faces.                                                M
*                                                                            *
* PARAMETERS:                                                                M
*   Pt:         Center location of Base of CYLINder.                         M
*   Dir:        Direction and distance from Pt to other base of cylinder.    M
*   R:          Radius of Base of the cylinder.                              M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:    A CYLINDER Primitive.                               M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenCYLINObject, cylinder, primitives                                 M
*****************************************************************************/
IPObjectStruct *PrimGenCYLINObject(VectorType Pt, VectorType Dir, RealType R)
{
    int i;
    RealType Angle, AngleStep;
    PointType LastCirclePt, CirclePt, TLastCirclePt, TCirclePt, TPt, Dummy;
    VectorType ForwardDir, BackwardDir;
    NormalType LastCircleNrml, CircleNrml;
    MatrixType Mat;
    IPVertexStruct *VBase1, *VBase2, *PVertex;
    IPPolygonStruct *PBase1, *PBase2;
    IPObjectStruct *PCylin;

    CGGenTransMatrixZ2Dir(Mat, Pt, Dir, R);   /* Transform from unit circle. */

    PCylin = GenPolyObject("", NULL, NULL); /* Gen. the CYLIN object itself: */
    /* Also allocate the bases polygon header with first vertex on it: */
    PBase1 = IPAllocPolygon(0, 0,
			    VBase1 = IPAllocVertex(0, 0, NULL, NULL), NULL);
    PBase2 = IPAllocPolygon(0, 0,
			    VBase2 = IPAllocVertex(0, 0, NULL, NULL), NULL);

    PT_ADD(TPt, Pt, Dir);	       /* Translated circle center (by Dir). */

    /* Prepare the normal directions for the two bases: */
    PT_COPY(ForwardDir, Dir);
    PT_NORMALIZE(ForwardDir);
    PT_COPY(BackwardDir, ForwardDir);
    PT_SCALE(BackwardDir, -1.0);

    LastCirclePt[0] = 1.0;		/* First point is allways Angle = 0. */
    LastCirclePt[1] = 0.0;
    LastCirclePt[2] = 0.0;
    MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);

    UpdateVertexNormal(LastCircleNrml, LastCirclePt, Pt, FALSE, Dummy);

    PT_COPY(VBase1 -> Coord, LastCirclePt);/* Update 1st pt in base1 polygon.*/
    PT_COPY(VBase1 -> Normal, ForwardDir);
    PT_ADD(TLastCirclePt, LastCirclePt, Dir); /* Translated circle (by Dir). */
    PT_COPY(VBase2 -> Coord, TLastCirclePt);/*Update 1st pt in base2 polygon.*/
    PT_COPY(VBase2 -> Normal, BackwardDir);

    AngleStep = M_PI * 2 / GlblResolution;

    for (i = 1; i <= GlblResolution; i++) {   /* Pass the whole base circle. */
	Angle = AngleStep * i;		     /* Prevent from additive error. */

	CirclePt[0] = cos(Angle);
	CirclePt[1] = sin(Angle);
	CirclePt[2] = 0.0;
	MatMultVecby4by4(CirclePt, CirclePt, Mat);

	UpdateVertexNormal(CircleNrml, CirclePt, Pt, FALSE, Dummy);

	PT_ADD(TCirclePt, CirclePt, Dir);     /* Translated circle (by Dir). */

	PCylin -> U.Pl = GenPolygon4Vrtx(TLastCirclePt, TCirclePt, CirclePt,
					 LastCirclePt, Pt, PCylin -> U.Pl);
	/* Update the normals for this cylinder side polygon vertices: */
	PVertex = PCylin -> U.Pl -> PVertex;
	PT_COPY(PVertex -> Normal, LastCircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);
	PVertex = PVertex -> Pnext;
	PT_COPY(PVertex -> Normal, CircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);
	PVertex = PVertex -> Pnext;
	PT_COPY(PVertex -> Normal, CircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);
	PVertex = PVertex -> Pnext;
	PT_COPY(PVertex -> Normal, LastCircleNrml);
	IP_SET_NORMAL_VRTX(PVertex);

	/* And add this vertices to the two cylinder bases: 		     */
	/* Note Base1 is build forward, while Base2 is build backward so it  */
	/* will be consistent - cross product of 2 consecutive edges will    */
	/* point into the model.					     */
	if (i == GlblResolution) {     /* Its last point - make it circular. */
	    VBase1 -> Pnext = PBase1 -> PVertex;
	    VBase2 -> Pnext = PBase2 -> PVertex;
	}
	else {
	    VBase1 -> Pnext = IPAllocVertex(0, 0, NULL, NULL);
	    VBase1 = VBase1 -> Pnext;
	    PT_COPY(VBase1 -> Coord, CirclePt);
	    PT_COPY(VBase1 -> Normal, ForwardDir);
	    PBase2 -> PVertex = IPAllocVertex(0, 0, NULL, PBase2 -> PVertex);
	    PT_COPY(PBase2 -> PVertex -> Coord, TCirclePt);
	    PT_COPY(PBase2 -> PVertex -> Normal, BackwardDir);
	}

	PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
	PT_COPY(TLastCirclePt, TCirclePt);
	PT_COPY(LastCircleNrml, CircleNrml);
    }

    /* Update base polygon plane equation. */
    IritPrsrUpdatePolyPlane2(PBase1, TPt);
    IP_SET_CONVEX_POLY(PBase1);		       /* Mark it as convex polygon. */
    PBase1 -> Pnext = PCylin -> U.Pl;   /* And stick it into cylin polygons. */
    PCylin -> U.Pl = PBase1;
    /* Update base polygon plane equation. */
    IritPrsrUpdatePolyPlane2(PBase2, Pt);
    IP_SET_CONVEX_POLY(PBase2);		       /* Mark it as convex polygon. */
    PBase2 -> Pnext = PCylin -> U.Pl;   /* And stick it into cylin polygons. */
    PCylin -> U.Pl = PBase2;

    return PCylin;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a SPHERE geometric object defined by Center, the       M
* center of the sphere and R, its radius.				     M
*   See also PrimSetResolution on fineness control of approximation of the   M
* primitive using flat faces.                                                M
*                                                                            *
* PARAMETERS:                                                                M
*   Center:   Center location of SPHERE.                                     M
*   R         Radius of sphere.                                              M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:    A SPHERE Primitive.                                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenSPHEREObject, sphere, primitives                                  M
*****************************************************************************/
IPObjectStruct *PrimGenSPHEREObject(VectorType Center, RealType R)
{
    int i, j, k;
    RealType TetaAngle, TetaAngleStep, FeeAngle, FeeAngleStep,
	CosFeeAngle1, SinFeeAngle1, CosFeeAngle2, SinFeeAngle2;
    PointType LastCircleLastPt, LastCirclePt, CirclePt, CircleLastPt, Dummy;
    IPVertexStruct *PVertex;
    IPObjectStruct *PSphere;

    PSphere = GenPolyObject("", NULL, NULL);/* Gen the SPHERE object itself: */

    TetaAngleStep = M_PI * 2.0 / GlblResolution;     /* Runs from 0 to 2*PI. */
    FeeAngleStep = M_PI * 2.0 / GlblResolution;	/* Runs from -PI/2 yo +PI/2. */

    /* Generate the lowest (south pole) triangular polygons: */
    FeeAngle = (-M_PI/2.0) + FeeAngleStep; /* First circle above south pole. */
    CosFeeAngle1 = cos(FeeAngle) * R;
    SinFeeAngle1 = sin(FeeAngle) * R;
    PT_COPY(LastCirclePt, Center);		/* Calculate the south pole. */
    LastCirclePt[2] -= R;
    PT_COPY(CircleLastPt, Center);    /* Calc. last point on current circle. */
    CircleLastPt[0] += CosFeeAngle1;
    CircleLastPt[2] += SinFeeAngle1;

    for (i = 1; i <= GlblResolution; i++) {/* Pass whole (horizontal) circle.*/
	TetaAngle = TetaAngleStep * i;	     /* Prevent from additive error. */

	PT_COPY(CirclePt, Center); /* Calc. current point on current circle. */
	CirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
	CirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
	CirclePt[2] += SinFeeAngle1;

	PSphere -> U.Pl = GenPolygon3Vrtx(LastCirclePt, CircleLastPt,
					CirclePt, Center, PSphere -> U.Pl);
	/* Update normals: */
	for (j = 0, PVertex = PSphere -> U.Pl -> PVertex;
	     j < 3;
	     j++, PVertex = PVertex -> Pnext) {
	    UpdateVertexNormal(PVertex -> Normal, PVertex -> Coord, Center,
								FALSE, Dummy);
	    IP_SET_NORMAL_VRTX(PVertex);
	}

	PT_COPY(CircleLastPt, CirclePt);/* Save pt in last pt for next time. */
    }

    /* Generate the middle rectangular polygons: */
    for (i = 1; i < GlblResolution / 2 - 1; i++) {/* For all horiz. circles. */
	FeeAngle = (-M_PI/2.0) + FeeAngleStep * i;
	CosFeeAngle1 = cos(FeeAngle) * R;
	SinFeeAngle1 = sin(FeeAngle) * R;
	FeeAngle = (-M_PI/2.0) + FeeAngleStep * (i + 1);
	CosFeeAngle2 = cos(FeeAngle) * R;
	SinFeeAngle2 = sin(FeeAngle) * R;
	PT_COPY(CircleLastPt, Center);/* Calc. last point on current circle. */
	CircleLastPt[0] += CosFeeAngle2;
	CircleLastPt[2] += SinFeeAngle2;
	PT_COPY(LastCircleLastPt, Center);/* Calc. last point on last circle.*/
	LastCircleLastPt[0] += CosFeeAngle1;
	LastCircleLastPt[2] += SinFeeAngle1;

	for (j = 1; j <= GlblResolution; j++) {/* Pass whole (horiz.) circle.*/
	    TetaAngle = TetaAngleStep * j;   /* Prevent from additive error. */

	    PT_COPY(CirclePt, Center);/* Calc. current pt on current circle. */
	    CirclePt[0] += cos(TetaAngle) * CosFeeAngle2;
	    CirclePt[1] += sin(TetaAngle) * CosFeeAngle2;
	    CirclePt[2] += SinFeeAngle2;
	    PT_COPY(LastCirclePt, Center);/* Calc. current pt on last circle.*/
	    LastCirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
	    LastCirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
	    LastCirclePt[2] += SinFeeAngle1;

	    PSphere -> U.Pl = GenPolygon4Vrtx(LastCirclePt, LastCircleLastPt,
			CircleLastPt, CirclePt, Center, PSphere -> U.Pl);
	    /* Update normals: */
	    for (k = 0, PVertex = PSphere -> U.Pl -> PVertex;
		 k < 4;
		 k++, PVertex = PVertex -> Pnext) {
		UpdateVertexNormal(PVertex -> Normal, PVertex -> Coord, Center,
								FALSE, Dummy);
		IP_SET_NORMAL_VRTX(PVertex);
	    }

	    PT_COPY(CircleLastPt, CirclePt);	      /* Save pt in last pt. */
	    PT_COPY(LastCircleLastPt, LastCirclePt);
	}
    }

    /* Generate the upper most (north pole) triangular polygons: */
    FeeAngle = (M_PI/2.0) - FeeAngleStep;  /* First circle below north pole. */
    CosFeeAngle1 = cos(FeeAngle) * R;
    SinFeeAngle1 = sin(FeeAngle) * R;
    PT_COPY(LastCirclePt, Center);		/* Calculate the north pole. */
    LastCirclePt[2] += R;
    PT_COPY(CircleLastPt, Center);    /* Calc. last point on current circle. */
    CircleLastPt[0] += CosFeeAngle1;
    CircleLastPt[2] += SinFeeAngle1;

    for (i = 1; i <= GlblResolution; i++) {/* Pass whole (horizontal) circle.*/
	TetaAngle = TetaAngleStep * i;	     /* Prevent from additive error. */

	PT_COPY(CirclePt, Center); /* Calc. current point on current circle. */
	CirclePt[0] += cos(TetaAngle) * CosFeeAngle1;
	CirclePt[1] += sin(TetaAngle) * CosFeeAngle1;
	CirclePt[2] += SinFeeAngle1;

	PSphere -> U.Pl = GenPolygon3Vrtx(LastCirclePt, CirclePt, CircleLastPt,
					    Center, PSphere -> U.Pl);

	/* Update normals: */
	for (j = 0, PVertex = PSphere -> U.Pl -> PVertex;
	     j < 3;
	     j++, PVertex = PVertex -> Pnext) {
	    UpdateVertexNormal(PVertex -> Normal, PVertex -> Coord, Center,
								FALSE, Dummy);
	    IP_SET_NORMAL_VRTX(PVertex);
	}

	PT_COPY(CircleLastPt, CirclePt);/* Save pt in last pt for next time. */
    }

    return PSphere;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a TORUS geometric object defined by Center - torus 3d  M
* center point, the main torus plane normal Normal, major radius Rmajor and  M
* minor radius Rminor (Tube radius).					     M
*   Teta runs on the major circle, Fee on the minor one. Then                M
* X = (Rmajor + Rminor * cos(Fee)) * cos(Teta)				     V
* Y = (Rmajor + Rminor * cos(Fee)) * sin(Teta)				     V
* Z = Rminor * sin(Fee)							     V
*   See also PrimSetResolution on fineness control of approximation of the   M
* primitive using flat faces.                                                M
*                                                                            *
* PARAMETERS:                                                                M
*   Center:  Center location of the TORUS primitive.                         M
*   Normal:  Normal to the major plane of the torus.                         M
*   Rmajor:  Major radius of torus.                                          M
*   Rminor:  Minor radius of torus.                                          M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:    A TOURS Primitive.                                  M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenTORUSObject, torus, primitives                                    M
*****************************************************************************/
IPObjectStruct *PrimGenTORUSObject(VectorType Center,
				   VectorType Normal,
				   RealType Rmajor,
				   RealType Rminor)
{
    int i, j;
    RealType TetaAngle, TetaAngleStep, FeeAngle, FeeAngleStep,
	CosFeeAngle1, SinFeeAngle1, CosFeeAngle2, SinFeeAngle2;
    PointType LastCircleLastPt, LastCirclePt, CirclePt, CircleLastPt,
	LastInPt, InPt, Dummy;
    MatrixType Mat;
    IPVertexStruct *PVertex;
    IPObjectStruct *PTorus;

    CGGenTransMatrixZ2Dir(Mat, Center, Normal, 1.0);/* Trans from unit circ. */

    PTorus = GenPolyObject("", NULL, NULL); /* Gen. the Torus object itself: */

    TetaAngleStep = M_PI * 2.0 / GlblResolution;     /* Runs from 0 to 2*PI. */
    FeeAngleStep = M_PI * 2.0 / GlblResolution;	     /* Runs from 0 to 2*PI. */

    for (i = 1; i <= GlblResolution; i++) {
	FeeAngle = FeeAngleStep * (i - 1);
	CosFeeAngle1 = cos(FeeAngle) * Rminor;
	SinFeeAngle1 = sin(FeeAngle) * Rminor;
	FeeAngle = FeeAngleStep * i;
	CosFeeAngle2 = cos(FeeAngle) * Rminor;
	SinFeeAngle2 = sin(FeeAngle) * Rminor;
	LastCircleLastPt[0] = Rmajor + CosFeeAngle1;
	LastCircleLastPt[1] = 0.0;
	LastCircleLastPt[2] = SinFeeAngle1;
	MatMultVecby4by4(LastCircleLastPt, LastCircleLastPt, Mat);
	LastCirclePt[0] = Rmajor + CosFeeAngle2;
	LastCirclePt[1] = 0.0;
	LastCirclePt[2] = SinFeeAngle2;
	MatMultVecby4by4(LastCirclePt, LastCirclePt, Mat);
	/* Point inside the object relative to this polygon: */
	LastInPt[0] = Rmajor;
	LastInPt[1] = 0.0;
	LastInPt[2] = 0.0;
	MatMultVecby4by4(LastInPt, LastInPt, Mat);

	for (j = 1; j <= GlblResolution; j++) {
	    TetaAngle = TetaAngleStep * j;   /* Prevent from additive error. */

	    CircleLastPt[0] = (Rmajor + CosFeeAngle1) * cos(TetaAngle);
	    CircleLastPt[1] = (Rmajor + CosFeeAngle1) * sin(TetaAngle);
	    CircleLastPt[2] = SinFeeAngle1;
	    MatMultVecby4by4(CircleLastPt, CircleLastPt, Mat);
	    CirclePt[0] = (Rmajor + CosFeeAngle2) * cos(TetaAngle);
	    CirclePt[1] = (Rmajor + CosFeeAngle2) * sin(TetaAngle);
	    CirclePt[2] = SinFeeAngle2;
	    MatMultVecby4by4(CirclePt, CirclePt, Mat);
	    /* Point inside the object relative to this polygon: */
	    InPt[0] = Rmajor * cos(TetaAngle);
	    InPt[1] = Rmajor * sin(TetaAngle);
	    InPt[2] = 0.0;
	    MatMultVecby4by4(InPt, InPt, Mat);

	    PTorus -> U.Pl = GenPolygon4Vrtx(CirclePt, CircleLastPt,
		LastCircleLastPt, LastCirclePt, InPt, PTorus -> U.Pl);

	    /* Update normals: */
	    PVertex = PTorus -> U.Pl -> PVertex;
	    UpdateVertexNormal(PVertex -> Normal, PVertex -> Coord, InPt,
								FALSE, Dummy);
	    IP_SET_NORMAL_VRTX(PVertex);
	    PVertex = PVertex -> Pnext;
	    UpdateVertexNormal(PVertex -> Normal, PVertex -> Coord, InPt,
								FALSE, Dummy);
	    IP_SET_NORMAL_VRTX(PVertex);
	    PVertex = PVertex -> Pnext;
	    UpdateVertexNormal(PVertex -> Normal, PVertex -> Coord, LastInPt,
								FALSE, Dummy);
	    IP_SET_NORMAL_VRTX(PVertex);
	    PVertex = PVertex -> Pnext;
	    UpdateVertexNormal(PVertex -> Normal, PVertex -> Coord, LastInPt,
								FALSE, Dummy);
	    IP_SET_NORMAL_VRTX(PVertex);

	    PT_COPY(LastCirclePt, CirclePt);	      /* Save pt in last pt. */
	    PT_COPY(LastCircleLastPt, CircleLastPt);
	    PT_COPY(LastInPt, InPt);
	}
    }

    return PTorus;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a POLYDISK geometric object defined by the normal N    M
* and the translation vector T. The object is a planar disk (a circle of     M
* GlblResolution points in it...) and its radius is equal to R.		     M
*   The normal direction is assumed to point to the inside of the object.    M
*   See also PrimSetResolution on fineness control of approximation of the   M
* primitive using flat faces.                                                M
*                                                                            *
* PARAMETERS:                                                                M
*   N:         Normal to the plane this disk included in.                    M
*   T:         A translation factor of the center of the disk.               M
*   R:         Radius of teh disk.                                           M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:   A single polygon object - a disk.                    M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenPOLYDISKObject, disk, primitives                                  M
*****************************************************************************/
IPObjectStruct *PrimGenPOLYDISKObject(VectorType N, VectorType T, RealType R)
{
    int i;
    RealType Angle, AngleStep;
    PointType CirclePt;
    MatrixType Mat;
    IPVertexStruct *V;
    IPPolygonStruct *PCirc;
    IPObjectStruct *PPDisk;

    CGGenTransMatrixZ2Dir(Mat, T, N, R);      /* Transform from unit circle. */
    PT_NORMALIZE(N);

    PPDisk = GenPolyObject("", NULL, NULL); /* Gen. the PLANE object itself: */
    PCirc = IPAllocPolygon(0, 0, V = IPAllocVertex(0, 0, NULL, NULL), NULL);
    PPDisk -> U.Pl = PCirc;

    CirclePt[0] = 1.0;			/* First point is allways Angle = 0. */
    CirclePt[1] = 0.0;
    CirclePt[2] = 0.0;
    MatMultVecby4by4(CirclePt, CirclePt, Mat);
    PT_COPY(V -> Coord, CirclePt);     /* Update first pt in circle polygon. */
    PT_COPY(V -> Normal, N);

    AngleStep = M_PI * 2 / GlblResolution;

    for (i = 1; i <= GlblResolution; i++) {   /* Pass the whole base circle. */
	Angle = AngleStep * i;		     /* Prevent from additive error. */

	CirclePt[0] = cos(Angle);
	CirclePt[1] = sin(Angle);
	CirclePt[2] = 0.0;

	MatMultVecby4by4(CirclePt, CirclePt, Mat);

	/* And add this vertices to the two cylinder bases: */
	if (i == GlblResolution) {     /* Its last point - make it circular. */
	    V -> Pnext = PCirc -> PVertex;
	}
	else {
	    V -> Pnext = IPAllocVertex(0, 0, NULL, NULL);
	    V = V -> Pnext;
	    PT_COPY(V -> Coord, CirclePt);
	    PT_COPY(V -> Normal, N);
	}
    }

    PT_ADD(CirclePt, CirclePt, N);   /* Make a point "IN" the circle object. */
    /* Update base polygon plane equation. */
    IritPrsrUpdatePolyPlane2(PCirc, CirclePt);
    IP_SET_CONVEX_POLY(PCirc);		       /* Mark it as convex polygon. */

    return PPDisk;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to create a POLYGON/LINE directly from its specified vertices.   M
*   The validity of the elements in the provided list is tested to make sure M
* they are vectors or points.                                                M
*   No test is made to make sure all vertices are on one plane, and that no  M
* two vertices are similar.						     M
*                                                                            *
* PARAMETERS:                                                                M
*   PObjList:     List of vertices/points to construct as a polygon/line.    M
*   IsPolyline:   If TRUE, make a polyline, otherwise a polygon.             M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:   A polygon/line constructed from PObjList.            M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenPOLYGONObject, polygon, polyline, primitives                      M
*****************************************************************************/
IPObjectStruct *PrimGenPOLYGONObject(IPObjectStruct *PObjList, int IsPolyline)
{
    int i,
	NumVertices = 0;
    IPVertexStruct *V, *VHead,
	*VTail = NULL;
    IPPolygonStruct *PPoly;
    IPObjectStruct *PObj, *PObjPoly;

    if (!IP_IS_OLST_OBJ(PObjList))
	IritFatalError("GenPOLYObject: Not object list object!");

    i = 0;
    while ((PObj = ListObjectGet(PObjList, i++)) != NULL) {
	if (!IP_IS_VEC_OBJ(PObj) &&
	    !IP_IS_POINT_OBJ(PObj) &&
	    !IP_IS_CTLPT_OBJ(PObj)) {
	    IritWarningError("None vector object found in list, empty object result.");
	    return NULL;
	}
	NumVertices++;
    }

    if (NumVertices < 2 + !IsPolyline) {
	IritWarningError("Too few vertices, empty object result.");
	return NULL;
    }

    PPoly = IPAllocPolygon(0, 0, VHead = NULL, NULL);
    i = 0;
    while ((PObj = ListObjectGet(PObjList, i++)) != NULL) {
	V = IPAllocVertex(0, 0, NULL, NULL);
	if (IP_IS_VEC_OBJ(PObj))
	    PT_COPY(V -> Coord, PObj -> U.Vec);
	else if (IP_IS_POINT_OBJ(PObj))
	    PT_COPY(V -> Coord, PObj -> U.Pt);
	else if (IP_IS_CTLPT_OBJ(PObj)) {
	    IPObjectStruct
	        *PVec = IritPrsrCoerceObjectTo(PObj, IP_OBJ_VECTOR);

	    PT_COPY(V -> Coord, PVec -> U.Vec);
	    IPFreeObject(PVec);
	}

	if (VHead == NULL) {
	    PPoly -> PVertex = VHead = VTail = V;
	}
	else {
	    VTail -> Pnext = V;
	    VTail = V;
	}
    }

    PObjPoly = GenPolyObject("", NULL, NULL);    /* Gen. POLY object itself: */
    PObjPoly -> U.Pl = PPoly;
    if (IsPolyline) {
	IP_SET_POLYLINE_OBJ(PObjPoly);
    }
    else {
	IP_SET_POLYGON_OBJ(PObjPoly);

	VTail -> Pnext = VHead;		      /* Close the vertex list loop. */

	/* Update polygon plane equation and vertices normals. */
	CGPlaneFrom3Points(PPoly -> Plane, VHead -> Coord,
			   VHead -> Pnext -> Coord,
			   VHead -> Pnext -> Pnext -> Coord);
	V = VHead;
	do {
	    PT_COPY(V -> Normal, PPoly -> Plane);
	    V = V -> Pnext;
	}
	while (V != VHead);
    }

    return PObjPoly;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Routine to create an OBJECT directly from set of sspecified polygons.      m
*   No test is made for the validity of the model in any sense.		     M
*                                                                            *
* PARAMETERS:                                                                M
*   PObjList:      List of polygonal objects.                                M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:   A single object containing all polygons in all       M
*                       provided objects, by a simple union.                 M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenObjectFromPolyList, primitives                                    M
*****************************************************************************/
IPObjectStruct *PrimGenObjectFromPolyList(IPObjectStruct *PObjList)
{
    int i;
    IPPolygonStruct *PPoly,
	*PTail = NULL;
    IPObjectStruct *PObj, *PObjPoly;

    if (!IP_IS_OLST_OBJ(PObjList))
	IritFatalError("GenObjectFromPolyList: Not object list object!");

    i = 0;
    while ((PObj = ListObjectGet(PObjList, i++)) != NULL) {
	if (!IP_IS_POLY_OBJ(PObj)) {
	    IritWarningError("None polygon object found in list, empty object result.");
	    return NULL;
	}
    }

    PObjPoly = GenPolyObject("", NULL, NULL);/* Gen. the POLY object itself: */
    i = 0;
    while ((PObj = ListObjectGet(PObjList, i)) != NULL) {
	if (i++ == 0) {
	    if (IP_IS_POLYLINE_OBJ(PObj))
	        IP_SET_POLYLINE_OBJ(PObjPoly);
	    else
		IP_SET_POLYGON_OBJ(PObjPoly);
	}
	else {
	    if ((IP_IS_POLYLINE_OBJ(PObj) && !IP_IS_POLYLINE_OBJ(PObjPoly)) ||
		(IP_IS_POLYGON_OBJ(PObj) && !IP_IS_POLYGON_OBJ(PObjPoly))) {
		IritWarningError("Polygons mixed with polylines.");
		return NULL;
	    }	      
	}

	PPoly = CopyPolygonList(PObj -> U.Pl);

	if (PTail == NULL) {
	    PObjPoly -> U.Pl = PPoly;
	}
	else {
	    PTail -> Pnext = PPoly;
	}
	for (PTail = PPoly; PTail -> Pnext != NULL; PTail = PTail -> Pnext);
    }

    return PObjPoly;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Not supported.                                                             *
*                                                                            *
* PARAMETERS:                                                                *
*   PObj:                                                                    *
*                                                                            *
* RETURN VALUE:                                                              *
*   IPObjectStruct *:                                                        *
*****************************************************************************/
IPObjectStruct *PrimGenCROSSECObject(IPObjectStruct *PObj)
{
    if (PObj && !IP_IS_POLY_OBJ(PObj))
	IritFatalError("CrossSec: operation on non polygonal object");

    IritWarningError("GenCrossSecObject not implemented");

    return NULL;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to a create surface of revolution by rotating the given cross    M
* section along the Z axis.						     M
*   Input can either be a polygon/line or a freefrom curve object.           M
*   If input is a polyline/gon, it must never be coplanar with the Z axis.   M
*   See also PrimSetResolution on fineness control of approximation of the   M
* primitive using flat faces.                                                M
*                                                                            *
* PARAMETERS:                                                                M
*   Cross:     To rotate around the Z axis forming a surface of revolution.  M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:    A surface of revolution.                            M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenSURFREVObject, surface of revolution, primitives                  M
*****************************************************************************/
IPObjectStruct *PrimGenSURFREVObject(IPObjectStruct *Cross)
{
    int i, j;
    RealType XYSize;
    MatrixType Mat;			   /* Rotation Matrix around Z axes. */
    IPVertexStruct *V, *V1, *V1Head, *V2, *V2Head, *VIn, *VInHead;
    IPPolygonStruct *Pl1, *Pl2, *PlIn, *PlNew = NULL;
    IPObjectStruct *PSurfRev;

    if (!IP_IS_POLY_OBJ(Cross) && !IP_IS_CRV_OBJ(Cross)) {
	IritWarningError("Cross section is not poly/crv. Empty object result");
	return NULL;
    }

    if (IP_IS_POLY_OBJ(Cross)) {
	if (APX_EQ(Cross -> U.Pl -> Plane[0], 0.0) &&
	    APX_EQ(Cross -> U.Pl -> Plane[1], 0.0)) {
	    IritWarningError("Cross-section perpendicular to Z. Empty object result");
	    return NULL;
	}

	Pl1 = IPAllocPolygon(0, 0,
		   V1Head = CopyVertexList(Cross -> U.Pl -> PVertex), NULL);
	PLANE_COPY(Pl1 -> Plane, Cross -> U.Pl -> Plane);
	Pl2 = IPAllocPolygon(0, 0,
		   V2Head = CopyVertexList(Cross -> U.Pl -> PVertex), NULL);
	PLANE_COPY(Pl2 -> Plane, Cross -> U.Pl -> Plane);
	PlIn = GenInsidePoly(Pl1);
	VInHead = PlIn -> PVertex;
	MatGenMatRotZ1(2.0 * M_PI / GlblResolution, Mat);

	for (i = 0; i < GlblResolution; i++)
	{
	    V2 = V2Head;
	    do {
		MatMultVecby4by4(V2 -> Coord, V2 -> Coord , Mat);
		V2 = V2 -> Pnext;
	    }
	    while (V2 != NULL && V2 != V2Head);

	    V1 = V1Head;
	    if (i < GlblResolution - 1) /* If this is last loop use original */
	        V2 = V2Head; /* poly as we might accumulate error during the */
	    else			/* transformations along the circle. */
		V2 = Cross -> U.Pl -> PVertex;
	    VIn = VInHead;

	    do {
		PlNew = GenPolygon4Vrtx(V1 -> Coord, V1 -> Pnext -> Coord,
				        V2 -> Pnext -> Coord, V2 -> Coord,
				        VIn -> Coord, PlNew);

	        /* Update normals: */
	        for (j = 0, V = PlNew -> PVertex; j < 4; j++, V = V -> Pnext) {
		    V -> Normal[0] = V -> Coord[0];
		    V -> Normal[1] = V -> Coord[1];
		    V -> Normal[2] = 0.0;
		    /* Make sure normal does not point in opposite direction.*/
		    if (DOT_PROD(V -> Normal, PlNew -> Plane) < 0.0)
		        PT_SCALE(V -> Normal, -1.0);

		    /* Since Z normal component should be fixed for all normals: */
		    V -> Normal[2] = PlNew -> Plane[2];
		    XYSize = APX_EQ(ABS(PlNew -> Plane[2]), 1.0) ?
					0.0 : 1 - SQR(PlNew -> Plane[2]);
		    XYSize = sqrt(XYSize / (SQR(V -> Coord[0]) +
					    SQR(V -> Coord[1])));
		    V -> Normal[0] *= XYSize;
		    V -> Normal[1] *= XYSize;
	        }

	        VIn = VIn -> Pnext;
	        V1 = V1 -> Pnext;
	        V2 = V2 -> Pnext;
	    }
	    while (V1 -> Pnext != NULL && V1 != V1Head);

	    V1 = V1Head;
	    do {
	        MatMultVecby4by4(V1 -> Coord, V1 -> Coord , Mat);
	        V1 = V1 -> Pnext;
	    }
	    while (V1 != NULL && V1 != V1Head);
	    VIn = VInHead;
	    do {
	        MatMultVecby4by4(VIn -> Coord, VIn -> Coord , Mat);
	        VIn = VIn -> Pnext;
	    }
	    while (VIn != NULL && VIn != VInHead);
        }

        IPFreePolygonList(PlIn);
        IPFreePolygonList(Pl1);
        IPFreePolygonList(Pl2);

        PSurfRev = GenPolyObject("", NULL, NULL);
	PSurfRev -> U.Pl = PlNew;

        return PSurfRev;
    }
    else if (IP_IS_CRV_OBJ(Cross)) {
	if (CAGD_NUM_OF_PT_COORD(Cross -> U.Crvs -> PType) < 3) {
	    IritWarningError("Cross-section perpendicular to Z. Empty object result");
	    return NULL;
	}

        PSurfRev = GenSRFObject(CagdSurfaceRev(Cross -> U.Crvs));
	return PSurfRev;
    }
    else
        return NULL;
}

/*****************************************************************************
* DESCRIPTION:                                                               M
*   Routine to a create an extrusion surface out of the given cross section  M
* and the given direction.                                                   M
*   Input can either be a polygon/line or a freefrom curve object.           M
*   If input is a polyline/gon, it must never be coplanar with Dir.          M
*   See also PrimSetResolution on fineness control of approximation of the   M
* primitive using flat faces.                                                M
*									     M
* PARAMETERS:                                                                M
*   Cross:     To extrude in direction Dir.                                  M
*   Dir:       Direction and magnitude of extrusion.                         M
*                                                                            *
* RETURN VALUE:                                                              M
*   IPObjectStruct *:     An extrusion surface.                              M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimGenEXTRUDEObject, extrusion surface, primitives                      M
*****************************************************************************/
IPObjectStruct *PrimGenEXTRUDEObject(IPObjectStruct *Cross, VectorType Dir)
{
    int i;
    RealType R;
    CagdVecStruct CagdDir;
    IPVertexStruct *V1, *V1Head, *V2, *VIn;
    IPPolygonStruct *PBase1, *PBase2, *Pl, *PlIn;
    IPObjectStruct *PExtrude;

    if (!IP_IS_POLY_OBJ(Cross) && !IP_IS_CRV_OBJ(Cross)) {
	IritWarningError("Cross section is not poly/crv. Empty object result");
	return NULL;
    }

    if (IP_IS_POLY_OBJ(Cross)) {
	R = DOT_PROD(Cross -> U.Pl -> Plane, Dir);
	if (APX_EQ(R, 0.0)) {
	    IritWarningError("Extrusion direction in cross-section plane. Empty object result");
	    return NULL;
	}

	/* Prepare two bases (update their plane normal to point INDISE): */
	PBase1 = IPAllocPolygon(0, 0,
			      CopyVertexList(Cross -> U.Pl -> PVertex), NULL);
	Pl = PBase2 = IPAllocPolygon(0, 0,
			     CopyVertexList(Cross -> U.Pl -> PVertex), PBase1);
	V1 = V1Head = PBase2 -> PVertex;
	do {
	    PT_ADD(V1 -> Coord, Dir, V1 -> Coord);
	    V1 = V1 -> Pnext;
	}
	while (V1 != NULL && V1 != V1Head);
	if (R > 0.0) {
	    PLANE_COPY(PBase1 -> Plane, Cross -> U.Pl -> Plane);
	    for (i = 0; i < 3; i++)
		PBase2 -> Plane[i] = (-Cross -> U.Pl -> Plane[i]);
	    PBase2 -> Plane[3] = (-DOT_PROD(PBase2 -> Plane,
					    PBase2 -> PVertex -> Coord));
	}
	else {
	    for (i = 0; i < 4; i++)
		PBase1 -> Plane[i] = (-Cross -> U.Pl -> Plane[i]);
	    PLANE_COPY(PBase2 -> Plane, Cross -> U.Pl -> Plane);
	    PBase2 -> Plane[3] = (-DOT_PROD(PBase2 -> Plane,
					    PBase2 -> PVertex -> Coord));
	}

	/* Now generate all the 4 corner polygon between the two bases: */
	V1 = V1Head = PBase1 -> PVertex;
	V2 = PBase2 -> PVertex;
	PlIn = GenInsidePoly(PBase1);
	VIn = PlIn -> PVertex;
	do {
	    Pl = GenPolygon4Vrtx(V1 -> Coord, V1 -> Pnext -> Coord,
		         V2 -> Pnext -> Coord, V2 -> Coord, VIn -> Coord, Pl);
	    VIn = VIn -> Pnext;
	    V1 = V1 -> Pnext;
	    V2 = V2 -> Pnext;
	}
	while (V1 -> Pnext != NULL && V1 != V1Head);

	IPFreePolygonList(PlIn);

	PExtrude = GenPolyObject("", NULL, NULL);
	PExtrude -> U.Pl = Pl;

	/* Update all the polygon vertices normals. */
	for (Pl = PExtrude -> U.Pl; Pl != NULL; Pl = Pl -> Pnext) {
	    V1 = V1Head = Pl -> PVertex;
    	    do {
    	        PT_COPY(V1 -> Normal, Pl -> Plane);
    	        V1 = V1 -> Pnext;
    	    }
    	    while (V1 != NULL && V1 != V1Head);
	}

	return PExtrude;
    }
    else if (IP_IS_CRV_OBJ(Cross)) {
        for (i = 0; i < 3; i++)
	    CagdDir.Vec[i] = Dir[i];
        PExtrude = GenSRFObject(CagdExtrudeSrf(Cross -> U.Crvs, &CagdDir));

        return PExtrude;
    }
    else
        return NULL;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Routine to create a pseudo polygon out of a given polygon such that each *
* vertex Vi is in the inside side of the corresponding edge ViVi+1 in the    *
* given polygon. Used in polygon generation for EXTRUDE/SURFREV operations.  *
*                                                                            *
* PARAMETERS:                                                                *
*   Pl:        Input polygon                                                 *
*                                                                            *
* RETURN VALUE:                                                              *
*   IPPolygonStruct *:  Pseudo one.                                          *
*****************************************************************************/
static IPPolygonStruct *GenInsidePoly(IPPolygonStruct *Pl)
{
    int Axes;
    RealType Dx, Dy;
    PointType Pt;
    MatrixType Mat;
    IPPolygonStruct *PlIn;
    IPVertexStruct *VHead, *V, *Vnext, *VInHead,
	*VIn = NULL;

    PlIn = IPAllocPolygon(0, 0, VInHead = NULL, NULL);

    /* Generate transformation matrix to bring polygon to a XY parallel      */
    /* plane, and transform a copy of the polygon to that plane.	     */
    GenRotateMatrix(Mat, Pl -> Plane);
    /* We dont want to modify original! */
    VHead = V = CopyVertexList(Pl -> PVertex);
    Pl = IPAllocPolygon(0, 0, VHead, NULL);
    do {
	MatMultVecby4by4(V -> Coord, V -> Coord, Mat);
	V = V -> Pnext;
    }
    while (V != NULL && V != VHead);

    V = VHead;
    do {
	Vnext = V -> Pnext;
	Dx = ABS(V -> Coord[0] - Vnext -> Coord[0]);
	Dy = ABS(V -> Coord[1] - Vnext -> Coord[1]);
	/* Prepare middle point. */
	Pt[0] = (V -> Coord[0] + Vnext -> Coord[0]) / 2.0;
	Pt[1] = (V -> Coord[1] + Vnext -> Coord[1]) / 2.0;
	Pt[2] = V -> Coord[2];
	/* If Dx > Dy fire ray in +Y direction, otherwise in +X direction    */
	/* and if number of intersections is even (excluding the given point */
	/* itself) then that direction is the outside, otherwise, its inside.*/
	Axes = (Dx > Dy ? 1 : 0);
	if (CGPolygonRayInter(Pl, Pt, Axes) % 2 == 0) {
	    /* The amount we move along Axes is not of a big meaning as long */
	    /* as it is not zero, so MAX(Dx, Dy) guarantee non zero value... */
	    Pt[Axes] -= MAX(Dx, Dy);
	}
	else {
	    Pt[Axes] += MAX(Dx, Dy);
	}

	/* Now Pt holds point which is in the inside part of vertex V, Vnext.*/
	/* Put it in the pseudo inside polygon PlIn:			     */
	if (VInHead) {
	    VIn -> Pnext = IPAllocVertex(0, 0, NULL, NULL);
	    VIn = VIn -> Pnext;
	}
	else {
	    PlIn -> PVertex = VInHead = VIn = IPAllocVertex(0, 0, NULL, NULL);
	}
	PT_COPY(VIn -> Coord, Pt);

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

    IPFreePolygonList(Pl);	      /* Free copied (and trans.) vrtx list. */

    /* Transform PlIn to the plane where original Pl is... */
    if (!MatInverseMatrix(Mat, Mat))		 /* Find the inverse matrix. */
	IritFatalError("GenInsidePoly: Inverse matrix does not exits");
    VIn = VInHead;
    do {
	MatMultVecby4by4(VIn -> Coord, VIn -> Coord, Mat);
	VIn = VIn -> Pnext;
    }
    while (VIn != NULL && VIn != VInHead);

    return PlIn;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Routine to create a polygon out of a list of 4 vertices V1/2/3/4.	     *
*   The fifth vertex is inside (actually, this is not true, as this point    *
* will be in the positive part of the plane, which only locally in the       *
* object...) the object, so the polygon's normal direction can be evaluated  *
* uniquely.								     *
*   No test is made to make sure the 4 points are co-planar...		     *
*   The points are placed in order.                      		     *
*                                                                            *
* PARAMETERS:                                                                *
*   V1, V2, V3, V4:    Four vertices of the constructed polygon.             *
*   Vin:               A vertex that can be assumed to be inside the         *
*                      object for normal evaluation of the plane of polygon. *
*   Pnext:             Next is chain of polygons, in linked list.            *
*                                                                            *
* RETURN VALUE:                                                              *
*   IPPolygonStruct *: The constructed polygon.                              *
*****************************************************************************/
static IPPolygonStruct *GenPolygon4Vrtx(VectorType V1,
					VectorType V2,
					VectorType V3,
					VectorType V4,
					VectorType Vin,
					IPPolygonStruct *Pnext)
{
    IPPolygonStruct *PPoly;
    IPVertexStruct *V;

    PPoly = IPAllocPolygon(0, 0, V = IPAllocVertex(0, 0, NULL, NULL), Pnext);
    PT_COPY(V -> Coord, V1);

    V -> Pnext = IPAllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
    PT_COPY(V -> Coord, V2);

    V -> Pnext = IPAllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
    PT_COPY(V -> Coord, V3);

    V -> Pnext = IPAllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
    PT_COPY(V -> Coord, V4);

    V -> Pnext = PPoly -> PVertex;	   /* Make the Vertex list circular. */

    IritPrsrUpdatePolyPlane2(PPoly, Vin);	   /* Update plane equation. */

    IP_SET_CONVEX_POLY(PPoly);		       /* Mark it as convex polygon. */

    return PPoly;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
* Routine to create a polygon out of a list of 3vertices V1/2/3.  	     *
*   The fifth vertex is inside (actually, this is not true, as this point    *
* will be in the positive part of the plane, which only locally in the       *
* object...) the object, so the polygon's normal direction can be evaluated  *
* uniquely.								     *
*   No test is made to make sure the 4 points are co-planar...		     *
*   The points are placed in order.					     *
*                                                                            *
* PARAMETERS:                                                                *
*   V1, V2, V3:    Three vertices of the constructed polygon.                *
*   Vin:           A vertex that can be assumed to be inside the             *
*                  object for normal evaluation of the plane of polygon.     *
*   Pnext:         Next is chain of polygons, in linked list.                *
*                                                                            *
* RETURN VALUE:                                                              *
*   IPPolygonStruct *: The constructed polygon.                              *
*****************************************************************************/
static IPPolygonStruct *GenPolygon3Vrtx(VectorType V1,
					VectorType V2,
					VectorType V3,
					VectorType Vin,
					IPPolygonStruct *Pnext)
{
    IPPolygonStruct *PPoly;
    IPVertexStruct *V;

    PPoly = IPAllocPolygon(0, 0, V = IPAllocVertex(0, 0, NULL, NULL), Pnext);
    PT_COPY(V -> Coord, V1);

    V -> Pnext = IPAllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
    PT_COPY(V -> Coord, V2);

    V -> Pnext = IPAllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
    PT_COPY(V -> Coord, V3);

    V -> Pnext = PPoly -> PVertex;	   /* Make the Vertex list circular. */

    IritPrsrUpdatePolyPlane2(PPoly, Vin);	   /* Update plane equation. */

    IP_SET_CONVEX_POLY(PPoly);		       /* Mark it as convex polygon. */

    return PPoly;
}

/*****************************************************************************
* DESCRIPTION:                                                               *
*   Routine to update the Vertex Pt normal equation. The normal should point *
* InPt but should be perpendicular to the line Pt-PerpPt if Perpendicular is *
* TRUE. THe normal is normalized to a unit length.			     *
*                                                                            *
* PARAMETERS:                                                                *
*   Normal:        To update.                                                *
*   Pt:            The normal belongs to this location.                      *
*   InPt:          To define the normal direction as InPt - Pt.              *
*   Perpendicular: If TRUE, Normal also perpedicular to PerpPt - Pt.         *
*   PerpPt:        To define perpendicular relation as PerpPt - Pt.          *
*                                                                            *
* RETURN VALUE:                                                              *
*   void                                                                     *
*****************************************************************************/
static void UpdateVertexNormal(NormalType Normal,
			       PointType Pt,
			       PointType InPt,
			       int Perpendicular,
			       PointType PerpPt)
{
    VectorType V1, V2;

    PT_SUB(Normal, InPt, Pt);

    if (Perpendicular) {
	PT_SUB(V1, PerpPt, Pt);
	GMVecCrossProd(V2, V1, Normal);
	GMVecCrossProd(Normal, V2, V1);
    }

    PT_NORMALIZE(Normal);
}

/*****************************************************************************
* DESCRIPTION:                                                               M
* Routine to set the polygonal resolution (fineness).                        M
*   Resolution sroutghly the number of edges a circular primitive will have  M
* along the entire circle.                                                   M
*                                                                            *
* PARAMETERS:                                                                M
*   Resolution:    To set as new resolution for all primitve constructors.   M
*                                                                            *
* RETURN VALUE:                                                              M
*   void                                                                     M
*                                                                            *
* KEYWORDS:                                                                  M
*   PrimSetResolution, primitives, resolution                                M
*****************************************************************************/
void PrimSetResolution(int Resolution)
{
    GlblResolution = MAX(Resolution, MIN_RESOLUTION);
}
