# Computing surface normal, floating point arithmetic

If I have a $n$-gon in $\mathbb{R^3}$, and I want to compute the surface normal, how can I get a value that minimizes error in a floating-point system? For example:

• Would I gain accuracy by first translating all vertices to the origin before performing further processing? Or would something be "lost in translation"?
• A typical way of finding the surface normal is to form two vectors from three points and compute the cross product. Which are the "best" points to choose?
-
How are you generating the vertices? If they're stored as floating point numbers $x_i\gg1$, then you have accuracy issues which will be magnified no matter how you compute the normal. If you have access to some arbitrary precision representation of the vertices, I would probably use this to come up with unit vectors along the edges and compute the normal using those. docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html – dls Jan 22 '12 at 0:02
This may be a question better suited for: scicomp.stackexchange.com – dls Jan 22 '12 at 0:04
The vertices are exported from a modeling program, and stored as 32 bit floats. I understand that error is unavoidable, but I was wondering if there is a method to ensure that the error is as small as possible - it seems that some choices ought to be better than others? – user_123abc Jan 22 '12 at 0:43
It really depends on why you think the n-gon provided really does lie in a plane. Perfection in that is impossible short of certain extensions of rational arithmetic. If I were you, I would be taking random samples from the program that supplies the n-gons, for each one calculate the $n(n-1)/2$ pairwise cross products, normalize those at unit length and adjust $\pm$ signs as necessary, and do some statistics on error. Or, find the center of gravity of the vertices and work that in... – Will Jagy Jan 22 '12 at 1:04
@Clinton: I assume you didn't intend to flag your own question as off-topic? – Zev Chonoles Jan 22 '12 at 20:06

## 1 Answer

Since your input points are given in single precision, and you can almost certainly do your computation with double-precision intermediates, your main source of error will be uncertainty in the exact position of the input points.

Ideally, I think you should want to compute the normal vector as the cross product of two diagonals: $(v_i-v_j)\times(v_k-v_l)$ with $i,j,k,l$ chosen to maximize the magnitude of the cross product. This at once selects against using very short diagonals (whose direction is more uncertain than long ones), and against using nearly parallel diagonals (where small errors in the directions of each diagonal could lead to large errors in the final normal). However, there are on the order of $n^4$ possible combinations of diagonals, so unless you have unlimited time or $n$ is always small, you don't want to try all combinations.

One possible heuristic would be to choose $i$ and $j$ as the endpoints of the longest diagonal (which can be found in $O(n^2)$ time) and then try cross products with all combinations of $k$ and $l$, also in $O(n^2)$ time.

Faster but more crudely, choose $i$ arbitrarily, then choose $j$ so the length of $v_i-v_j$ is maximal. Choose $k$ to maximize the distance from $v_k$ to the diagonal $v_iv_j$, and finally choose $l$ to maximize the cross product.

In any case, don't attempt to translate to polygon anywhere before your computations. That will just lead to intermediates of the form $(a+t)-(b+t)$ rather than $a-b$, and the former cannot possibly give better results than the latter. At best it is harmless.

-
Thanks for the reply! I asked around at some other locations, and several people pointed me to this algorithm: cs.haifa.ac.il/~gordon/plane.pdf, You might be interested in having a look. There's an additional explanation here: cs.berkeley.edu/~ug/slide/pipeline/assignments/… – user_123abc Jan 23 '12 at 21:39