home *** CD-ROM | disk | FTP | other *** search
/ Programmer 7500 / MAX_PROGRAMMERS.iso / INFO / IRIT / IRITS.ZIP / CONVEX.C < prev    next >
Encoding:
C/C++ Source or Header  |  1990-05-05  |  15.9 KB  |  467 lines
  1. /*****************************************************************************
  2. *   "Irit" - the 3d polygonal solid modeller.                     *
  3. *                                         *
  4. * Written by:  Gershon Elber                Ver 0.2, Mar. 1990   *
  5. ******************************************************************************
  6. *   Module to handle convexity of polygons: Test for convexity, and split    *
  7. * none convex polygons into convex ones.                     *
  8. *****************************************************************************/
  9.  
  10. #define DEBUG2         /* Defines some printing/debugging routines. */
  11.  
  12. #ifdef __MSDOS__
  13. #include <alloc.h>
  14. #endif /* __MSDOS__ */
  15.  
  16. #include <stdio.h>
  17. #include <ctype.h>
  18. #include <math.h>
  19. #include <string.h>
  20. #include "program.h"
  21. #include "geomat3d.h"
  22. #include "booleang.h"
  23. #include "objects.h"
  24. #include "allocatg.h"
  25. #include "priorqg.h"
  26. #include "convexg.h"
  27. #include "convexl.h"
  28.  
  29. static void SplitPolyIntoTwo(PolygonStruct *Pl, VertexStruct *V,
  30.                 PolygonStruct **Pl1, PolygonStruct **Pl2);
  31. static VertexStruct *FindRayPolyInter(PolygonStruct *Pl,
  32.         VertexStruct *VRay, PointType RayDir, PointType PInter);
  33. static void TestConvexityDir(PolygonStruct *Pl);
  34.  
  35. #ifdef DEBUG2
  36. static void PrintVrtxList(VertexStruct *V);
  37. #endif /* DEBUG2 */
  38.  
  39. /*****************************************************************************
  40. *   Routine to prepare transformation martix to rotate such that Dir is         *
  41. * parallel to the Z axes. Used by the convex decomposition to rotate the     *
  42. * polygons to be XY plane parallel.                         *
  43. *    Algorithm: form a 4 by 4 matrix from Dir as follows:             *
  44. *                |  Tx  Ty  Tz  0 |   A transformation which takes the coord *
  45. *                |  Bx  By  Bz  0 |  system into T, N & B as required.         *
  46. * [X  Y  Z  1] * |  Nx  Ny  Nz  0 |                         *
  47. *                |  0   0   0   1 |                         *
  48. * N is exactly Dir, but we got freedom on T & B - T & B must be on         *
  49. * a plane perpendicular to N and perpendicular between them but thats all!   *
  50. * T is therefore selected using this (heuristic ?) algorithm:             *
  51. * Let P be the axis of which the absolute N coefficient is the smallest.     *
  52. * Let B be (N cross P) and T be (B cross N).                     *
  53. *****************************************************************************/
  54. void GenRotateMatrix(MatrixType Mat, VectorType Dir)
  55. {
  56.     int i, j;
  57.     RealType R;
  58.     VectorType DirN, T, B, P;
  59.  
  60.     PT_COPY(DirN, Dir);
  61.     PT_NORMALIZE(DirN);
  62.     PT_CLEAR(P);
  63.     for (i=1, j=0, R = ABS(DirN[0]); i<3; i++) if (R > ABS(DirN[i])) {
  64.     R = DirN[i];
  65.     j = i;
  66.     }
  67.     P[j] = 1.0;/* Now P is set to the axis with the biggest angle from DirN. */
  68.  
  69.     VecCrossProd(B, DirN, P);                  /* calc the bi-normal. */
  70.     PT_NORMALIZE(B);
  71.     VecCrossProd(T, B, DirN);                /* calc the tangent. */
  72.  
  73.     MatGenUnitMat(Mat);
  74.     for (i=0; i<3; i++) {
  75.     Mat[i][0] = T[i];
  76.     Mat[i][1] = B[i];
  77.     Mat[i][2] = DirN[i];
  78.     }
  79. }
  80.  
  81. /*****************************************************************************
  82. *   Routine to test all polygons in a given object for convexity, and split  *
  83. * non convex ones, non destructively - the original object is not modified.  *
  84. *****************************************************************************/
  85. ObjectStruct *ConvexPolyObjectN(ObjectStruct *PObj)
  86. {
  87.     ObjectStruct *PObjCopy;
  88.  
  89.     PObjCopy = CopyObject(NULL, PObj, FALSE);
  90.     ConvexPolyObject(PObjCopy);
  91.  
  92.  
  93.     return PObjCopy;
  94. }
  95.  
  96. /*****************************************************************************
  97. *   Routine to test all the polygons in a given object for convexity, and    *
  98. * split the non convex ones.                             *
  99. *****************************************************************************/
  100. void ConvexPolyObject(ObjectStruct *PObj)
  101. {
  102.     PolygonStruct *Pl, *PlPrev = NULL, *PlSplit, *PlTemp;
  103.  
  104.     if (!IS_GEOM_OBJ(PObj))
  105.     FatalError("ConvexPolyObject: parameter given is not geometric object");
  106.  
  107.     Pl = PObj -> U.Pl;
  108.     while (Pl != NULL) {
  109.     if (!ConvexPolygon(Pl)) {
  110.         PlSplit = SplitNonConvexPoly(Pl);
  111.         CleanUpPolygonList(&PlSplit);       /* Zero length edges etc. */
  112.         if (PlSplit != NULL) {     /* Something is wrong here, ignore. */
  113.         if (Pl == PObj -> U.Pl) PObj -> U.Pl = PlSplit;    /* First. */
  114.         else PlPrev -> Pnext = PlSplit;
  115.  
  116.         PlTemp = PlSplit;
  117.         while (PlTemp -> Pnext != NULL) PlTemp = PlTemp -> Pnext;
  118.         PlTemp -> Pnext = Pl -> Pnext;
  119.         PlPrev = PlTemp;
  120.  
  121.         Pl -> Pnext = NULL;
  122.         MyFree((char *) Pl, POLYGON_TYPE);    /* Free old polygon. */
  123.         Pl = PlPrev -> Pnext;
  124.         }
  125.         else {
  126.         if (Pl == PObj -> U.Pl) PObj -> U.Pl = Pl -> Pnext;
  127.         PlPrev = Pl -> Pnext;
  128.         Pl -> Pnext = NULL;
  129.         MyFree((char *) Pl, POLYGON_TYPE);    /* Free old polygon. */
  130.         Pl = PlPrev;
  131.         }
  132.     }
  133.     else {
  134.         PlPrev = Pl;
  135.         Pl = Pl -> Pnext;
  136.     }
  137.     }
  138. }
  139.  
  140. /*****************************************************************************
  141. *   Routine to test if the given polygon is convex or not.             *
  142. * Algorithm: The polygon is convex iff the normals generated from cross      *
  143. * products of two consecutive edges points to the same direction. he same    *
  144. * direction is tested by dot product with the polygon plane normal which     *
  145. * should produce consistent sign for all normals, in order the polygon to    *
  146. * be convex.                                     *
  147. *   The routine returns TRUE iff the polygon is convex. In addition the      *
  148. * polygon CONVEX tag (see Irit.h) is also updated.                 *
  149. *   Note that if the polygon IS marked as convex, nothing is tested!         *
  150. *   Also convex polygons are tested so that the vertices are ordered such    *
  151. * that cross product of two adjacent edges will be in the normal direction.  *
  152. *****************************************************************************/
  153. int ConvexPolygon(PolygonStruct *Pl)
  154. {
  155.     int FirstTime = TRUE;
  156.     RealType NormalSign;
  157.     VectorType V1, V2, PolyNormal, Normal;
  158.     VertexStruct *V = Pl -> V, *VNext;
  159.  
  160.     if (IS_CONVEX_POLY(Pl)) return TRUE;       /* Nothing to do around here. */
  161.  
  162.     /* Copy only A, B, C from Ax+By+Cz+D = 0: */
  163.     PT_COPY(PolyNormal, Pl -> Plane);
  164.  
  165.     do {
  166.     VNext = V -> Pnext;
  167.  
  168.     PT_SUB(V1, VNext -> Pt, V -> Pt);
  169.     PT_SUB(V2, VNext -> Pnext -> Pt, VNext -> Pt);
  170.     VecCrossProd(Normal, V1, V2);
  171.     if (PT_LENGTH(Normal) < EPSILON) {
  172.         V = VNext;
  173.         continue;                    /* Skip colinear points. */
  174.     }
  175.  
  176.     if (FirstTime) {
  177.         FirstTime = FALSE;
  178.         NormalSign = DOT_PROD(Normal, PolyNormal);
  179.     }
  180.     else if (NormalSign * DOT_PROD(Normal, PolyNormal) < 0.0) {
  181.         RST_CONVEX_POLY(Pl);
  182.         return FALSE;          /* Different signs --> not convex. */
  183.     }
  184.  
  185.     V = VNext;
  186.     }
  187.     while (V != Pl -> V);
  188.  
  189.     SET_CONVEX_POLY(Pl);
  190.  
  191.     if (NormalSign < 0.0) ReverseVrtxList(Pl);
  192.  
  193.     return TRUE;    /* All signs are the same --> the polygon is convex. */
  194. }
  195.  
  196. /*****************************************************************************
  197. *   Routine to split non convex polygon to list of convex ones.             *
  198. * This algorithm is much simpler than the one used before:             *
  199. * 1. Remove a polygon from GlblList. If non exists stop.             *
  200. * 2. Search for non convex corner. If not found stop - polygon is convex.    *
  201. *    Otherwise let the non convex polygon found be V(i).             *
  202. * 3. Fire a ray from V(i) in the opposite direction to V(i-1). Find the      *
  203. *    closest intersection of the ray with polygon boundary P.             *
  204. * 4. Split the polygon into two at V(i)-P edge and push the two new polygons *
  205. *    on the GlblList.                                 *
  206. * 5. Goto 1.                                     *
  207. *****************************************************************************/
  208. PolygonStruct *SplitNonConvexPoly(PolygonStruct *Pl)
  209. {
  210.     int IsConvex;
  211.     PolygonStruct *GlblList, *GlblSplitPl = NULL, *Pl1, *Pl2;
  212.     VectorType V1, V2, PolyNormal, Normal;
  213.     VertexStruct *V, *VNext;
  214.  
  215.     TestConvexityDir(Pl);
  216.  
  217.     GlblList = AllocPolygon(0, 0, CopyVList(Pl -> V), NULL);
  218.     PLANE_COPY(GlblList -> Plane, Pl -> Plane);
  219.  
  220.     /* Copy only A, B, C from Ax+By+Cz+D = 0 plane equation: */
  221.     PT_COPY(PolyNormal, Pl -> Plane);
  222.  
  223.     while (GlblList != NULL) {
  224.     Pl = GlblList;
  225.     GlblList = GlblList -> Pnext;
  226.     Pl -> Pnext = NULL;
  227.  
  228.     IsConvex = TRUE;
  229.     V = Pl -> V;
  230.     do {
  231.         VNext = V -> Pnext;
  232.  
  233.         PT_SUB(V1, VNext -> Pt, V -> Pt);
  234.         PT_SUB(V2, VNext -> Pnext -> Pt, VNext -> Pt);
  235.         VecCrossProd(Normal, V1, V2);
  236.         if (PT_LENGTH(Normal) < EPSILON) {
  237.         V = VNext;
  238.         continue;                /* Skip colinear points. */
  239.         }
  240.  
  241.         if (DOT_PROD(Normal, PolyNormal) < 0.0) {
  242.         SplitPolyIntoTwo(Pl, V, &Pl1, &Pl2);
  243.  
  244.         Pl -> V = NULL;          /* Dont free vertices - used in Pl1/2. */
  245.         Pl -> Pnext = NULL;
  246.         MyFree((char *) Pl, POLYGON_TYPE);
  247.  
  248.         Pl1 -> Pnext = GlblList;   /* Push the polygons on GlblList. */
  249.         GlblList = Pl1;
  250.         Pl2 -> Pnext = GlblList;
  251.         GlblList = Pl2;
  252.  
  253.         IsConvex = FALSE;
  254.         break;
  255.         }
  256.  
  257.         V = VNext;
  258.     }
  259.     while (V != Pl -> V);
  260.  
  261.     if (IsConvex) {
  262.         SET_CONVEX_POLY(Pl);
  263.         Pl -> Pnext = GlblSplitPl;
  264.         GlblSplitPl = Pl;
  265.     }
  266.     }
  267.  
  268.     return GlblSplitPl;
  269. }
  270.  
  271. /*****************************************************************************
  272. *   Split polygon at the vertex specified by V -> Pnext, given V, into two,  *
  273. * by firing a ray from V -> Pnext in the opposite direction to V and finding *
  274. * closest intersection with the polygon P. (V -> Pnext, P) is the edge to    *
  275. * split the polygon at.                                 *
  276. *****************************************************************************/
  277. static void SplitPolyIntoTwo(PolygonStruct *Pl, VertexStruct *V,
  278.                 PolygonStruct **Pl1, PolygonStruct **Pl2)
  279. {
  280.     PointType Vl, PInter;
  281.     VertexStruct *VInter, *VNew1, *VNew2;
  282.  
  283.     PT_SUB(Vl, V -> Pnext -> Pt, V -> Pt);
  284.     VInter = FindRayPolyInter(Pl, V -> Pnext, Vl, PInter);
  285.     V = V -> Pnext;
  286.  
  287.     /* Make the two polygon vertices lists. */
  288.     VNew1 = AllocVertex(V -> Count, V -> Tags, NULL, V -> Pnext);
  289.     PT_COPY(VNew1 -> Pt, V -> Pt);
  290.     SET_INTERNAL_EDGE(V);
  291.     if (PT_EQ(VInter -> Pt, PInter)) {
  292.     /* Intersection points is close to VInter point. */
  293.     VNew2 = AllocVertex(VInter -> Count, VInter -> Tags, NULL,
  294.                             VInter -> Pnext);
  295.     PT_COPY(VNew2 -> Pt, VInter -> Pt);
  296.     VInter -> Pnext = VNew1;
  297.     SET_INTERNAL_EDGE(VInter);
  298.     V -> Pnext = VNew2;
  299.     }
  300.     else
  301.     if (PT_EQ(VInter -> Pnext -> Pt, PInter)) {
  302.     /* Intersection points is close to VInter -> Pnext point. */
  303.     VNew2 = AllocVertex(VInter -> Pnext -> Count, VInter -> Pnext -> Tags,
  304.                     NULL, VInter -> Pnext -> Pnext);
  305.     PT_COPY(VNew2 -> Pt, VInter -> Pnext -> Pt);
  306.     VInter -> Pnext -> Pnext = VNew1;
  307.     SET_INTERNAL_EDGE(VInter -> Pnext);
  308.     V -> Pnext = VNew2;
  309.     }
  310.     else {
  311.     /* PInter is in the middle of (VInter, VInter -> Pnext) edge: */
  312.     VNew2 = AllocVertex(VInter -> Count, VInter -> Tags, NULL,
  313.                             VInter -> Pnext);
  314.     PT_COPY(VNew2 -> Pt, PInter);
  315.     VInter -> Pnext = AllocVertex(0, 0, NULL, VNew1);
  316.     PT_COPY(VInter -> Pnext -> Pt, PInter);
  317.     SET_INTERNAL_EDGE(VInter -> Pnext);
  318.     V -> Pnext = VNew2;
  319.     }
  320.  
  321.     *Pl1 = AllocPolygon(0, 0, VNew1, NULL);
  322.     PLANE_COPY((*Pl1) -> Plane, Pl -> Plane);
  323.     *Pl2 = AllocPolygon(0, 0, VNew2, NULL);
  324.     PLANE_COPY((*Pl2) -> Plane, Pl -> Plane);
  325. }
  326.  
  327. /*****************************************************************************
  328. *   Find where a ray first intersect a given polygon. The ray starts at one  *
  329. * of the polygon vertices and so distance less than EPSILON is ignored.      *
  330. *   Returns the vertex V in which (V, V -> Pnext) has the closest         *
  331. * intersection with the ray PRay, DRay at Inter.                 *
  332. *****************************************************************************/
  333. static VertexStruct *FindRayPolyInter(PolygonStruct *Pl,
  334.             VertexStruct *VRay, PointType RayDir, PointType PInter)
  335. {
  336.     RealType t1, t2, MinT = INFINITY;
  337.     PointType Vl, Ptemp1, Ptemp2, PRay;
  338.     VertexStruct *V = Pl -> V, *VNext, *VInter = NULL;
  339.  
  340.     PT_COPY(PRay, VRay -> Pt);
  341.  
  342.     do {
  343.     VNext = V -> Pnext;
  344.     if (V != VRay && VNext != VRay) {
  345.         PT_SUB(Vl, V -> Pnext -> Pt, V -> Pt);
  346.         if (CGDistPointLine(PRay, V -> Pt, Vl) > EPSILON) {
  347.         /* Only if the point the ray is shoot from is not on line: */
  348.         CG2PointsFromLineLine(PRay, RayDir, V -> Pt, Vl, Ptemp1, &t1,
  349.                                  Ptemp2, &t2);
  350.         if (CGDistPointPoint(Ptemp1, Ptemp2) < EPSILON * 10.0 &&
  351.             t1 > EPSILON && t1 < MinT &&
  352.             (t2 <= 1.0 || APX_EQ(t2, 1.0)) &&
  353.             (t2 >= 0.0 || APX_EQ(t2, 0.0))) {
  354.             PT_COPY(PInter, Ptemp2);
  355.             VInter = V;
  356.             MinT = t1;
  357.         }
  358.         }
  359.     }
  360.     V = VNext;
  361.     }
  362.     while (V != Pl -> V);
  363.  
  364.     if (VInter == NULL)
  365.     FatalError("FindRayPolyInter: No Ray intersection in Convex decomp.");
  366.  
  367.     return VInter;
  368. }
  369.  
  370. /*****************************************************************************
  371. *   Test convexity direction - a cross product of two edges of a convex      *
  372. * corner of the polygon should point in the normal direction. if this is not *
  373. * the case - the polygon vertices are reveresed.                 *
  374. *****************************************************************************/
  375. static void TestConvexityDir(PolygonStruct *Pl)
  376. {
  377.     int Coord = 0;
  378.     VectorType V1, V2, Normal;
  379.     VertexStruct *V, *VExtrem;
  380.  
  381.     /* Find the minimum component direction of the normal and used that axes */
  382.     /* to find an extremum point on the polygon - that extrmum point must be */
  383.     /* a convex corner so we can find the normal direction of convex corner. */
  384.     if (ABS(Pl -> Plane[1]) < ABS(Pl -> Plane[Coord])) Coord = 1;
  385.     if (ABS(Pl -> Plane[2]) < ABS(Pl -> Plane[Coord])) Coord = 2;
  386.     V = VExtrem = Pl -> V;
  387.     do {
  388.     if (V -> Pt[Coord] > VExtrem -> Pt[Coord]) VExtrem = V;
  389.     V = V -> Pnext;
  390.     }
  391.     while (V != Pl -> V);
  392.  
  393.     /* Make sure next vertex is not at the extremum value: */
  394.     while (APX_EQ(VExtrem -> Pt[Coord], VExtrem -> Pnext -> Pt[Coord]))
  395.                         VExtrem = VExtrem -> Pnext;
  396.  
  397.     /* O.K. V form a convex corner - evaluate its two edges cross product:   */
  398.     for (V = Pl -> V; V -> Pnext != VExtrem; V = V -> Pnext);  /* Find Prev. */
  399.     PT_SUB(V1, VExtrem -> Pt, V -> Pt);
  400.     PT_SUB(V2, VExtrem -> Pnext -> Pt, VExtrem -> Pt);
  401.     VecCrossProd(Normal, V1, V2);
  402.  
  403.     if (DOT_PROD(Normal, Pl -> Plane) < 0.0) ReverseVrtxList(Pl);
  404. }
  405.  
  406. /*****************************************************************************
  407. *   Reverse the vertex list of a given polygon. This is used mainly to       *
  408. * reverse polygons such that cross product of consecutive edges which form   *
  409. * a convex corner will point in the polygon normal direction.             *
  410. *****************************************************************************/
  411. void ReverseVrtxList(PolygonStruct *Pl)
  412. {
  413.     ByteType Tags, Count;
  414.     VertexStruct *V = Pl -> V, *VNext = V -> Pnext, *VNextNext;
  415.  
  416.     do {
  417.     VNextNext = VNext -> Pnext;
  418.     VNext -> Pnext = V;                 /* Reverse the pointer! */
  419.  
  420.     V = VNext;               /* Advance all 3 pointers by one. */
  421.     VNext = VNextNext;
  422.     VNextNext = VNextNext -> Pnext;
  423.     }
  424.     while (V != Pl -> V);
  425.  
  426.     V = Pl -> V;      /* Move the Tags/Count by one - to the right edge. */
  427.     Tags = V -> Tags;
  428.     Count = V -> Count;
  429.     do {
  430.     if (V -> Pnext == Pl -> V) {
  431.         V -> Tags = Tags;
  432.         V -> Count = Count;
  433.     }
  434.     else {
  435.         V -> Tags = V -> Pnext -> Tags;
  436.         V -> Count = V -> Pnext -> Count;
  437.     }
  438.  
  439.     V = V -> Pnext;
  440.     }
  441.     while (V != Pl -> V);
  442. }
  443.  
  444. #ifdef DEBUG2
  445.  
  446. /*****************************************************************************
  447. *   Print the content of the given vertex list, to standard output.         *
  448. *****************************************************************************/
  449. static void PrintVrtxList(VertexStruct *V)
  450. {
  451.     VertexStruct *VFirst = V;
  452.  
  453.     if (V == NULL) return;
  454.  
  455.     printf("VERTEX LIST:\n");
  456.     do {
  457.     printf("@%p:  %12lf %12lf %12lf, Tags = %02x\n",
  458.         V, V -> Pt[0], V -> Pt[1], V -> Pt[2], V -> Tags);
  459.     V = V -> Pnext;
  460.     }
  461.     while (V != NULL && V != VFirst);
  462.  
  463.     if (V == NULL) printf("Loop is not closed (NULL terminated)\n");
  464. }
  465.  
  466. #endif /* DEBUG2 */
  467.