Floating-point error when checking for coplanar 3D points

Question

I am looking for an algorithm to check if a point is coplanar with a given 3D plane, defined out of three vertices, while minimizing floating point errors.

I would like to minimize the amount of multiplications and division to mitigate floating point errors.

My implementation uses floats, I cannot go double.

I cannot use an external library.

---

My current method suffers from these errors:

I have code defining a plane using the general form of the plane equation:

ax + by + cz + d = 0

I compute these coefficients using three 3D vertices v0, v1 and v2 as follow:

// Pseudo-code to define a plane (with class Vector3 defining a vector in 3D)
Vector3 A = v1 - v0;
Vector3 B = v2 - v0;
Vector3 N = cross_product(A,B); // Normal vector
N.Normalize();                  // Unit normal vector storing coefs. a, b, c
float d = dot_product(N,v0);

To check if another vertex p is coplanar, I plug the point into the plane equation and check if the result is 0:

// Pseudo-code for coplanar test:
bool is_coplanar()
{
    float res = N.x()*p.x() + N.y()*p.y() + N.z()*p.z() - d;
    return true if res is "almost" null; // "almost" is: abs(res)<EPSILON
}

---

My code fails in this case:

v0 = [-8.50001907, 0, 323]
v1 = [8.49998093, 0, 323]
v2 = [-8.50001907, 1.49999976, 322.598083]

Then the plane coefficients are:

N = [-0, 0.258814692, 0.965926945]
d = 311.994415

And when I plug the point v2, I find a result "far" from 0 (although v2 was used to define the plane):

res = -3.05175781e-05

My EPSILON is currently 1e-5.

Tested on compiler qcc 4.4.2 (QNX Momentics, similar to gcc). With no optimization -O0.

Answer 1

Such geometric predicates suffer in a lot of ways from floating point errors. The only industrial strength solution is to use adaptable arithmetic filtering (provided that a robust implementation of the coplanar test is not covering you).

Luckily such implementations (that would take quite some time to write) are already available. In the previous link the orient3d predicate does what you need: Given 3 plane forming points, decide whether a 4th one lies above,below or on the plane

If such an implementation is an overkill, check the simple one. It offers 4 in total:

orient3dfast() Approximate 3D orientation test. Nonrobust.
orient3dexact() Exact 3D orientation test. Robust.
orient3dslow() Another exact 3D orientation test. Robust.
orient3d() Adaptive exact 3D orientation test. Robust.

Disclaimer: The code listing is provided as a tutorial of the mathematical concepts and programming techniques needed to reach a robust solution. I'm neither suggesting nor implying copy-pasting anything.