/*****************************************************************************
* "Irit" - the 3d 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 "program.h"
#include "allocate.h"
#include "attribut.h"
#include "convex.h"
#include "geomat3d.h"
#include "graphgen.h"
#include "objects.h"
#include "primitiv.h"
#include "windows.h"
#define MIN_RESOLUTION 4
static PolygonStruct *GenInsidePoly(PolygonStruct *Pl);
static PolygonStruct *GenPolygon4Vrtx(VectorType V1, VectorType V2,
VectorType V3, VectorType V4, VectorType Vin, PolygonStruct *Pnext);
static PolygonStruct *GenPolygon3Vrtx(VectorType V1, VectorType V2,
VectorType V3, VectorType Vin, PolygonStruct *Pnext);
static void GenTransformMatrix(MatrixType Mat, VectorType Trans,
VectorType Dir, RealType Scale);
static void UpdateVertexNormal(NormalType Normal, PointType Pt, PointType InPt,
int Perpendicular, PointType PerpPt);
/*****************************************************************************
* Routine to create a BOX geometric object defined by Pt - the minimun *
* 3d point, and Width - Dx Dy & Dz vector. 4 *
* Order of vertices is as 5 7 *
* follows in the picture: | 6 | *
* | | | *
* (Note vertex 0 is hidden behind edge 2-6) | | | *
* 1 | 3 *
* 2 *
*****************************************************************************/
ObjectStruct * GenBOXObject(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 GenGBOXObject(Pt, Dir1, Dir2, Dir3);
}
/*****************************************************************************
* Routine to create a GBOX geometric object defined by Pt - the minimun *
* 3d point, and 3 direction Vectors Dir1, Dir2, Dir3. If two of the *
* direction vectors are parallel the GBOX converges to zero volume. A NULL *
* pointer is returned in that case! *
* 4 *
* Order of vertices is as 5 7 *
* follows in the picture: | 6 | *
* | | | *
* (Note vertex 0 is hidden behind edge 2-6) | | | *
* 1 | 3 *
* 2 *
*****************************************************************************/
ObjectStruct * GenGBOXObject(VectorType Pt,
VectorType Dir1, VectorType Dir2, VectorType Dir3)
{
int i;
VectorType Temp;
VectorType V[8]; /* Hold 8 vertices of BOX. */
VertexStruct *PVertex;
PolygonStruct *PPolygon;
ObjectStruct *PBox;
VecCrossProd(Temp, Dir1, Dir2);
if (APX_EQ(PT_LENGTH(Temp), 0.0)) return NULL;
VecCrossProd(Temp, Dir2, Dir3);
if (APX_EQ(PT_LENGTH(Temp), 0.0)) return NULL;
VecCrossProd(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.P = GenPolygon4Vrtx(V[0], V[1], V[3], V[2], V[4], PBox -> U.Pl.P);
PBox -> U.Pl.P = GenPolygon4Vrtx(V[6], V[7], V[5], V[4], V[0], PBox -> U.Pl.P);
PBox -> U.Pl.P = GenPolygon4Vrtx(V[4], V[5], V[1], V[0], V[2], PBox -> U.Pl.P);
PBox -> U.Pl.P = GenPolygon4Vrtx(V[5], V[7], V[3], V[1], V[0], PBox -> U.Pl.P);
PBox -> U.Pl.P = GenPolygon4Vrtx(V[7], V[6], V[2], V[3], V[1], PBox -> U.Pl.P);
PBox -> U.Pl.P = GenPolygon4Vrtx(V[6], V[4], V[0], V[2], V[3], PBox -> U.Pl.P);
/* Update the vertices normals using the polygon plane equation: */
for (PPolygon = PBox -> U.Pl.P;
PPolygon != NULL;
PPolygon = PPolygon -> Pnext) {
PVertex = PPolygon -> V;
do {
PT_COPY(PVertex -> Normal, PPolygon -> Plane);
PVertex = PVertex -> Pnext;
}
while (PVertex != PPolygon -> V);
}
SetObjectColor(PBox, GlblPrimColor); /* Set its default color. */
return PBox;
}
/*****************************************************************************
* Routine to fetch the resolution parameter from the RESOLUTION object. *
* If ClipToMin TRUE, the Resolution is clipped to not be below MIN_RES. *
*****************************************************************************/
int GetResolution(int ClipToMin)
{
int Resolution;
ObjectStruct *PObj = GetObject("RESOLUTION");
if (PObj == NULL || !IS_NUM_OBJ(PObj)) {
WndwInputWindowPutStr("No numeric object name RESOLUTION is defined");
Resolution = DEFAULT_RESOLUTION;
}
else
Resolution = ClipToMin ? MAX(((int) (PObj -> U.R)), MIN_RESOLUTION)
: (int) (PObj -> U.R);
Resolution = (Resolution / 2) * 2; /* Make sure its an even number. */
return Resolution;
}
/*****************************************************************************
* Routine to create a CONE geometric object defined by Pt - the base *
* 3d center point, Dir - the cone direction and length, and base radius R. *
*****************************************************************************/
ObjectStruct *GenCONEObject(VectorType Pt, VectorType Dir, RealType *R)
{
int i, Resolution;
RealType Angle, AngleStep;
PointType LastCirclePt, CirclePt, ApexPt;
NormalType LastCircleNrml, CircleNrml, ApexNrml;
MatrixType Mat;
VertexStruct *VBase, *PVertex;
PolygonStruct *PBase;
ObjectStruct *PCone;
Resolution = GetResolution(TRUE); /* Get refinement factor of object. */
GenTransformMatrix(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 = AllocPolygon(0, 0, VBase = AllocVertex(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 -> Pt, LastCirclePt); /* Update first pt in base polygon. */
PT_COPY(VBase -> Normal, Dir);
AngleStep = M_PI * 2 / Resolution;
for (i = 1; i <= Resolution; 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.P = GenPolygon3Vrtx(LastCirclePt, ApexPt,
CirclePt, Pt, PCone -> U.Pl.P);
/* Update the normals for this cone side polygon vertices: */
PVertex = PCone -> U.Pl.P -> V;
PT_COPY(PVertex -> Normal, LastCircleNrml);
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);
PVertex = PVertex -> Pnext;
PT_COPY(PVertex -> Normal, CircleNrml);
/* And add this vertex to base polygon: */
if (i == Resolution) /* Its last point - make it circular. */
VBase -> Pnext = PBase -> V;
else {
VBase -> Pnext = AllocVertex(0, 0, NULL, NULL);
VBase = VBase -> Pnext;
PT_COPY(VBase -> Normal, Dir);
PT_COPY(VBase -> Pt, CirclePt);
}
PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
PT_COPY(LastCircleNrml, CircleNrml);
}
UpdatePolyPlane(PBase, ApexPt); /* Update base polygon plane equation. */
SET_CONVEX_POLY(PBase); /* Mark it as convex polygon. */
PBase -> Pnext = PCone -> U.Pl.P;/* And stick it into the cone polygons. */
PCone -> U.Pl.P = PBase;
SetObjectColor(PCone, GlblPrimColor); /* Set its default color. */
return PCone;
}
/*****************************************************************************
* Routine to create a CONE2 geometric object defined by Pt - main base *
* 3d center point, Dir - the cone direction and length, and bases R1 and R2. *
*****************************************************************************/
ObjectStruct *GenCONE2Object(VectorType Pt, VectorType Dir, RealType *R1,
RealType *R2)
{
int i, Resolution;
RealType Angle, AngleStep;
PointType LastCirclePt, CirclePt, ApexPt, LastApexPt1, ApexPt1;
NormalType LastCircleNrml, CircleNrml;
VectorType InvDir;
MatrixType Mat1, Mat2;
VertexStruct *VBase1, *VBase2, *PVertex;
PolygonStruct *PBase1, *PBase2;
ObjectStruct *PCone;
Resolution = GetResolution(TRUE); /* Get refinement factor of object. */
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);
GenTransformMatrix(Mat1, Pt, Dir, *R1); /* Transform from unit circle. */
GenTransformMatrix(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 = AllocPolygon(0, 0, VBase1 = AllocVertex(0, 0, NULL, NULL), NULL);
PBase2 = AllocPolygon(0, 0, VBase2 = AllocVertex(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 -> Pt, LastCirclePt);/* Update first pt in base1 polygon. */
PT_COPY(VBase1 -> Normal, Dir);
PT_COPY(VBase2 -> Pt, LastApexPt1); /* Update first pt in base2 polygon. */
PT_COPY(VBase2 -> Normal, InvDir);
AngleStep = M_PI * 2 / Resolution;
for (i = 1; i <= Resolution; 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.P = GenPolygon4Vrtx(LastCirclePt, LastApexPt1, ApexPt1,
CirclePt, Pt, PCone -> U.Pl.P);
/* Update the normals for this cone side polygon vertices: */
PVertex = PCone -> U.Pl.P -> V;
PT_COPY(PVertex -> Normal, LastCircleNrml);
PVertex = PVertex -> Pnext;
PT_COPY(PVertex -> Normal, LastCircleNrml );
PVertex = PVertex -> Pnext;
PT_COPY(PVertex -> Normal, CircleNrml);
PVertex = PVertex -> Pnext;
PT_COPY(PVertex -> Normal, CircleNrml);
/* And add these vertices to base polygons: */
if (i == Resolution) { /* Its last point - make it circular. */
VBase1 -> Pnext = PBase1 -> V;
VBase2 -> Pnext = PBase2 -> V;
}
else {
VBase1 -> Pnext = AllocVertex(0, 0, NULL, NULL);
VBase1 = VBase1 -> Pnext;
PT_COPY(VBase1 -> Pt, CirclePt);
PT_COPY(VBase1 -> Normal, Dir);
VBase2 -> Pnext = AllocVertex(0, 0, NULL, NULL);
VBase2 = VBase2 -> Pnext;
PT_COPY(VBase2 -> Pt, 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);
}
UpdatePolyPlane(PBase1, ApexPt); /* Update base polygon plane equation. */
SET_CONVEX_POLY(PBase1); /* Mark it as convex polygon. */
UpdatePolyPlane(PBase2, Pt); /* Update base polygon plane equation. */
SET_CONVEX_POLY(PBase2); /* Mark it as convex polygon. */
PBase1 -> Pnext = PCone -> U.Pl.P; /* And stick into the cone polygons. */
PCone -> U.Pl.P = PBase1;
PBase2 -> Pnext = PCone -> U.Pl.P;
PCone -> U.Pl.P = PBase2;
SetObjectColor(PCone, GlblPrimColor); /* Set its default color. */
return PCone;
}
/*****************************************************************************
* Routine to create a CYLINder geometric object defined by Pt - the base *
* 3d center point, Dir - the cylin direction and length, and base radius R. *
* The second base is defined from first one by translating it by vector *
* Dir, and its points are prefixed with T. *
*****************************************************************************/
ObjectStruct *GenCYLINObject(VectorType Pt, VectorType Dir, RealType *R)
{
int i, Resolution;
RealType Angle, AngleStep;
PointType LastCirclePt, CirclePt, TLastCirclePt, TCirclePt, TPt, Dummy;
VectorType ForwardDir, BackwardDir;
NormalType LastCircleNrml, CircleNrml;
MatrixType Mat;
VertexStruct *VBase1, *VBase2, *PVertex;
PolygonStruct *PBase1, *PBase2;
ObjectStruct *PCylin;
Resolution = GetResolution(TRUE); /* Get refinement factor of object. */
GenTransformMatrix(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 = AllocPolygon(0, 0, VBase1 = AllocVertex(0, 0, NULL, NULL), NULL);
PBase2 = AllocPolygon(0, 0, VBase2 = AllocVertex(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 -> Pt, LastCirclePt);/* Update first pt in base1 polygon. */
PT_COPY(VBase1 -> Normal, ForwardDir);
PT_ADD(TLastCirclePt, LastCirclePt, Dir); /* Translated circle (by Dir). */
PT_COPY(VBase2 -> Pt, TLastCirclePt);/* Update first pt in base2 polygon.*/
PT_COPY(VBase2 -> Normal, BackwardDir);
AngleStep = M_PI * 2 / Resolution;
for (i = 1; i <= Resolution; 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.P = GenPolygon4Vrtx(TLastCirclePt, TCirclePt, CirclePt,
LastCirclePt, Pt, PCylin -> U.Pl.P);
/* Update the normals for this cylinder side polygon vertices: */
PVertex = PCylin -> U.Pl.P -> V;
PT_COPY(PVertex -> Normal, LastCircleNrml);
PVertex = PVertex -> Pnext;
PT_COPY(PVertex -> Normal, CircleNrml);
PVertex = PVertex -> Pnext;
PT_COPY(PVertex -> Normal, CircleNrml);
PVertex = PVertex -> Pnext;
PT_COPY(PVertex -> Normal, LastCircleNrml);
/* 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 == Resolution) { /* Its last point - make it circular. */
VBase1 -> Pnext = PBase1 -> V;
VBase2 -> Pnext = PBase2 -> V;
}
else {
VBase1 -> Pnext = AllocVertex(0, 0, NULL, NULL);
VBase1 = VBase1 -> Pnext;
PT_COPY(VBase1 -> Pt, CirclePt);
PT_COPY(VBase1 -> Normal, ForwardDir);
PBase2 -> V = AllocVertex(0, 0, NULL, PBase2 -> V);
PT_COPY(PBase2 -> V -> Pt, TCirclePt);
PT_COPY(PBase2 -> V -> Normal, BackwardDir);
}
PT_COPY(LastCirclePt, CirclePt);/* Save pt in last pt for next time. */
PT_COPY(TLastCirclePt, TCirclePt);
PT_COPY(LastCircleNrml, CircleNrml);
}
UpdatePolyPlane(PBase1, TPt); /* Update base polygon plane equation. */
SET_CONVEX_POLY(PBase1); /* Mark it as convex polygon. */
PBase1 -> Pnext = PCylin -> U.Pl.P; /* And stick it into cylin polygons. */
PCylin -> U.Pl.P = PBase1;
UpdatePolyPlane(PBase2, Pt); /* Update base polygon plane equation. */
SET_CONVEX_POLY(PBase2); /* Mark it as convex polygon. */
PBase2 -> Pnext = PCylin -> U.Pl.P; /* And stick it into cylin polygons. */
PCylin -> U.Pl.P = PBase2;
SetObjectColor(PCylin, GlblPrimColor); /* Set its default color. */
return PCylin;
}
/*****************************************************************************
* Routine to create a SPHERE geometric object defined by Center - sphere *
* 3d center point, and R its radius. *
* Note polygons on the poles are triangles, while others are rectangles *
* Teta is horizontal circle angle, Fee is the vertical (spherical coords.) *
* The vertical axes here is assumed to be the Z axes. *
*****************************************************************************/
ObjectStruct *GenSPHEREObject(VectorType Center, RealType *R)
{
int i, j, k, Resolution;
RealType TetaAngle, TetaAngleStep, FeeAngle, FeeAngleStep,
CosFeeAngle1, SinFeeAngle1, CosFeeAngle2, SinFeeAngle2;
PointType LastCircleLastPt, LastCirclePt, CirclePt, CircleLastPt, Dummy;
VertexStruct *PVertex;
ObjectStruct *PSphere;
Resolution = GetResolution(TRUE); /* Get refinement factor of object. */
PSphere = GenPolyObject("", NULL, NULL);/* Gen the SPHERE object itself: */
TetaAngleStep = M_PI * 2.0 / Resolution; /* Runs from 0 to 2*PI. */
FeeAngleStep = M_PI * 2.0 / Resolution; /* 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 <= Resolution; 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.P = GenPolygon3Vrtx(LastCirclePt, CircleLastPt,
CirclePt, Center, PSphere -> U.Pl.P);
/* Update normals: */
for (j = 0, PVertex = PSphere -> U.Pl.P -> V;
j < 3;
j++, PVertex = PVertex -> Pnext)
UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, Center,
FALSE, Dummy);
PT_COPY(CircleLastPt, CirclePt);/* Save pt in last pt for next time. */
}
/* Generate the middle rectangular polygons: */
for (i = 1; i < Resolution/2-1; i++) { /* For all horizontal circles do. */
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 <= Resolution; j++) {/* Pass whole (horizontal) 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.P = GenPolygon4Vrtx(LastCirclePt, LastCircleLastPt,
CircleLastPt, CirclePt, Center, PSphere -> U.Pl.P);
/* Update normals: */
for (k = 0, PVertex = PSphere -> U.Pl.P -> V;
k < 4;
k++, PVertex = PVertex -> Pnext)
UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, Center,
FALSE, Dummy);
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 <= Resolution; 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.P = GenPolygon3Vrtx(LastCirclePt, CirclePt, CircleLastPt,
Center, PSphere -> U.Pl.P);
/* Update normals: */
for (j = 0, PVertex = PSphere -> U.Pl.P -> V;
j < 3;
j++, PVertex = PVertex -> Pnext)
UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, Center,
FALSE, Dummy);
PT_COPY(CircleLastPt, CirclePt);/* Save pt in last pt for next time. */
}
SetObjectColor(PSphere, GlblPrimColor); /* Set its default color. */
return PSphere;
}
/*****************************************************************************
* Routine to create a TORUS geometric object defined by Center - torus 3d *
* center point, the main torus plane normal Normal, major radius Rmajor and *
* minor radius Rminor (Tube radius). *
* Teta runs on the major circle, Fee on the minor one (Before Mat trans.): *
* X = (Rmajor + Rminor * cos(Fee)) * cos(Teta) *
* Y = (Rmajor + Rminor * cos(Fee)) * sin(Teta) *
* Z = Rminor * sin(Fee) *
*****************************************************************************/
ObjectStruct *GenTORUSObject(VectorType Center, VectorType Normal,
RealType *Rmajor, RealType *Rminor)
{
int i, j, Resolution;
RealType TetaAngle, TetaAngleStep, FeeAngle, FeeAngleStep,
CosFeeAngle1, SinFeeAngle1, CosFeeAngle2, SinFeeAngle2;
PointType LastCircleLastPt, LastCirclePt, CirclePt, CircleLastPt,
LastInPt, InPt, Dummy;
MatrixType Mat;
VertexStruct *PVertex;
ObjectStruct *PTorus;
Resolution = GetResolution(TRUE); /* Get refinement factor of object. */
GenTransformMatrix(Mat, Center, Normal, 1.0); /* Trans from unit circle. */
PTorus = GenPolyObject("", NULL, NULL); /* Gen. the Torus object itself: */
TetaAngleStep = M_PI * 2.0 / Resolution; /* Runs from 0 to 2*PI. */
FeeAngleStep = M_PI * 2.0 / Resolution; /* Runs from 0 to 2*PI. */
for (i = 1; i <= Resolution; 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 <= Resolution; 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.P = GenPolygon4Vrtx(CirclePt, CircleLastPt,
LastCircleLastPt, LastCirclePt, InPt, PTorus -> U.Pl.P);
/* Update normals: */
PVertex = PTorus -> U.Pl.P -> V;
UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, InPt,
FALSE, Dummy);
PVertex = PVertex -> Pnext;
UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, InPt,
FALSE, Dummy);
PVertex = PVertex -> Pnext;
UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, LastInPt,
FALSE, Dummy);
PVertex = PVertex -> Pnext;
UpdateVertexNormal(PVertex -> Normal, PVertex -> Pt, LastInPt,
FALSE, Dummy);
PT_COPY(LastCirclePt, CirclePt); /* Save pt in last pt. */
PT_COPY(LastCircleLastPt, CircleLastPt);
PT_COPY(LastInPt, InPt);
}
}
SetObjectColor(PTorus, GlblPrimColor); /* Set its default color. */
return PTorus;
}
/*****************************************************************************
* Routine to create a PLANE geometric object defined by the normal N and *
* the translation vector T. The object is a FINITE plane (a circle of *
* RESOLUTION point in it...) and its radius is equal to R. *
* The normal direction is assumed to point to the inside of the object. *
*****************************************************************************/
ObjectStruct *GenPLANEObject(VectorType N, VectorType T, RealType *R)
{
int i, Resolution;
RealType Angle, AngleStep;
PointType CirclePt;
MatrixType Mat;
VertexStruct *V;
PolygonStruct *PCirc;
ObjectStruct *PPlane;
Resolution = GetResolution(TRUE); /* Get refinement factor of object. */
GenTransformMatrix(Mat, T, N, *R); /* Transform from unit circle. */
PT_NORMALIZE(N);
PPlane = GenPolyObject("", NULL, NULL); /* Gen. the PLANE object itself: */
PCirc = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), NULL);
PPlane -> U.Pl.P = 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 -> Pt, CirclePt); /* Update first pt in circle polygon. */
PT_COPY(V -> Normal, N);
AngleStep = M_PI * 2 / Resolution;
for (i = 1; i <= Resolution; 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 == Resolution) { /* Its last point - make it circular. */
V -> Pnext = PCirc -> V;
}
else {
V -> Pnext = AllocVertex(0, 0, NULL, NULL);
V = V -> Pnext;
PT_COPY(V -> Pt, CirclePt);
PT_COPY(V -> Normal, N);
}
}
PT_ADD(CirclePt, CirclePt, N); /* Make a point "IN" the circle object. */
UpdatePolyPlane(PCirc, CirclePt); /* Update base polygon plane equation. */
SET_CONVEX_POLY(PCirc); /* Mark it as convex polygon. */
SetObjectColor(PPlane, GlblPrimColor); /* Set its default color. */
return PPlane;
}
/*****************************************************************************
* Routine to create a POLYGON directly from its specified vertices. *
* The validity of list elements is tested to make sure they are vectors as *
* nobody else checked it till now (lists can hold any object type!). *
* No test is made to make sure all vertices are on one plane, and that no *
* two vertices are similar. *
*****************************************************************************/
ObjectStruct *GenPOLYGONObject(ObjectStruct *PObjList)
{
int i, NumVertices = 0;
VertexStruct *V, *VHead, *VTail;
PolygonStruct *PPoly;
ObjectStruct *PObj, *PObjPoly;
if (!IS_OLST_OBJ(PObjList))
FatalError("GenPOLYObject: Not object list object!");
i = 0;
while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
if (!IS_VEC_OBJ(PObj)) {
WndwInputWindowPutStr("POLY: None vector object found in list, empty object result.");
return NULL;
}
NumVertices++;
}
if (NumVertices < 3) {
WndwInputWindowPutStr("POLY: Less than 3 vertices, empty object result.");
return NULL;
}
PPoly = AllocPolygon(0, 0, VHead = NULL, NULL);
i = 0;
while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
V = AllocVertex(0, 0, NULL, NULL);
PT_COPY(V -> Pt, PObj -> U.Vec);
if (VHead == NULL) {
PPoly -> V = VHead = VTail = V;
}
else {
VTail -> Pnext = V;
VTail = V;
}
}
VTail -> Pnext = VHead; /* Close the vertex list loop. */
/* Update polygon plane equation and vertices normals. */
CGPlaneFrom3Points(PPoly -> Plane, VHead -> Pt, VHead -> Pnext -> Pt,
VHead -> Pnext -> Pnext -> Pt);
V = VHead;
do {
PT_COPY(V -> Normal, PPoly -> Plane);
V = V -> Pnext;
}
while (V != VHead);
PObjPoly = GenPolyObject("", NULL, NULL);/* Gen. the POLY object itself: */
PObjPoly -> U.Pl.P = PPoly;
SetObjectColor(PObjPoly, GlblPrimColor); /* Set its default color. */
return PObjPoly;
}
/*****************************************************************************
* Routine to create an OBJECT directly from set of sspecified polygons. *
* No test is made for the validity of the model in any sense. *
*****************************************************************************/
ObjectStruct *GenObjectFromPolyList(ObjectStruct *PObjList)
{
int i;
PolygonStruct *PPoly,
*PTail = NULL;
ObjectStruct *PObj, *PObjPoly;
if (!IS_OLST_OBJ(PObjList))
FatalError("GenObjectFromPolyList: Not object list object!");
i = 0;
while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
if (!IS_POLY_OBJ(PObj)) {
WndwInputWindowPutStr("POLYTOOBJ: None polygon object found in list, empty object result.");
return NULL;
}
}
PObjPoly = GenPolyObject("", NULL, NULL);/* Gen. the POLY object itself: */
SetObjectColor(PObjPoly, GlblPrimColor); /* Set its default color. */
i = 0;
while ((PObj = PObjList -> U.PObjList[i]) != NULL && i++ < MAX_OBJ_LIST) {
PPoly = CopyPolygonList(PObj -> U.Pl.P);
if (PTail == NULL) {
PObjPoly -> U.Pl.P = PPoly;
}
else {
PTail -> Pnext = PPoly;
}
for (PTail = PPoly; PTail -> Pnext != NULL; PTail = PTail -> Pnext);
}
return PObjPoly;
}
/*****************************************************************************
* Routine to a cross section polygon, by interactively allowing the user *
* to create it. This polygon is returned as a one polygon object. This *
* polygon is mainly for Surface of Revolution (SURFREV) and extrusion *
* (EXTRUDE) operations, although it may be used as an open object by itself. *
*****************************************************************************/
ObjectStruct *GenCROSSECObject(ObjectStruct *PObj)
{
if (PObj && !IS_POLY_OBJ(PObj))
FatalError("CrossSec: operation on non polygonal object");
WndwInputWindowPutStr("GenCrossSecObject not implemented");
return NULL;
}
/*****************************************************************************
* Routine to a create surface of revolution by rotating the given cross *
* section along Z axes. *
* If input is a polyline/gon, it must not be perpendicular to Z. *
*****************************************************************************/
ObjectStruct *GenSURFREVObject(ObjectStruct *Cross)
{
int Resolution, i, j;
RealType XYSize;
MatrixType Mat; /* Rotation Matrix around Z axes. */
VertexStruct *V, *V1, *V1Head, *V2, *V2Head, *VIn, *VInHead;
PolygonStruct *Pl1, *Pl2, *PlIn, *PlNew = NULL;
ObjectStruct *PSurfRev;
if (!IS_POLY_OBJ(Cross) && !IS_CRV_OBJ(Cross)) {
WndwInputWindowPutStr("SurfRev: cross section is not poly/crv. Empty object result");
return NULL;
}
if (IS_POLY_OBJ(Cross)) {
if (APX_EQ(Cross -> U.Pl.P -> Plane[0], 0.0) &&
APX_EQ(Cross -> U.Pl.P -> Plane[1], 0.0)) {
WndwInputWindowPutStr("SurfRev: cross-section perpendicular to Z. Empty object result");
return NULL;
}
Pl1 = AllocPolygon(0, 0, V1Head = CopyVList(Cross -> U.Pl.P -> V), NULL);
PLANE_COPY(Pl1 -> Plane, Cross -> U.Pl.P -> Plane);
Pl2 = AllocPolygon(0, 0, V2Head = CopyVList(Cross -> U.Pl.P -> V), NULL);
PLANE_COPY(Pl2 -> Plane, Cross -> U.Pl.P -> Plane);
PlIn = GenInsidePoly(Pl1);
VInHead = PlIn -> V;
Resolution = GetResolution(TRUE);/* Get refinement factor of object. */
MatGenMatRotZ1(2.0 * M_PI / Resolution, Mat);
for (i = 0; i < Resolution; i++)
{
V2 = V2Head;
do {
MatMultVecby4by4(V2 -> Pt, V2 -> Pt , Mat);
V2 = V2 -> Pnext;
}
while (V2 != NULL && V2 != V2Head);
V1 = V1Head;
if (i < Resolution - 1) /* If this is last loop use the original */
V2 = V2Head; /* poly as we might accumulate error during the */
else /* transformations along the circle. */
V2 = Cross -> U.Pl.P -> V;
VIn = VInHead;
do {
PlNew = GenPolygon4Vrtx(V1 -> Pt, V1 -> Pnext -> Pt,
V2 -> Pnext -> Pt, V2 -> Pt,
VIn -> Pt, PlNew);
/* Update normals: */
for (j = 0, V = PlNew -> V; j < 4; j++, V = V -> Pnext) {
V -> Normal[0] = V -> Pt[0];
V -> Normal[1] = V -> Pt[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 -> Pt[0]) + SQR(V -> Pt[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 -> Pt, V1 -> Pt , Mat);
V1 = V1 -> Pnext;
}
while (V1 != NULL && V1 != V1Head);
VIn = VInHead;
do {
MatMultVecby4by4(VIn -> Pt, VIn -> Pt , Mat);
VIn = VIn -> Pnext;
}
while (VIn != NULL && VIn != VInHead);
}
MyFree((char *) PlIn, ALLOC_POLYGON);
MyFree((char *) Pl1, ALLOC_POLYGON);
MyFree((char *) Pl2, ALLOC_POLYGON);
PSurfRev = GenPolyObject("", NULL, NULL);
PSurfRev -> U.Pl.P = PlNew;
SetObjectColor(PSurfRev, GlblPrimColor); /* Set its default color. */
return PSurfRev;
}
else if (IS_CRV_OBJ(Cross)) {
if (CAGD_NUM_OF_PT_COORD(Cross->U.Crv.Crv -> PType) < 3) {
WndwInputWindowPutStr("SurfRev: cross-section perpendicular to Z. Empty object result");
return NULL;
}
PSurfRev = GenSrfObject("", CagdSurfaceRev(Cross->U.Crv.Crv), NULL);
SetObjectColor(PSurfRev, GlblPrimColor); /* Set its default color. */
return PSurfRev;
}
else
return NULL;
}
/*****************************************************************************
* Routine to a create surface of extrusion out of the given cross section *
* and the given direction. A NULL object will be returned, if the direction *
* vector is in the cross section plane. *
*****************************************************************************/
ObjectStruct *GenEXTRUDEObject(ObjectStruct *Cross, VectorType Dir)
{
int i;
RealType R;
CagdVecStruct CagdDir;
VertexStruct *V1, *V1Head, *V2, *VIn;
PolygonStruct *PBase1, *PBase2, *Pl, *PlIn;
ObjectStruct *PExtrude;
if (!IS_POLY_OBJ(Cross) && !IS_CRV_OBJ(Cross)) {
WndwInputWindowPutStr("Extrude: cross section is not poly/crv. Empty object result");
return NULL;
}
if (IS_POLY_OBJ(Cross)) {
R = DOT_PROD(Cross -> U.Pl.P -> Plane, Dir);
if (APX_EQ(R, 0.0)) {
WndwInputWindowPutStr("Extrusion direction in cross-section plane. Empty object result");
return NULL;
}
/* Prepare two bases (update their plane normal to point INDISE): */
PBase1 = AllocPolygon(0, 0, CopyVList(Cross -> U.Pl.P -> V), NULL);
Pl = PBase2 = AllocPolygon(0, 0, CopyVList(Cross -> U.Pl.P -> V), PBase1);
V1 = V1Head = PBase2 -> V;
do {
PT_ADD(V1 -> Pt, Dir, V1 -> Pt);
V1 = V1 -> Pnext;
}
while (V1 != NULL && V1 != V1Head);
if (R > 0.0) {
PLANE_COPY(PBase1 -> Plane, Cross -> U.Pl.P -> Plane);
for (i = 0; i < 3; i++)
PBase2 -> Plane[i] = (-Cross -> U.Pl.P -> Plane[i]);
PBase2 -> Plane[3] = (-DOT_PROD(PBase2 -> Plane, PBase2 -> V -> Pt));
}
else {
for (i = 0; i < 4; i++)
PBase1 -> Plane[i] = (-Cross -> U.Pl.P -> Plane[i]);
PLANE_COPY(PBase2 -> Plane, Cross -> U.Pl.P -> Plane);
PBase2 -> Plane[3] = (-DOT_PROD(PBase2 -> Plane, PBase2 -> V -> Pt));
}
/* Now generate all the 4 corner polygon between the two bases: */
V1 = V1Head = PBase1 -> V;
V2 = PBase2 -> V;
PlIn = GenInsidePoly(PBase1);
VIn = PlIn -> V;
do {
Pl = GenPolygon4Vrtx(V1 -> Pt, V1 -> Pnext -> Pt,
V2 -> Pnext -> Pt, V2 -> Pt, VIn -> Pt, Pl);
VIn = VIn -> Pnext;
V1 = V1 -> Pnext;
V2 = V2 -> Pnext;
}
while (V1 -> Pnext != NULL && V1 != V1Head);
MyFree((char *) PlIn, ALLOC_POLYGON);
PExtrude = GenPolyObject("", NULL, NULL);
PExtrude -> U.Pl.P = Pl;
SetObjectColor(PExtrude, GlblPrimColor); /* Set its default color. */
/* Update all the polygon vertices normals. */
for (Pl = PExtrude -> U.Pl.P; Pl != NULL; Pl = Pl -> Pnext) {
V1 = V1Head = Pl -> V;
do {
PT_COPY(V1 -> Normal, Pl -> Plane);
V1 = V1 -> Pnext;
}
while (V1 != NULL && V1 != V1Head);
}
return PExtrude;
}
else if (IS_CRV_OBJ(Cross)) {
for (i = 0; i < 3; i++) CagdDir.Vec[i] = Dir[i];
PExtrude = GenSrfObject("", CagdExtrudeSrf(Cross->U.Crv.Crv, &CagdDir),
NULL);
SetObjectColor(PExtrude, GlblPrimColor); /* Set its default color. */
return PExtrude;
}
else
return NULL;
}
/*****************************************************************************
* 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. *
*****************************************************************************/
static PolygonStruct *GenInsidePoly(PolygonStruct *Pl)
{
int Axes;
RealType Dx, Dy;
PointType Pt;
MatrixType Mat;
PolygonStruct *PlIn;
VertexStruct *VHead, *V, *Vnext, *VInHead, *VIn = NULL;
PlIn = AllocPolygon(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);
VHead = V = CopyVList(Pl -> V); /* We dont want to modify original! */
Pl = AllocPolygon(0, 0, VHead, NULL);
do {
MatMultVecby4by4(V -> Pt, V -> Pt, Mat);
V = V -> Pnext;
}
while (V != NULL && V != VHead);
V = VHead;
do {
Vnext = V -> Pnext;
Dx = ABS(V -> Pt[0] - Vnext -> Pt[0]);
Dy = ABS(V -> Pt[1] - Vnext -> Pt[1]);
Pt[0] = (V -> Pt[0] + Vnext -> Pt[0]) / 2.0;/* Prepare middle point. */
Pt[1] = (V -> Pt[1] + Vnext -> Pt[1]) / 2.0;
Pt[2] = V -> Pt[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 = AllocVertex(0, 0, NULL, NULL);
VIn = VIn -> Pnext;
}
else {
PlIn -> V = VInHead = VIn = AllocVertex(0, 0, NULL, NULL);
}
PT_COPY(VIn ->Pt, Pt);
V = Vnext;
}
while (V != NULL && V != VHead);
VIn -> Pnext = VInHead;
MyFree((char *) Pl, ALLOC_POLYGON);/* Free copied (and trans.) vrtx list.*/
/* Transform PlIn to the plane where original Pl is... */
if (!MatInverseMatrix(Mat, Mat)) /* Find the inverse matrix. */
FatalError("GenInsidePoly: Inverse matrix does not exits");
VIn = VInHead;
do {
MatMultVecby4by4(VIn -> Pt, VIn -> Pt, Mat);
VIn = VIn -> Pnext;
}
while (VIn != NULL && VIn != VInHead);
return PlIn;
}
/*****************************************************************************
* 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 normal direction can be evaluated uniquely. *
* No test is made to make sure the 4 points are co-planar... *
* The points are placed into a linear line in order. *
*****************************************************************************/
static PolygonStruct *GenPolygon4Vrtx(VectorType V1, VectorType V2,
VectorType V3, VectorType V4, VectorType Vin, PolygonStruct *Pnext)
{
PolygonStruct *PPoly;
VertexStruct *V;
PPoly = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), Pnext);
PT_COPY(V -> Pt, V1);
V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
PT_COPY(V -> Pt, V2);
V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
PT_COPY(V -> Pt, V3);
V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
PT_COPY(V -> Pt, V4);
V -> Pnext = PPoly -> V; /* Make the Vertex list circular. */
UpdatePolyPlane(PPoly, Vin); /* Update plane equation. */
SET_CONVEX_POLY(PPoly); /* Mark it as convex polygon. */
return PPoly;
}
/*****************************************************************************
* Routine to create a polygon out of a list of 3 vertices V1/2/3 *
* The forth 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 normal direction can be evaluated uniquely. *
* The points are placed into a linear line in order. *
*****************************************************************************/
static PolygonStruct *GenPolygon3Vrtx(VectorType V1, VectorType V2,
VectorType V3, VectorType Vin, PolygonStruct *Pnext)
{
PolygonStruct *PPoly;
VertexStruct *V;
PPoly = AllocPolygon(0, 0, V = AllocVertex(0, 0, NULL, NULL), Pnext);
PT_COPY(V -> Pt, V1);
V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
PT_COPY(V -> Pt, V2);
V -> Pnext = AllocVertex(0, 0, NULL, NULL); V = V -> Pnext;
PT_COPY(V -> Pt, V3);
V -> Pnext = PPoly -> V; /* Make the Vertex list circular. */
UpdatePolyPlane(PPoly, Vin); /* Update plane equation. */
SET_CONVEX_POLY(PPoly); /* Mark it as convex polygon. */
return PPoly;
}
/*****************************************************************************
* Routine to preper transformation martix to do the following (in this *
* order): scale by Scale, rotate such that the Z axis is in direction Dir *
* and then translate by Trans. *
* Algorithm: given the Trans vector, it forms the 4th line of Mat. Dir is *
* used to form the second line (the first 3 lines set the rotation), and *
* finally Scale is used to scale first 3 lines/columns to the needed scale: *
* | Tx Ty Tz 0 | A transformation which takes the coord *
* | Bx By Bz 0 | system into T, N & B as required and *
* [X Y Z 1] * | Nx Ny Nz 0 | then translate it to C. T, N, B are *
* | Cx Cy Cz 1 | scaled by Scale. *
* N is exactly Dir (unit vec) but we got freedom on T & B - T & B must be on *
* a plane perpendicular to N and perpendicular between them but thats all! *
* T is therefore selected using this (heuristic ?) algorithm: *
* Let P be the axis of which the absolute N coefficient is the smallest. *
* Let B be (N cross P) and T be (B cross N). *
*****************************************************************************/
static void GenTransformMatrix(MatrixType Mat, VectorType Trans,
VectorType Dir, RealType Scale)
{
int i, j;
RealType R;
VectorType DirN, T, B, P;
MatrixType TempMat;
PT_COPY(DirN, Dir);
PT_NORMALIZE(DirN);
PT_CLEAR(P);
for (i = 1, j = 0, R = ABS(DirN[0]); i < 3; i++) if (R > ABS(DirN[i])) {
R = DirN[i];
j = i;
}
P[j] = 1.0;/* Now P is set to the axis with the biggest angle from DirN. */
VecCrossProd(B, DirN, P); /* Calc the bi-normal. */
VecCrossProd(T, B, DirN); /* Calc the tangent. */
MatGenUnitMat(Mat);
for (i = 0; i < 3; i++) {
Mat[0][i] = T[i];
Mat[1][i] = B[i];
Mat[2][i] = DirN[i];
}
MatGenMatScale(Scale, Scale, Scale, TempMat);
MatMultTwo4by4(Mat, TempMat, Mat);
MatGenMatTrans(Trans[0], Trans[1], Trans[2], TempMat);
MatMultTwo4by4(Mat, Mat, TempMat);
}
/*****************************************************************************
* Routine to update the Plane equation of the given polygon such that the *
* Vin vertex will be in the positive side of it. *
* It is assumed the polygon has at list 3 points... *
*****************************************************************************/
void UpdatePolyPlane(PolygonStruct *PPoly, VectorType Vin)
{
int i;
VectorType V1, V2;
RealType Len;
VertexStruct *V;
V = PPoly -> V; PT_SUB(V1, V -> Pt, V -> Pnext -> Pt);
V = V -> Pnext; PT_SUB(V2, V -> Pt, V -> Pnext -> Pt);
VecCrossProd(PPoly -> Plane, V1, V2); /* Find Plane Normal. */
/* Normalize the plane such that the normal has length of 1: */
Len = PT_LENGTH(PPoly -> Plane);
for (i = 0; i < 3; i++) PPoly -> Plane[i] /= Len;
PPoly -> Plane[3] = (-DOT_PROD(PPoly -> Plane, PPoly -> V -> Pt));
if (DOT_PROD(PPoly -> Plane, Vin) + PPoly -> Plane[3] < 0) {
/* Flip plane normal and reverse the vertex list. */
ReverseVrtxList(PPoly);
for (i = 0; i < 4; i++) PPoly -> Plane[i] = (-PPoly -> Plane[i]);
}
}
/*****************************************************************************
* 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. *
*****************************************************************************/
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);
VecCrossProd(V2, V1, Normal);
VecCrossProd(Normal, V2, V1);
}
PT_NORMALIZE(Normal);
}
