7
$\begingroup$

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.

$\endgroup$
5
  • 1
    $\begingroup$ Perhaps cs.cmu.edu/~quake/robust.html could be of some use? $\endgroup$ – Marcin Łoś Feb 20 '14 at 23:07
  • $\begingroup$ you only need two points to define a plane. $\endgroup$ – TKM Feb 20 '14 at 23:27
  • $\begingroup$ @TKM If I project p on one of the vector (I assume you mean v1-v0 or v2-v0), its components will change in all cases (on the plane or not), right? $\endgroup$ – Julien-L Feb 20 '14 at 23:54
  • 1
    $\begingroup$ yes. Another method is to project p onto each of the two vectors and then add the projections to see if you get p. $\endgroup$ – TKM Feb 21 '14 at 2:38
  • 1
    $\begingroup$ Why do you think that multiplication is the source of your errors? I'm not an expert in numerical analysis, but one of the first things I learnt is that subtraction is one of the biggest sources of error. Your choice of which vertex to call $v_0$ might have a big effect. $\endgroup$ – Peter Taylor Feb 21 '14 at 8:30
7
$\begingroup$

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.

$\endgroup$
4
  • $\begingroup$ I'd love to know more about this. I still don't understand the nature of determinants. Is there some explanation or proof of this property? $\endgroup$ – Dani Barca Casafont Apr 21 '17 at 9:22
  • 2
    $\begingroup$ @Dani: A point in the plane spanned by $P_1$, $P_2$, $P_3$ (in $\mathbb R^3$ and not collinear) can be described as $P_1+\mu_2(P_2-P_1)+\mu_3(P_3-P_1)$ or $\mu_1P_1+\mu_2P_2+\mu_3P_3$ with $\mu_1=1-\mu_2-\mu_3$. So if a point $P_4$ is to be on that plane, you have to solve four equations: $\mu_1x_1+\mu_2x_2+\mu_3x_3=x_4$, same for $y$ and $z$, and $\mu_1+\mu_2+\mu_3=1$. This is an overdetermined inhomogeneous system of equations. Moving the right hand side to the left you obtain $4$ equations in $4$ variables, with $(0,0,0,0)$ always a solution and a non-zero solution only if the det is zero. $\endgroup$ – MvG Apr 23 '17 at 16:04
  • $\begingroup$ And does this work for flat planes with 4D points? $\endgroup$ – Teg Louis Dec 31 '17 at 23:19
  • $\begingroup$ @TegLouis: You can use the same approach to check whether 5 points in 4d lie on a common hyperplane (i.e. 3-dimensional affine subspace). For 4 points in 4d lying on a 2-dimensional plane this does not work, as it treats the columns as homogeneous coordinates. $\endgroup$ – MvG Jan 2 '18 at 0:33

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.