home *** CD-ROM | disk | FTP | other *** search
/ Photo CD Demo 1 / Demo.bin / gems / gemsiii / planests.c < prev    next >
C/C++ Source or Header  |  1992-03-16  |  8KB  |  256 lines
  1. /*
  2. GROUPING NEARLY COPLANAR POLYGONS INTO COPLANAR SETS
  3.  
  4.  
  5. David Salesin
  6. Filippo Tampieri
  7.  
  8. Cornell University
  9. */
  10.  
  11.  
  12.  
  13.  
  14.  
  15.  
  16.  
  17.  
  18.  
  19. /*
  20.     This code partitions a given set of arbitrary 3D polygons into
  21.     subsets of coplanar polygons.
  22.     The input polygons are organized in a linked list and the
  23.     resulting subsets of coplanar polygons are returned in the form
  24.     of a binary tree; each node of the binary tree contains a
  25.     different subset (once again stored as a linked list) and its
  26.     plane equation.  An inorder traversal of the binary tree will
  27.     return the sets of coplanar polygons sorted by plane equations
  28.     according to the total ordering imposed by the relation
  29.     implemented by the routine 'comparePlaneEqs'.
  30. */
  31.  
  32. #include <stdio.h>
  33. #include <math.h>
  34.  
  35. #define X   0
  36. #define Y   1
  37. #define Z   2
  38. #define D   3
  39.  
  40. #define VZERO(v)    (v[X] = v[Y] = v[Z] = 0.0)
  41. #define VNORM(v)    (sqrt(v[X] * v[X] + v[Y] * v[Y] + v[Z] * v[Z]))
  42. #define VDOT(u, v)  (u[0] * v[0] + u[1] * v[1] + u[2] * v[2])
  43. #define VINCR(u, v) (u[X] += v[X], u[Y] += v[Y], u[Z] += v[Z])
  44.  
  45. typedef float Vector[3];
  46. typedef Vector Point;
  47. typedef Vector Normal;
  48. typedef float Plane[4];
  49.  
  50. /*
  51.     Polygon--a polygon is stored as an array 'vertices' of size
  52.     'numVertices'.  Pointer 'next' is used to implement sets of
  53.     polygons through linked lists.
  54. */
  55. typedef struct polygon {
  56.     Point *vertices;
  57.     int numVertices;
  58.     struct polygon *next;
  59. } Polygon;
  60.  
  61. /*
  62.     Node--each node stores a set of coplanar polygons. The set is
  63.     implemented as a linked list pointed to by 'plist', and the
  64.     plane equation of the set is stored in 'plane'. Pointers 'left'
  65.     and 'right' point to the subtrees containing sets of coplanar
  66.     polygons with plane equations respectively less than and
  67.     greater than that stored in 'plane'.
  68. */
  69. typedef struct node {
  70.         Plane plane;
  71.         Polygon *plist;
  72.         struct node *left, *right;
  73. } Node;
  74.  
  75. static float linEps;  /* tolerances used by 'comparePlaneEq' to */
  76. static float cosEps;  /* account for numerical errors           */
  77.  
  78. /*
  79.     coplanarSets--takes as input a linked list of polygons 'plist',
  80.     and two tolerances 'linearEps' and 'angularEps' and returns a
  81.     binary tree of sets of coplanar polygons.
  82.     The tolerances are used to deal with floating-point arithmetic
  83.     and numerical errors when comparing plane equations; two plane
  84.     equations are considered equal if the angle between their
  85.     respective normals is less than or equal to angularEps (in
  86.     degrees) and the distance between the two planes is less than
  87.     or equal to linearEps.
  88. */
  89. Node *coplanarSets(plist, linearEps, angularEps)
  90. Polygon *plist;
  91. float linearEps, angularEps;
  92. {
  93.     Node *tree;
  94.     Plane plane;
  95.     void computePlaneEq();
  96.     Node *pset, *insertPlaneEq();
  97.     Polygon *polygon, *nextPolygon;
  98.  
  99.     /* initialize the tolerances used by comparePlaneEqs() */
  100.     linEps = linearEps;
  101.     cosEps = cos(angularEps * M_PI / 180.0);
  102.  
  103.     /* place each input polygon in the appropriate set
  104.        of coplanar polygons
  105.     */
  106.     tree = NULL;
  107.     polygon = plist;
  108.     while(polygon != NULL) {
  109.         nextPolygon = polygon->next;
  110.         /* first, compute the plane equation of the polygon */
  111.         computePlaneEq(polygon, plane);
  112.         /* then, find the set of coplanar polygons with this plane
  113.            equation or create a new, empty one if none exist
  114.         */
  115.         tree = insertPlaneEq(tree, plane, &pset);
  116.         /* finally, add the polygon to the selected set of
  117.            coplanar polygons
  118.         */
  119.         polygon->next = pset->plist;
  120.         pset->plist = polygon;
  121.         /* go to the next polygon in the input list. Note that the
  122.            'next' field in the polygon structure is reused to
  123.            assemble lists of coplanar polygons; thus the necessity
  124.            for 'nextPolygon'
  125.         */
  126.         polygon = nextPolygon;
  127.     }
  128.  
  129.     return(tree);
  130. }
  131.  
  132. /*
  133.     computePlaneEq--takes as input a pointer 'polygon' to an
  134.     arbitrary 3D polygon and returns in 'plane' the normalized
  135.     (unit normal) plane equation of the polygon.
  136.     Newell's method (see "Newell's Method for Computing the Plane
  137.     Equation of a Polygon" in this volume) is used for the
  138.     computation.
  139. */
  140. static void computePlaneEq(polygon, plane)
  141. Polygon *polygon;
  142. Plane plane;
  143. {
  144.     int i;
  145.     Point refpt;
  146.     Normal normal;
  147.     float *u, *v, len;
  148.  
  149.     /* first, compute the normal of the input polygon */
  150.     VZERO(normal);
  151.     VZERO(refpt);
  152.     for(i = 0; i < polygon->numVertices; i++) {
  153.         u = polygon->vertices[i];
  154.         v = polygon->vertices[(i + 1) % polygon->numVertices];
  155.         normal[X] += (u[Y] - v[Y]) * (u[Z] + v[Z]);
  156.         normal[Y] += (u[Z] - v[Z]) * (u[X] + v[X]);
  157.         normal[Z] += (u[X] - v[X]) * (u[Y] + v[Y]);
  158.         VINCR(refpt, u);
  159.     }
  160.     /* then, compute the normalized plane equation using the
  161.        arithmetic average of the vertices of the input polygon to
  162.        determine its last coefficient. Note that, for efficiency,
  163.        'refpt' stores the sum of the vertices rather than the
  164.        actual average; the division by 'polygon->numVertices' is
  165.        carried out together with the normalization when computing
  166.        'plane[D]'.
  167.     */
  168.     len = VNORM(normal);
  169.     plane[X] = normal[X] / len;
  170.     plane[Y] = normal[Y] / len;
  171.     plane[Z] = normal[Z] / len;
  172.     len *= polygon->numVertices;
  173.     plane[D] = -VDOT(refpt, normal) / len;
  174. }
  175.  
  176. /*
  177.     insertPlaneEq--takes as input a binary tree 'tree' of sets of
  178.     coplanar polygons, and a plane equation 'plane', and returns
  179.     a pointer 'pset' to a set of coplanar polygons with the given
  180.     plane equation. A new, empty set is created if none is found.
  181. */
  182.  
  183. static Node *insertPlaneEq(tree, plane, pset)
  184. Node *tree, **pset;
  185. Plane plane;
  186. {
  187.     int i, c, comparePlaneEqs();
  188.  
  189.     if(tree == NULL) {
  190.         /* the input plane is not in the tree.
  191.            Create a new set for it
  192.         */
  193.         tree = (Node *)malloc(sizeof(Node));
  194.         for(i = 0; i < 4; i++)
  195.             tree->plane[i] = plane[i];
  196.         tree->plist = NULL; /* no polygons in this set for now */
  197.         tree->left = tree->right = NULL;
  198.         *pset = tree;
  199.     } else {
  200.         /* compare the input plane equation with that
  201.            of the visited node
  202.         */
  203.         c = comparePlaneEqs(plane, tree->plane);
  204.         if(c < 0)
  205.             tree->left = insertPlaneEq(tree->left, plane, pset);
  206.         else if(c > 0)
  207.             tree->right = insertPlaneEq(tree->right, plane, pset);
  208.         else
  209.             /* the input plane is approximately equal to that
  210.                of this node
  211.             */
  212.             *pset = tree;
  213.     }
  214.  
  215.     return(tree);
  216. }
  217.  
  218. /*
  219.     comparePlaneEqs--compares two plane equations, 'p1' and 'p2',
  220.     and returns an integer less than, equal to, or greater than
  221.     zero, depending on whether 'p1' is less than, equal to, or
  222.     greater than 'p2'.  The two global variables, 'linEps' and
  223.     'cosEps' are tolerances used to account for numerical errors.
  224. */
  225. static int comparePlaneEqs(p1, p2)
  226. Plane p1, p2;
  227. {
  228.     int ret;
  229.     float cosTheta;
  230.  
  231.     /* compute the cosine of the angle between the normals
  232.        of the two planes
  233.     */
  234.     cosTheta = VDOT(p1, p2);
  235.  
  236.     if(cosTheta < 0.0)
  237.         ret = -1;
  238.     else if(cosTheta < cosEps)
  239.         ret = 1;
  240.     else {
  241.         /* the two planes are parallel and oriented in
  242.            the same direction
  243.         */
  244.         if(p1[3] + linEps < p2[3])
  245.             ret = -1;
  246.         else if(p2[3] + linEps < p1[3])
  247.             ret = 1;
  248.         else
  249.             /* the two planes are equal */
  250.             ret = 0;
  251.     }
  252.  
  253.     return ret;
  254. }
  255.  
  256.