Robust polygon normal calculation

Question

is there a good robust algorithm to calculate normal vector of a convex polygon (in 3D, of course)? For triangles, it is easy: one takes two of the triangle's edges and calculates the cross product:

vec3 u = point[0] - point[1], v = point[0] - point[2];
vec3 n = normalize(cross(u, v));

But this approach does not really extend to polygons very well. Some edges of the polygon can be nearly or "exactly" collinear (this will happen often in meshes where T-Junction removal took place), therefore it is necessary to choose a pair of edges, giving a "strong" normal (both edges are "long enough" and they hold "almost perpendicular" angle).

This approach will still not work for all polygons, though. Imagine a polygon in the shape of a disc. If it is very finely subdivided, all the edges will be very short and all of the consecutive edges will be almost collinear, regardless of the radius of the disc. At the same time, the normal is very well defined.

One solution could be to find the largest inscribed triangle and to calculate the normal of that. However, finding it will have complexity of O(n^2), which seems prohibitive.

A better solution could be to use SVD or Eigenvalue decomposition to calculate the normal, given all the polygon points, not just three or four.

Is there a standard algorithm for this? Anyone has a good way of doing it?

Answer 1

If you factorize the formula for a triangle, you'll get the following:

n ~ (p1 - p0) x (p2 - p0)
  = p0 x p1 + p1 x p2 + p2 x p0

You can generalize this formula for arbitrary polygons:

n ~ p0 x p1 + p1 x p2 + ... + pn x p0

So sum the cross product of consecutive edges. This is a robust algorithm and works for non-planar polygons.

If you can be sure that the polygon is planar, I would do the following (to save computation time):

Repeat k times
    Pick 3 random polygon vertices
    Calculate the normal of the according triangle
Choose the longest normal as the polygon's normal.

You might discard any normal that has a length <= epsilon.