home *** CD-ROM | disk | FTP | other *** search
/ Gold Fish 3 / goldfish_volume_3.bin / files / gfx / 3d / irit / geom_lib / convex.c < prev    next >
Encoding:
C/C++ Source or Header  |  1995-02-24  |  23.7 KB  |  609 lines
  1. /******************************************************************************
  2. * Convex.c - test convexity and converts polygons to convex ones.          *
  3. *******************************************************************************
  4. * Written by Gershon Elber, March 1990.                          *
  5. ******************************************************************************/
  6.  
  7. /* #define DEBUG2           Defines some printing/debugging routines. */
  8.  
  9. #include <stdio.h>
  10. #include <ctype.h>
  11. #include <math.h>
  12. #include <string.h>
  13. #include "allocate.h"
  14. #include "convex.h"
  15. #include "poly_cln.h"
  16. #include "geomat3d.h"
  17. #include "intrnrml.h"
  18. #include "priorque.h"
  19.  
  20. /* Used to hold edges (V | V -> Pnext) that intersect with level y = Ylevel. */
  21. typedef struct InterYVrtxList {
  22.     ByteType InterYType;
  23.     IPVertexStruct *V;
  24.     struct InterYVrtxList *Pnext;
  25. } InterYVrtxList;
  26.  
  27. #define CONVEX_EPSILON  1e-8       /* Colinearity of two normalized vectors. */
  28.  
  29. #define    INTER_Y_NONE    0               /* Y level intersection type. */
  30. #define    INTER_Y_START    1
  31. #define    INTER_Y_MIDDLE    2
  32.  
  33. #define LOOP_ABOVE_Y    0      /* Type of open loops extracted from polygon. */
  34. #define LOOP_BELOW_Y    1
  35.  
  36. /* Used to sort and combine the polygons above Ylevel together if possible.  */
  37. typedef struct SortPolysInX {
  38.     IPVertexStruct *VMinX, *VMaxX;
  39.     int Reverse;      /* If TRUE than VMinX vertex is AFTER VMaxX vertex. */
  40. } SortPolysInX;
  41.  
  42. /* The following are temporary flags used to mark vertices that were visited */
  43. /* by the loop tracing, at list once. As each vertex may be visited two at   */
  44. /* the most (as starting & as end point of open loop), this is enough.         */
  45. /* INTER_TAG is used to mark vertices that created brand new with intersected*/
  46. /* with the line y = Ylevel. Those are used to detect INTERNAL edges - if    */
  47. /* at list one end of it is INTER_TAG, that edge is INTERNAL.             */
  48. #define INTER_TAG   0x40
  49. #define VISITED_TAG 0x80
  50.  
  51. #define    IS_INTER_VRTX(Vrtx)    ((Vrtx)->Tags & INTER_TAG)
  52. #define    SET_INTER_VRTX(Vrtx)    ((Vrtx)->Tags |= INTER_TAG)
  53. #define    RST_INTER_VRTX(Vrtx)    ((Vrtx)->Tags &= ~INTER_TAG)
  54.  
  55. #define    IS_VISITED_VRTX(Vrtx)    ((Vrtx)->Tags & VISITED_TAG)
  56. #define    SET_VISITED_VRTX(Vrtx)    ((Vrtx)->Tags |= VISITED_TAG)
  57. #define    RST_VISITED_VRTX(Vrtx)    ((Vrtx)->Tags &= ~VISITED_TAG)
  58.  
  59. static int SplitPolyIntoTwo(IPPolygonStruct *Pl,
  60.                 IPVertexStruct *V,
  61.                 IPPolygonStruct **Pl1,
  62.                 IPPolygonStruct **Pl2);
  63. static IPVertexStruct *FindRayPolyInter(IPPolygonStruct *Pl,
  64.                     IPVertexStruct *VRay,
  65.                     PointType RayDir,
  66.                     PointType PInter);
  67. static void TestConvexityDir(IPPolygonStruct *Pl);
  68.  
  69. #ifdef DEBUG2
  70. static void PrintVrtxList(IPVertexStruct *V);
  71. static void PrintPoly(IPPolygonStruct *P);
  72. #endif /* DEBUG2 */
  73.  
  74. /*****************************************************************************
  75. * DESCRIPTION:                                                               M
  76. *   Routine to prepare a transformation martix to rotate such that Dir is    M
  77. * parallel to the Z axes. Used by the convex decomposition to rotate the     M
  78. * polygons to be XY plane parallel.                         M
  79. *    Algorithm: form a 4 by 4 matrix from Dir as follows:             M
  80. *                |  Tx  Ty  Tz  0 |   A transformation which takes the coord V
  81. *                |  Bx  By  Bz  0 |  system into T, N & B as required.         V
  82. * [X  Y  Z  1] * |  Nx  Ny  Nz  0 |                         V
  83. *                |  0   0   0   1 |                         V
  84. *   N is exactly Dir, but we got freedom on T & B which must be on         M
  85. * a plane perpendicular to N and perpendicular between them but thats all!   M
  86. *   T is therefore selected using this (heuristic ?) algorithm:             M
  87. *   Let P be the axis of which the absolute N coefficient is the smallest.   M
  88. *   Let B be (N cross P) and T be (B cross N).                     M
  89. *                                                                            *
  90. * PARAMETERS:                                                                M
  91. *   Mat:     To place the constructed homogeneous transformation.            M
  92. *   Dir:     To derive a transformation such that Dir goes to Z axis.        M
  93. *                                                                            *
  94. * RETURN VALUE:                                                              M
  95. *   void                                                                     M
  96. *                                                                            *
  97. * KEYWORDS:                                                                  M
  98. *   GenRotateMatrix, transformations                                         M
  99. *****************************************************************************/
  100. void GenRotateMatrix(MatrixType Mat, VectorType Dir)
  101. {
  102.     int i, j;
  103.     RealType R;
  104.     VectorType DirN, T, B, P;
  105.  
  106.     PT_COPY(DirN, Dir);
  107.     PT_NORMALIZE(DirN);
  108.     PT_CLEAR(P);
  109.     for (i = 1, j = 0, R = ABS(DirN[0]); i < 3; i++)
  110.     if (R > ABS(DirN[i])) {
  111.         R = DirN[i];
  112.         j = i;
  113.     }
  114.     P[j] = 1.0;/* Now P is set to the axis with the biggest angle from DirN. */
  115.  
  116.     GMVecCrossProd(B, DirN, P);                  /* calc the bi-normal. */
  117.     PT_NORMALIZE(B);
  118.     GMVecCrossProd(T, B, DirN);                /* calc the tangent. */
  119.  
  120.     MatGenUnitMat(Mat);
  121.     for (i = 0; i < 3; i++) {
  122.     Mat[i][0] = T[i];
  123.     Mat[i][1] = B[i];
  124.     Mat[i][2] = DirN[i];
  125.     }
  126. }
  127.  
  128. /*****************************************************************************
  129. * DESCRIPTION:                                                               M
  130. *   Routine to test all polygons in a given object for convexity, and split  M
  131. * non convex ones, non destructively - the original object is not modified.  M
  132. *   This function will introduce new vertices to the split polygons.         M
  133. *                                                                            *
  134. * PARAMETERS:                                                                M
  135. *   PObj:       To test for convexity of its polygons.                       M
  136. *                                                                            *
  137. * RETURN VALUE:                                                              M
  138. *   IPObjectStruct *:  A duplocate of PObj, but with convex polygons only.   M
  139. *                                                                            *
  140. * KEYWORDS:                                                                  M
  141. *   ConvexPolyObjectN, convexity, convex polygon                             M
  142. *****************************************************************************/
  143. IPObjectStruct *ConvexPolyObjectN(IPObjectStruct *PObj)
  144. {
  145.     IPObjectStruct
  146.     *PObjCopy= CopyObject(NULL, PObj, FALSE);
  147.  
  148.     ConvexPolyObject(PObjCopy);
  149.  
  150.     return PObjCopy;
  151. }
  152.  
  153. /*****************************************************************************
  154. * DESCRIPTION:                                                               M
  155. *   Routine to test all polygons in a given object for convexity, and split  M
  156. * non convex ones, in place.                                                 M
  157. *   This function will introduce new vertices to the split polygons.         M
  158. *                                                                            *
  159. * PARAMETERS:                                                                M
  160. *   PObj:       To test for convexity of its polygons, and split into convex M
  161. *               polygons non convex polygons found, in place. Either a       M
  162. *        polygonal object or a list of polygonal objects.         M
  163. *                                                                            *
  164. * RETURN VALUE:                                                              M
  165. *   void                                                                 M
  166. *                                                                            *
  167. * KEYWORDS:                                                                  M
  168. *   ConvexPolyObject, convexity, convex polygon                              M
  169. *****************************************************************************/
  170. void ConvexPolyObject(IPObjectStruct *PObj)
  171. {
  172.     IPPolygonStruct *Pl, *PlSplit, *PlTemp,
  173.     *PlPrev = NULL;
  174.  
  175.     if (IP_IS_OLST_OBJ(PObj)) {
  176.     int i = 0;
  177.     IPObjectStruct *Obj;
  178.  
  179.     while ((Obj = ListObjectGet(PObj, i++)) != NULL)
  180.         ConvexPolyObject(Obj);
  181.  
  182.     return;
  183.     }
  184.  
  185.     if (!IP_IS_POLY_OBJ(PObj) || IP_IS_POLYLINE_OBJ(PObj))
  186.         return;
  187.  
  188.     Pl = PObj -> U.Pl;
  189.     while (Pl != NULL) {
  190.     if (!ConvexPolygon(Pl)) {
  191.         PlSplit = SplitNonConvexPoly(Pl);
  192.         CleanUpPolygonList(&PlSplit);       /* Zero length edges etc. */
  193.         if (PlSplit != NULL) {     /* Something is wrong here, ignore. */
  194.         if (Pl == PObj -> U.Pl)
  195.             PObj -> U.Pl = PlSplit;               /* First. */
  196.         else
  197.             PlPrev -> Pnext = PlSplit;
  198.  
  199.         PlTemp = PlSplit;
  200.         while (PlTemp -> Pnext != NULL)
  201.             PlTemp = PlTemp -> Pnext;
  202.         PlTemp -> Pnext = Pl -> Pnext;
  203.         PlPrev = PlTemp;
  204.         IPFreePolygon(Pl);            /* Free old polygon. */
  205.         Pl = PlPrev -> Pnext;
  206.         }
  207.         else {
  208.         if (Pl == PObj -> U.Pl)
  209.             PObj -> U.Pl = Pl -> Pnext;
  210.         PlPrev = Pl -> Pnext;
  211.         IPFreePolygon(Pl);            /* Free old polygon. */
  212.         Pl = PlPrev;
  213.         }
  214.     }
  215.     else {
  216.         PlPrev = Pl;
  217.         Pl = Pl -> Pnext;
  218.     }
  219.     }
  220. }
  221.  
  222. /*****************************************************************************
  223. * DESCRIPTION:                                                               M
  224. *   Routine to test if the given polygon is convex or not.             M
  225. * Algorithm: The polygon is convex iff the normals generated from cross      M
  226. * products of two consecutive edges points to the same direction.            M
  227. *   Note a 5 star polygon satisfies this constraint but it is self           M
  228. * intersectingand we assume given polygon is not selft intersecting.         M
  229. *   The computed direction is alos verified against the polygon's plane      M
  230. * normal.                                     M
  231. *   The routine returns TRUE iff the polygon is convex. In addition the      M
  232. * polygon CONVEX tag (see IPPolygonStruct) is also updated.             M
  233. *   If the polygon is already marked as convex, nothing is tested!         M
  234. *                                                                            *
  235. * PARAMETERS:                                                                M
  236. *   Pl:        To test its convexity condition.                              M
  237. *                                                                            *
  238. * RETURN VALUE:                                                              M
  239. *   int:       TRUE if convex, FALSE otherwise.                              M
  240. *                                                                            *
  241. * KEYWORDS:                                                                  M
  242. *   ConvexPolygon, convexity, convex polygon                                 M
  243. *****************************************************************************/
  244. int ConvexPolygon(IPPolygonStruct *Pl)
  245. {
  246.     int FirstTime = TRUE;
  247.     RealType Size,
  248.     NormalSign = 0.0;
  249.     VectorType V1, V2, PolyNormal, Normal;
  250.     IPVertexStruct *VNext,
  251.     *V = Pl -> PVertex;
  252.  
  253.     if (IP_IS_CONVEX_POLY(Pl))
  254.     return TRUE;                   /* Nothing to do around here. */
  255.  
  256.     /* Copy only A, B, C from Ax+By+Cz+D = 0: */
  257.     PT_COPY(PolyNormal, Pl -> Plane);
  258.  
  259.     do {
  260.     VNext = V -> Pnext;
  261.  
  262.     PT_SUB(V1, VNext -> Coord, V -> Coord);
  263.     if ((Size = PT_LENGTH(V1)) > IRIT_EPSILON) {
  264.         Size = 1.0 / Size;
  265.         PT_SCALE(V1, Size);
  266.     }
  267.     PT_SUB(V2, VNext -> Pnext -> Coord, VNext -> Coord);
  268.     if ((Size = PT_LENGTH(V2)) > IRIT_EPSILON) {
  269.         Size = 1.0 / Size;
  270.         PT_SCALE(V2, Size);
  271.     }
  272.     GMVecCrossProd(Normal, V1, V2);
  273.  
  274.     if (PT_LENGTH(Normal) < CONVEX_EPSILON) {
  275.         V = VNext;
  276.         continue;                    /* Skip colinear points. */
  277.     }
  278.  
  279.     if (FirstTime) {
  280.         FirstTime = FALSE;
  281.         NormalSign = DOT_PROD(Normal, PolyNormal);
  282.     }
  283.     else if (NormalSign * DOT_PROD(Normal, PolyNormal) < 0.0) {
  284.         IP_RST_CONVEX_POLY(Pl);
  285.         return FALSE;          /* Different signs --> not convex. */
  286.     }
  287.  
  288.     V = VNext;
  289.     }
  290.     while (V != Pl -> PVertex);
  291.  
  292.     IP_SET_CONVEX_POLY(Pl);
  293.  
  294.     if (NormalSign < 0.0)
  295.         IritPrsrReverseVrtxList(Pl);
  296.  
  297.     return TRUE;    /* All signs are the same --> the polygon is convex. */
  298. }
  299.  
  300. /*****************************************************************************
  301. * DESCRIPTION:                                                               M
  302. *   Routine to split non convex polygon into a list of convex ones.         M
  303. * 1. Remove a polygon from GlblList. If non exists stop.             M
  304. * 2. Search for non convex corner. If not found stop - polygon is convex.    M
  305. *    Otherwise let the non convex polygon found be V(i).             M
  306. * 3. Fire a ray from V(i) in the opposite direction to V(i-1). Find the      M
  307. *    closest intersection of the ray with polygon boundary P.             M
  308. * 4. Split the polygon into two at V(i)-P edge and push the two new polygons M
  309. *    on the GlblList.                                 M
  310. * 5. Goto 1.                                     M
  311. *                                                                            *
  312. * PARAMETERS:                                                                M
  313. *   Pl:        Non convex polygon to split into convex ones.                 M
  314. *                                                                            *
  315. * RETURN VALUE:                                                              M
  316. *   IPPolygonStruct *:  A list of convex polygons resulting from splitting   M
  317. *                       up Pl.                                               M
  318. *                                                                            *
  319. * KEYWORDS:                                                                  M
  320. *   SplitNonConvexPoly, convexity, convex polygon                            M
  321. *****************************************************************************/
  322. IPPolygonStruct *SplitNonConvexPoly(IPPolygonStruct *Pl)
  323. {
  324.     int IsConvex;
  325.     RealType Size;
  326.     IPPolygonStruct *GlblList, *Pl1, *Pl2,
  327.     *GlblSplitPl = NULL;
  328.     VectorType V1, V2, PolyNormal, Normal;
  329.     IPVertexStruct *V, *VNext;
  330.  
  331.     TestConvexityDir(Pl);
  332.  
  333.     GlblList = IPAllocPolygon(0, 0, CopyVertexList(Pl -> PVertex), NULL);
  334.     PLANE_COPY(GlblList -> Plane, Pl -> Plane);
  335.  
  336.     /* Copy only A, B, C from Ax+By+Cz+D = 0 plane equation: */
  337.     PT_COPY(PolyNormal, Pl -> Plane);
  338.  
  339.     while (GlblList != NULL) {
  340.     Pl = GlblList;
  341.     GlblList = GlblList -> Pnext;
  342.     Pl -> Pnext = NULL;
  343.  
  344.     IsConvex = TRUE;
  345.     V = Pl -> PVertex;
  346.     do {
  347.         VNext = V -> Pnext;
  348.  
  349.         PT_SUB(V1, VNext -> Coord, V -> Coord);
  350.          if ((Size = PT_LENGTH(V1)) > IRIT_EPSILON) {
  351.         Size = 1.0 / Size;
  352.         PT_SCALE(V1, Size);
  353.         }
  354.         PT_SUB(V2, VNext -> Pnext -> Coord, VNext -> Coord);
  355.          if ((Size = PT_LENGTH(V2)) > IRIT_EPSILON) {
  356.         Size = 1.0 / Size;
  357.         PT_SCALE(V2, Size);
  358.         }
  359.         GMVecCrossProd(Normal, V1, V2);
  360.         if (PT_LENGTH(Normal) < CONVEX_EPSILON) {
  361.         V = VNext;
  362.         continue;                /* Skip colinear points. */
  363.         }
  364.  
  365.         if (DOT_PROD(Normal, PolyNormal) < 0.0 &&
  366.         SplitPolyIntoTwo(Pl, V, &Pl1, &Pl2)) {
  367.         Pl -> PVertex = NULL; /* Dont free vertices - used in Pl1/2. */
  368.         IPFreePolygon(Pl);
  369.  
  370.         Pl1 -> Pnext = GlblList;       /* Push polygons on GlblList. */
  371.         GlblList = Pl1;
  372.         Pl2 -> Pnext = GlblList;
  373.         GlblList = Pl2;
  374.  
  375.         IsConvex = FALSE;
  376.         break;
  377.         }
  378.  
  379.         V = VNext;
  380.     }
  381.     while (V != Pl -> PVertex);
  382.  
  383.     if (IsConvex) {
  384.         IP_SET_CONVEX_POLY(Pl);
  385.         Pl -> Pnext = GlblSplitPl;
  386.         GlblSplitPl = Pl;
  387.     }
  388.     }
  389.  
  390.     return GlblSplitPl;
  391. }
  392.  
  393. /*****************************************************************************
  394. * DESCRIPTION:                                                               *
  395. *   Split polygon Pl at the vertex specified by V -> Pnext, given V, into    *
  396. * two, by firing a ray from V -> Pnext in the opposite direction to V and    *
  397. * finding closest intersection with the polygon Pl. (V -> Pnext, P) is the   *
  398. * edge to split the polygon at.                             *
  399. *                                                                            *
  400. * PARAMETERS:                                                                *
  401. *   Pl:        Polygon to split at vertex V -> Pnext.                        *
  402. *   V:         One vertex before the vertex to split Pl at.                  *
  403. *   Pl1, Pl2:  Where the two new polygons should end up.                     *
  404. *                                                                            *
  405. * RETURN VALUE:                                                              *
  406. *   int:       TRUE if succesful, FALSE otehrwise.                           *
  407. *****************************************************************************/
  408. static int SplitPolyIntoTwo(IPPolygonStruct *Pl,
  409.                 IPVertexStruct *V,
  410.                 IPPolygonStruct **Pl1,
  411.                 IPPolygonStruct **Pl2)
  412. {
  413.     PointType Vl, PInter;
  414.     IPVertexStruct *VInter, *VNew1, *VNew2;
  415.  
  416.     PT_SUB(Vl, V -> Pnext -> Coord, V -> Coord);
  417.     VInter = FindRayPolyInter(Pl, V -> Pnext, Vl, PInter);
  418.     V = V -> Pnext;
  419.  
  420.     if (VInter == NULL ||
  421.     VInter == V ||
  422.     VInter -> Pnext == V)
  423.     return FALSE;
  424.  
  425.     /* Make the two polygon vertices lists. */
  426.     VNew1 = IPAllocVertex(V -> Count, V -> Tags, NULL, V -> Pnext);
  427.     PT_COPY(VNew1 -> Coord, V -> Coord);
  428.     PT_COPY(VNew1 -> Normal, V -> Normal);
  429.     IP_SET_INTERNAL_VRTX(V);
  430.     if (PT_APX_EQ(VInter -> Coord, PInter)) {
  431.     /* Intersection points is close to VInter point. */
  432.     VNew2 = IPAllocVertex(VInter -> Count, VInter -> Tags, NULL,
  433.                             VInter -> Pnext);
  434.     PT_COPY(VNew2 -> Coord, VInter -> Coord);
  435.     PT_COPY(VNew2 -> Normal, VInter -> Normal);
  436.     VInter -> Pnext = VNew1;
  437.     IP_SET_INTERNAL_VRTX(VInter);
  438.     V -> Pnext = VNew2;
  439.     }
  440.     else if (PT_APX_EQ(VInter -> Pnext -> Coord, PInter)) {
  441.     /* Intersection points is close to VInter -> Pnext point. */
  442.     VNew2 = IPAllocVertex(VInter -> Pnext -> Count,
  443.                   VInter -> Pnext -> Tags,
  444.                   NULL, VInter -> Pnext -> Pnext);
  445.     PT_COPY(VNew2 -> Coord, VInter -> Pnext -> Coord);
  446.     PT_COPY(VNew2 -> Normal, VInter -> Pnext -> Normal);
  447.     VInter -> Pnext -> Pnext = VNew1;
  448.     IP_SET_INTERNAL_VRTX(VInter -> Pnext);
  449.     V -> Pnext = VNew2;
  450.     }
  451.     else {
  452.     /* PInter is in the middle of (VInter, VInter -> Pnext) edge: */
  453.     VNew2 = IPAllocVertex(VInter -> Count, VInter -> Tags, NULL,
  454.                             VInter -> Pnext);
  455.     PT_COPY(VNew2 -> Coord, PInter);
  456.     VInter -> Pnext = IPAllocVertex(0, 0, NULL, VNew1);
  457.     PT_COPY(VInter -> Pnext -> Coord, PInter);
  458.     InterpNrmlBetweenTwo(VNew2, VInter, VNew2 -> Pnext);
  459.     PT_COPY(VInter -> Pnext -> Normal, VNew2 -> Normal);
  460.     IP_SET_INTERNAL_VRTX(VInter -> Pnext);
  461.     V -> Pnext = VNew2;
  462.     }
  463.  
  464.     *Pl1 = IPAllocPolygon(0, 0, VNew1, NULL);
  465.     PLANE_COPY((*Pl1) -> Plane, Pl -> Plane);
  466.     *Pl2 = IPAllocPolygon(0, 0, VNew2, NULL);
  467.     PLANE_COPY((*Pl2) -> Plane, Pl -> Plane);
  468.  
  469.     return TRUE;
  470. }
  471.  
  472. /*****************************************************************************
  473. * DESCRIPTION:                                                               *
  474. *   Finds where a ray first intersect a given polygon. The ray starts at one *
  475. * of the polygon's vertices and so distance less than EPSILON is ignored.    *
  476. *   Returns the vertex V in which (V, V -> Pnext) has the closest         *
  477. * intersection with the ray PRay, DRay at Inter.                 *
  478. *                                                                            *
  479. * PARAMETERS:                                                                *
  480. *   Pl:         Polygon to test for intersection with ray.                   *
  481. *   VRay:       Origin of ray.                                               *
  482. *   RayDir:     Direction of ray.                                            *
  483. *   PInter:     location of intersection is to be placed herein.             *
  484. *                                                                            *
  485. * RETURN VALUE:                                                              *
  486. *   IPVertexStruct *:   Vertex V such that edge (V, V -> Pnext) has the      *
  487. *                       needed intersection.                                 *
  488. *****************************************************************************/
  489. static IPVertexStruct *FindRayPolyInter(IPPolygonStruct *Pl,
  490.                     IPVertexStruct *VRay,
  491.                     PointType RayDir,
  492.                     PointType PInter)
  493. {
  494.     RealType t1, t2,
  495.     MinT = INFINITY;
  496.     PointType Vl, Ptemp1, Ptemp2, PRay;
  497.     IPVertexStruct *VNext,
  498.     *V = Pl -> PVertex,
  499.     *VInter = NULL;
  500.  
  501.     PT_COPY(PRay, VRay -> Coord);
  502.  
  503.     do {
  504.     VNext = V -> Pnext;
  505.     if (V != VRay && VNext != VRay) {
  506.         PT_SUB(Vl, VNext -> Coord, V -> Coord);
  507.         if (CGDistPointLine(PRay, V -> Coord, Vl) > IRIT_EPSILON) {
  508.         /* Only if the point the ray is shoot from is not on line: */
  509.         CG2PointsFromLineLine(PRay, RayDir, V -> Coord, Vl,
  510.                       Ptemp1, &t1, Ptemp2, &t2);
  511.         if (CGDistPointPoint(Ptemp1, Ptemp2) < IRIT_EPSILON * 10.0 &&
  512.             t1 > IRIT_EPSILON && t1 < MinT &&
  513.             (t2 <= 1.0 || APX_EQ(t2, 1.0)) &&
  514.             (t2 >= 0.0 || APX_EQ(t2, 0.0))) {
  515.             PT_COPY(PInter, Ptemp2);
  516.             VInter = V;
  517.             MinT = t1;
  518.         }
  519.         }
  520.     }
  521.     V = VNext;
  522.     }
  523.     while (V != Pl -> PVertex && V -> Pnext != NULL);
  524.  
  525.     return VInter;
  526. }
  527.  
  528. /*****************************************************************************
  529. * DESCRIPTION:                                                               *
  530. *   Test convexity direction - a cross product of two edges of a convex      *
  531. * corner of the polygon should point in the normal direction. if this is not *
  532. * the case - the polygon vertices are reveresed.                 *
  533. *                                                                            *
  534. * PARAMETERS:                                                                *
  535. *   Pl:         To test of convexity direction.                              *
  536. *                                                                            *
  537. * RETURN VALUE:                                                              *
  538. *   void                                                                     *
  539. *****************************************************************************/
  540. static void TestConvexityDir(IPPolygonStruct *Pl)
  541. {
  542.     int Coord = 0;
  543.     VectorType V1, V2, Normal;
  544.     IPVertexStruct *V, *VExtrem;
  545.  
  546.     /* Find the minimum component direction of the normal and used that axes */
  547.     /* to find an extremum point on the polygon - that extrmum point must be */
  548.     /* a convex corner so we can find the normal direction of convex corner. */
  549.     if (ABS(Pl -> Plane[1]) < ABS(Pl -> Plane[Coord]))
  550.     Coord = 1;
  551.     if (ABS(Pl -> Plane[2]) < ABS(Pl -> Plane[Coord]))
  552.     Coord = 2;
  553.     V = VExtrem = Pl -> PVertex;
  554.     do {
  555.     if (V -> Coord[Coord] > VExtrem -> Coord[Coord])
  556.         VExtrem = V;
  557.     V = V -> Pnext;
  558.     }
  559.     while (V != Pl -> PVertex && V != NULL);
  560.  
  561.     /* Make sure next vertex is not at the extremum value: */
  562.     while (APX_EQ(VExtrem -> Coord[Coord], VExtrem -> Pnext -> Coord[Coord]))
  563.     VExtrem = VExtrem -> Pnext;
  564.  
  565.     /* O.K. V form a convex corner - evaluate its two edges cross product:   */
  566.     for (V = Pl -> PVertex; V -> Pnext != VExtrem; V = V -> Pnext); /* Prev. */
  567.     PT_SUB(V1, VExtrem -> Coord, V -> Coord);
  568.     PT_SUB(V2, VExtrem -> Pnext -> Coord, VExtrem -> Coord);
  569.     GMVecCrossProd(Normal, V1, V2);
  570.  
  571.     if (DOT_PROD(Normal, Pl -> Plane) < 0.0)
  572.     IritPrsrReverseVrtxList(Pl);
  573. }
  574.  
  575. #ifdef DEBUG2
  576.  
  577. /*****************************************************************************
  578. * DESCRIPTION:                                                               *
  579. *   Print the content of the given vertex list, to standard output.         *
  580. *                                                                            *
  581. * PARAMETERS:                                                                *
  582. *   P:          Polygon to print its vertex list.                            *
  583. *                                                                            *
  584. * RETURN VALUE:                                                              *
  585. *   void                                                                     *
  586. *****************************************************************************/
  587. static void PrintPoly(IPPolygonStruct *P)
  588. {
  589.     IPVertexStruct
  590.     *VFirst = P -> PVertex,
  591.     *V = VFirst;
  592.  
  593.     if (V == NULL)
  594.     return;
  595.  
  596.     printf("VERTEX LIST:\n");
  597.     do {
  598.     printf("%12lg %12lg %12lg, Tags = %02x\n",
  599.         V -> Coord[0], V -> Coord[1], V -> Coord[2], V -> Tags);
  600.     V = V -> Pnext;
  601.     }
  602.     while (V != NULL && V != VFirst);
  603.  
  604.     if (V == NULL)
  605.     printf("Loop is not closed (NULL terminated)\n");
  606. }
  607.  
  608. #endif /* DEBUG2 */
  609.