/*
GROUPING NEARLY COPLANAR POLYGONS INTO COPLANAR SETS
David Salesin
Filippo Tampieri
Cornell University
*/
/*
This code partitions a given set of arbitrary 3D polygons into
subsets of coplanar polygons.
The input polygons are organized in a linked list and the
resulting subsets of coplanar polygons are returned in the form
of a binary tree; each node of the binary tree contains a
different subset (once again stored as a linked list) and its
plane equation. An inorder traversal of the binary tree will
return the sets of coplanar polygons sorted by plane equations
according to the total ordering imposed by the relation
implemented by the routine 'comparePlaneEqs'.
*/
#include
#include
#define X 0
#define Y 1
#define Z 2
#define D 3
#define VZERO(v) (v[X] = v[Y] = v[Z] = 0.0)
#define VNORM(v) (sqrt(v[X] * v[X] + v[Y] * v[Y] + v[Z] * v[Z]))
#define VDOT(u, v) (u[0] * v[0] + u[1] * v[1] + u[2] * v[2])
#define VINCR(u, v) (u[X] += v[X], u[Y] += v[Y], u[Z] += v[Z])
typedef float Vector[3];
typedef Vector Point;
typedef Vector Normal;
typedef float Plane[4];
/*
Polygon--a polygon is stored as an array 'vertices' of size
'numVertices'. Pointer 'next' is used to implement sets of
polygons through linked lists.
*/
typedef struct polygon {
Point *vertices;
int numVertices;
struct polygon *next;
} Polygon;
/*
Node--each node stores a set of coplanar polygons. The set is
implemented as a linked list pointed to by 'plist', and the
plane equation of the set is stored in 'plane'. Pointers 'left'
and 'right' point to the subtrees containing sets of coplanar
polygons with plane equations respectively less than and
greater than that stored in 'plane'.
*/
typedef struct node {
Plane plane;
Polygon *plist;
struct node *left, *right;
} Node;
static float linEps; /* tolerances used by 'comparePlaneEq' to */
static float cosEps; /* account for numerical errors */
/* declare local routines */
static void computePlaneEq() ;
static Node *insertPlaneEq() ;
static int comparePlaneEqs() ;
/*
coplanarSets--takes as input a linked list of polygons 'plist',
and two tolerances 'linearEps' and 'angularEps' and returns a
binary tree of sets of coplanar polygons.
The tolerances are used to deal with floating-point arithmetic
and numerical errors when comparing plane equations; two plane
equations are considered equal if the angle between their
respective normals is less than or equal to angularEps (in
degrees) and the distance between the two planes is less than
or equal to linearEps.
*/
Node *coplanarSets(plist, linearEps, angularEps)
Polygon *plist;
float linearEps, angularEps;
{
Node *tree;
Plane plane;
void computePlaneEq();
Node *pset, *insertPlaneEq();
Polygon *polygon, *nextPolygon;
/* initialize the tolerances used by comparePlaneEqs() */
linEps = linearEps;
cosEps = cos(angularEps * M_PI / 180.0);
/* place each input polygon in the appropriate set
of coplanar polygons
*/
tree = NULL;
polygon = plist;
while(polygon != NULL) {
nextPolygon = polygon->next;
/* first, compute the plane equation of the polygon */
computePlaneEq(polygon, plane);
/* then, find the set of coplanar polygons with this plane
equation or create a new, empty one if none exist
*/
tree = insertPlaneEq(tree, plane, &pset);
/* finally, add the polygon to the selected set of
coplanar polygons
*/
polygon->next = pset->plist;
pset->plist = polygon;
/* go to the next polygon in the input list. Note that the
'next' field in the polygon structure is reused to
assemble lists of coplanar polygons; thus the necessity
for 'nextPolygon'
*/
polygon = nextPolygon;
}
return(tree);
}
/*
computePlaneEq--takes as input a pointer 'polygon' to an
arbitrary 3D polygon and returns in 'plane' the normalized
(unit normal) plane equation of the polygon.
Newell's method (see "Newell's Method for Computing the Plane
Equation of a Polygon" in this volume) is used for the
computation.
*/
static void computePlaneEq(polygon, plane)
Polygon *polygon;
Plane plane;
{
int i;
Point refpt;
Normal normal;
float *u, *v, len;
/* first, compute the normal of the input polygon */
VZERO(normal);
VZERO(refpt);
for(i = 0; i < polygon->numVertices; i++) {
u = polygon->vertices[i];
v = polygon->vertices[(i + 1) % polygon->numVertices];
normal[X] += (u[Y] - v[Y]) * (u[Z] + v[Z]);
normal[Y] += (u[Z] - v[Z]) * (u[X] + v[X]);
normal[Z] += (u[X] - v[X]) * (u[Y] + v[Y]);
VINCR(refpt, u);
}
/* then, compute the normalized plane equation using the
arithmetic average of the vertices of the input polygon to
determine its last coefficient. Note that, for efficiency,
'refpt' stores the sum of the vertices rather than the
actual average; the division by 'polygon->numVertices' is
carried out together with the normalization when computing
'plane[D]'.
*/
len = VNORM(normal);
plane[X] = normal[X] / len;
plane[Y] = normal[Y] / len;
plane[Z] = normal[Z] / len;
len *= polygon->numVertices;
plane[D] = -VDOT(refpt, normal) / len;
}
/*
insertPlaneEq--takes as input a binary tree 'tree' of sets of
coplanar polygons, and a plane equation 'plane', and returns
a pointer 'pset' to a set of coplanar polygons with the given
plane equation. A new, empty set is created if none is found.
*/
static Node *insertPlaneEq(tree, plane, pset)
Node *tree, **pset;
Plane plane;
{
int i, c, comparePlaneEqs();
if(tree == NULL) {
/* the input plane is not in the tree.
Create a new set for it
*/
tree = (Node *)malloc(sizeof(Node));
for(i = 0; i < 4; i++)
tree->plane[i] = plane[i];
tree->plist = NULL; /* no polygons in this set for now */
tree->left = tree->right = NULL;
*pset = tree;
} else {
/* compare the input plane equation with that
of the visited node
*/
c = comparePlaneEqs(plane, tree->plane);
if(c < 0)
tree->left = insertPlaneEq(tree->left, plane, pset);
else if(c > 0)
tree->right = insertPlaneEq(tree->right, plane, pset);
else
/* the input plane is approximately equal to that
of this node
*/
*pset = tree;
}
return(tree);
}
/*
comparePlaneEqs--compares two plane equations, 'p1' and 'p2',
and returns an integer less than, equal to, or greater than
zero, depending on whether 'p1' is less than, equal to, or
greater than 'p2'. The two global variables, 'linEps' and
'cosEps' are tolerances used to account for numerical errors.
*/
static int comparePlaneEqs(p1, p2)
Plane p1, p2;
{
int ret;
float cosTheta;
/* compute the cosine of the angle between the normals
of the two planes
*/
cosTheta = VDOT(p1, p2);
if(cosTheta < 0.0)
ret = -1;
else if(cosTheta < cosEps)
ret = 1;
else {
/* the two planes are parallel and oriented in
the same direction
*/
if(p1[3] + linEps < p2[3])
ret = -1;
else if(p2[3] + linEps < p1[3])
ret = 1;
else
/* the two planes are equal */
ret = 0;
}
return ret;
}