home *** CD-ROM | disk | FTP | other *** search
/ AP Professional Graphics CD-ROM Library / AP Professional Graphics CD-ROM Library.iso / pc / unix / appendix / gemsiii / planeset.c < prev    next >
Encoding:
C/C++ Source or Header  |  1992-10-05  |  7.5 KB  |  261 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. /* declare local routines */
  79. static void computePlaneEq() ;
  80. static Node *insertPlaneEq() ;
  81. static int comparePlaneEqs() ;
  82.  
  83. /*
  84.     coplanarSets--takes as input a linked list of polygons 'plist',
  85.     and two tolerances 'linearEps' and 'angularEps' and returns a
  86.     binary tree of sets of coplanar polygons.
  87.     The tolerances are used to deal with floating-point arithmetic
  88.     and numerical errors when comparing plane equations; two plane
  89.     equations are considered equal if the angle between their
  90.     respective normals is less than or equal to angularEps (in
  91.     degrees) and the distance between the two planes is less than
  92.     or equal to linearEps.
  93. */
  94. Node *coplanarSets(plist, linearEps, angularEps)
  95. Polygon *plist;
  96. float linearEps, angularEps;
  97. {
  98.     Node *tree;
  99.     Plane plane;
  100.     void computePlaneEq();
  101.     Node *pset, *insertPlaneEq();
  102.     Polygon *polygon, *nextPolygon;
  103.  
  104.     /* initialize the tolerances used by comparePlaneEqs() */
  105.     linEps = linearEps;
  106.     cosEps = cos(angularEps * M_PI / 180.0);
  107.  
  108.     /* place each input polygon in the appropriate set
  109.        of coplanar polygons
  110.     */
  111.     tree = NULL;
  112.     polygon = plist;
  113.     while(polygon != NULL) {
  114.         nextPolygon = polygon->next;
  115.         /* first, compute the plane equation of the polygon */
  116.         computePlaneEq(polygon, plane);
  117.         /* then, find the set of coplanar polygons with this plane
  118.            equation or create a new, empty one if none exist
  119.         */
  120.         tree = insertPlaneEq(tree, plane, &pset);
  121.         /* finally, add the polygon to the selected set of
  122.            coplanar polygons
  123.         */
  124.         polygon->next = pset->plist;
  125.         pset->plist = polygon;
  126.         /* go to the next polygon in the input list. Note that the
  127.            'next' field in the polygon structure is reused to
  128.            assemble lists of coplanar polygons; thus the necessity
  129.            for 'nextPolygon'
  130.         */
  131.         polygon = nextPolygon;
  132.     }
  133.  
  134.     return(tree);
  135. }
  136.  
  137. /*
  138.     computePlaneEq--takes as input a pointer 'polygon' to an
  139.     arbitrary 3D polygon and returns in 'plane' the normalized
  140.     (unit normal) plane equation of the polygon.
  141.     Newell's method (see "Newell's Method for Computing the Plane
  142.     Equation of a Polygon" in this volume) is used for the
  143.     computation.
  144. */
  145. static void computePlaneEq(polygon, plane)
  146. Polygon *polygon;
  147. Plane plane;
  148. {
  149.     int i;
  150.     Point refpt;
  151.     Normal normal;
  152.     float *u, *v, len;
  153.  
  154.     /* first, compute the normal of the input polygon */
  155.     VZERO(normal);
  156.     VZERO(refpt);
  157.     for(i = 0; i < polygon->numVertices; i++) {
  158.         u = polygon->vertices[i];
  159.         v = polygon->vertices[(i + 1) % polygon->numVertices];
  160.         normal[X] += (u[Y] - v[Y]) * (u[Z] + v[Z]);
  161.         normal[Y] += (u[Z] - v[Z]) * (u[X] + v[X]);
  162.         normal[Z] += (u[X] - v[X]) * (u[Y] + v[Y]);
  163.         VINCR(refpt, u);
  164.     }
  165.     /* then, compute the normalized plane equation using the
  166.        arithmetic average of the vertices of the input polygon to
  167.        determine its last coefficient. Note that, for efficiency,
  168.        'refpt' stores the sum of the vertices rather than the
  169.        actual average; the division by 'polygon->numVertices' is
  170.        carried out together with the normalization when computing
  171.        'plane[D]'.
  172.     */
  173.     len = VNORM(normal);
  174.     plane[X] = normal[X] / len;
  175.     plane[Y] = normal[Y] / len;
  176.     plane[Z] = normal[Z] / len;
  177.     len *= polygon->numVertices;
  178.     plane[D] = -VDOT(refpt, normal) / len;
  179. }
  180.  
  181. /*
  182.     insertPlaneEq--takes as input a binary tree 'tree' of sets of
  183.     coplanar polygons, and a plane equation 'plane', and returns
  184.     a pointer 'pset' to a set of coplanar polygons with the given
  185.     plane equation. A new, empty set is created if none is found.
  186. */
  187.  
  188. static Node *insertPlaneEq(tree, plane, pset)
  189. Node *tree, **pset;
  190. Plane plane;
  191. {
  192.     int i, c, comparePlaneEqs();
  193.  
  194.     if(tree == NULL) {
  195.         /* the input plane is not in the tree.
  196.            Create a new set for it
  197.         */
  198.         tree = (Node *)malloc(sizeof(Node));
  199.         for(i = 0; i < 4; i++)
  200.             tree->plane[i] = plane[i];
  201.         tree->plist = NULL; /* no polygons in this set for now */
  202.         tree->left = tree->right = NULL;
  203.         *pset = tree;
  204.     } else {
  205.         /* compare the input plane equation with that
  206.            of the visited node
  207.         */
  208.         c = comparePlaneEqs(plane, tree->plane);
  209.         if(c < 0)
  210.             tree->left = insertPlaneEq(tree->left, plane, pset);
  211.         else if(c > 0)
  212.             tree->right = insertPlaneEq(tree->right, plane, pset);
  213.         else
  214.             /* the input plane is approximately equal to that
  215.                of this node
  216.             */
  217.             *pset = tree;
  218.     }
  219.  
  220.     return(tree);
  221. }
  222.  
  223. /*
  224.     comparePlaneEqs--compares two plane equations, 'p1' and 'p2',
  225.     and returns an integer less than, equal to, or greater than
  226.     zero, depending on whether 'p1' is less than, equal to, or
  227.     greater than 'p2'.  The two global variables, 'linEps' and
  228.     'cosEps' are tolerances used to account for numerical errors.
  229. */
  230. static int comparePlaneEqs(p1, p2)
  231. Plane p1, p2;
  232. {
  233.     int ret;
  234.     float cosTheta;
  235.  
  236.     /* compute the cosine of the angle between the normals
  237.        of the two planes
  238.     */
  239.     cosTheta = VDOT(p1, p2);
  240.  
  241.     if(cosTheta < 0.0)
  242.         ret = -1;
  243.     else if(cosTheta < cosEps)
  244.         ret = 1;
  245.     else {
  246.         /* the two planes are parallel and oriented in
  247.            the same direction
  248.         */
  249.         if(p1[3] + linEps < p2[3])
  250.             ret = -1;
  251.         else if(p2[3] + linEps < p1[3])
  252.             ret = 1;
  253.         else
  254.             /* the two planes are equal */
  255.             ret = 0;
  256.     }
  257.  
  258.     return ret;
  259. }
  260.  
  261.