/* NEWELL'S METHOD FOR COMPUTING THE PLANE EQUATION OF A POLYGON Filippo Tampieri Cornell University */ #include /* definition for the components of vectors and plane equations */ #define X 0 #define Y 1 #define Z 2 #define D 3 /* a few useful vector operations */ #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]) /* type definitions for vectors and plane equations */ typedef float Vector[3]; typedef Vector Point; typedef Vector Normal; typedef float Plane[4]; /* ** PlaneEquation--computes the plane equation of an arbitrary ** 3D polygon using Newell's method. ** ** Entry: ** verts - list of the vertices of the polygon ** nverts - number of vertices of the polygon ** Exit: ** plane - normalized (unit normal) plane equation */ PlaneEquation(verts, nverts, plane) Point *verts; int nverts; Plane plane; { int i; Point refpt; Normal normal; float *u, *v, len; /* compute the polygon normal and a reference point on the plane. Note that the actual reference point is refpt / nverts */ VZERO(normal); VZERO(refpt); for(i = 0; i < nverts; i++) { u = verts[i]; v = verts[(i + 1) % nverts]; 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); } /* normalize the polygon normal to obtain the first three coefficients of the plane equation */ len = VNORM(normal); plane[X] = normal[X] / len; plane[Y] = normal[Y] / len; plane[Z] = normal[Z] / len; /* compute the last coefficient of the plane equation */ len *= nverts; plane[D] = -VDOT(refpt, normal) / len; }