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?

Propane

Use Newell's algorithm for robustness: math.stackexchange.com/questions/2885839/… – Guillermoguilloche

khronos.org/opengl/wiki/Calculating_a_Surface_Normal – Guillermoguilloche

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.

Soubrette

Have you considered the disc example? I don't see it as very robust, since the cross products of the consecutive edges will be fairly small in there. But otherwise thanks, i guess this will work better most of the time. – Propane

The robustness comes with the summation. – Soubrette

Also, using epsilon guarantees the method not to be numerically robust. Consider having a mesh where there are narrow polygons. If the threshold is high, those polygons will not have a normal. If it is too low, the algorithm may accept a slightly non-planar edges or otherwise suboptimal edge pair, yielding imprecise normal even on polygons that have the normal well defined. – Propane

I'm not convinced about summation. Summing many small values does not work very well in limited floating point precision. – Propane

This is also known as Newell's Method (opengl.org/wiki/Calculating_a_Surface_Normal). – Propane

@theswine Are you sure this is Newell's Method, the formula at the link you provided does not look the same. – Shirtwaist

@LennyWhite khronos.org/opengl/wiki/… ?? It's a chapter in the page I sent. Are you seeing something else? – Propane

@theswine If we iterate through every two points of the polygon cur and next. The formula at the link is normal.x += (cur.y-next.y)*(cur.z+next.z); normal.y += (cur.z-next.z)*(cur.x+next.x); normal.z += (cur.x-next.x)*(cur.y+next.y);. While here the formula is normal += cross(cur, next). If you write the cross product in terms of the x,y,z components you don't get the same result. At least I didn't. – Shirtwaist

@Lenny: The individual summands are not equal, however, the sum is. If we take a look at the x-component, the expanded formula from Newell's method is: nx += cur.y * next.z - next.y * cur.z - next.y * next.z + cur.y * cur.z. The first two summands are from the cross product. The last two are additional. So every vertex adds y * z for itself and subtracts y * z for the next vertex. If you do this for all vertices, the additional terms cancel out and you get the sum of cross products. – Soubrette

This is also listed as "best fit normal" in paragraph 9.5.3 in 3D Math Primer for Graphics and Game Development" by Fletscher Dunn & Ian Parberry: books.google.ch/… – Extraordinary

Answer 2

start from an arbitrary vertex(lets call it vertex A), and move to the next vertex in the list(call it vertex B). calculate the perpendicular vector(call it vector P) to the AB vector. Then continue iterating in the vertex list to find the vertex that is perpendicularly the most distant from vector AB. So at each iteration take the dot product of the current element(take vertex B as the origin) with the vector P and take the one that has the greatest result in magnitude(take absolute value) and call it C. calculate the cross product of A B C vectors.

if the poly is convex you can stop iterating untill the perpendicular distances starts to get smaller in magnitude.

I came up with this idea, i do not know how efficient this method would be since I do not know any other algorithm to compare with.

Subcontinent

You cannot really calculate one perpendicular vector in 3-space. There are infinitely many of those, even if you limit yourself to unit length vectors. You can calculate distance from a point to a line passing through two points (AB). You can select such C that is the maximum distance from AB. This is very similar to choosing different ABCs and scoring them by the norm of their cross product. – Propane

In my solution, I assumed that the polygon is a flat surface, and all the points are coplanar. – Subcontinent

That sort of makes for an chicken and egg issue. You need to know that plane first, to score the points. But you are scoring the points, so that you can calculate that plane. I don't think that can work. – Propane

Answer 3

You can calculate the covariance matrix for all the points of the polygon (which will be a 3x3 matrix for 3D space). The normal to the polygon will be the Eigen vector corresponding to the smallest Eigen value.

Inhospitable

Any idea how can normal direction be determined? Because in calculating the covariance matrix, the order of the vertices is lost. Obviously, one can calculate a less precise normal using cross product and flip the precise normal to point in the more similar direction. Can you think of a more elegant solution? – Propane

You are right, covariance matrix will loose the connectivity and therefore normal vector may have a positive or negative sign depending on your convention. Sorry, I don't know any easy way, then the one you proposed, to get the sign. The method is very robust but expensive compared to other methods. – Inhospitable