[Math] Check if a point is on a plane? (Minimize the use of multiplications and divisions)

3dfloating pointgeometry

In $\mathbb R3$, given a plane $\mathcal P$ defined by three 3D points points $v_0, v_1, v_2$, I want to check if another point $p$ belongs to that plane, while avoiding the use of multiplications and divisions as much as possible.

The reason is to mitigate the floating point errors incurred by the computer representation of decimal numbers.

My current method is:

  1. Compute the general form of the plane equation $ax+by+cz+d=0$

    Where $a,b,c$ are the components of the plane's unit normal vector
    $N={(v_1-v_0)\times(v_2-v_0) \over \|(v_1-v_0)\times(v_2-v_0)\|}$

    And $d=N.v_0$

  2. Plug point $p$ into the plane equation: $res=a.p_x+b.p_y+c.p_z-d$

    If the result is null, the point lies on the plane

    Because of floating point errors, I actually check if the result is "almost" null: $|res|<\epsilon$

However, in certain cases, when I plug $v_2$ into the plane equation, I find that it does not belong to the plane, even though I used $v_2$ to compute the equation in the first place. (I obtain a result bigger than my $\epsilon$.)

This is due to floating point errors. See my question on Stack Overflow: https://stackoverflow.com/q/21916606/143504

So I am looking for an alternate method that would mitigate these errors.

Best Answer

Four points lie in a common plane if the determinant

$$\begin{vmatrix} x_1 & x_2 & x_3 & x_4 \\ y_1 & y_2 & y_3 & y_4 \\ z_1 & z_2 & z_3 & z_4 \\ 1 & 1 & 1 & 1 \end{vmatrix}$$

is zero. This has the benefit that it requires no divisions at all. It does require quite some multiplications, though.

You might however want to have a look at Shewchuk's page on Adaptive Precision Floating-Point Arithmetic and Fast Robust Predicates for Computational Geometry, where you can find papers on how to evaluate these predicates exactly, and also a C implementation doing just this.

As an alternative, you might want to use one of CGALs exact predicates kernels.

Note however that if your input points are double coordinates, then a point lying exactly on a plane spanned by three others is rather unlikely. So be sure that you don't make conceptual mistakes there.

Your statement about using $\epsilon$ when interpreting the result shows that you are aware of the numeric imprecision involved. However, I'm unsure how you'd decide on what values for $\epsilon$ are reasonable. That may depend on your application. Only then can you judge in which direction you'd rather err.

Related Question