plane intersection angles to determine if a point is in a plane, floating point and method questions

Question

I have been doing some work to analyze sets of 5 points in the general configuration.

[![5 points in 3 space][1]][1] [1]: https://i.stack.imgur.com/ln51Z.jpg

One question I am addressing is whether or not points A, B, E, and D are in the same plane. I am calculating this by the following method.

1. I have vectors vBA, vBD, and vBE already calculated.
   
2. vectors vBA and vBD form the plane pABD
   
3. vectors vBA and vBE form the plane pABE
   
4. planes pABD and pABE intersect along the vector vBA
   
   

I can calculate the angle at which pABD and pABE intersect along vBA by taking the vector rejection of vBD against vBA and the vector rejection of vBE against vBA and then calculating the angle between the rejection vectors.

Sine the rejection of vBD against vBA lies in the plane pABD and is perpendicular to vBA, the rejection of vBE against vBA lies in the plane pABE and is perpendicular to vBA, and both rejection vectors have their tail at the point B, the angle between these rejection vectors should be the angle at which the planes intersect.

I also calculate the plane intersection angles of pABE & pDBE along vBE and pEBA & pDBA along vBD. If all three of these plane intersection angles are equal to pi, then the points B,A,D,E are all in the same plane.

I originally tried this by calculating the cross products of the the vector combinations BAxBD, BAxBE, and BExBD. I then took that angles between each pair of cross product vectors and if the angles were equal to pi or 0, then the points were all in the same plane.

In programming, the concept of "equal to" is imprecise. It is difficult to evaluate if a real number is equal to 0, let alone pi. Realistically, you can say if it falls in some range around a number. At the moment, I am using 1x10^-9 as a tolerance, meaning that I call a value equal to pi if the value falls in the range of pi +/- 1x10^-9 inclusive. I found that with the cross product approach, I had difficulty creating real world examples (not integers) of 4 points in the same plane where all three cross product vectors were equal to pi or 0 given the above tolerance.

By using the plane intersection by vector rejections, there are fewer computational steps. When comparing the plane intersection angles, I am not going the last step to convert to radians but just taking the cosine as the dot product of the two vectors divided by the product of the vector magnitudes. This way I just have to ask if the plane intersection cosine is equal to -1. With this method, my real world examples are within 2x10^-16 of -1.

This is one such example, id x y z B -0.394666 -0.0990807 -0.610038 A 1.90877 1.16538 -0.268284 D -0.672164 -0.83449 -1.06511 E -0.707941588 0.02009757 -0.449844072 C -15.359866 -23.1923807 21.479662

ABD_EBD plane intersection cosine = -1 ABE_DBE plane intersection cosine = -1.000000000000000222 EBA_DBA plane intersection cosine = -0.999999999999999889

I have some questions regarding this methodology.

1. Is the method of using plane intersection angles valid and supportable?
   
2. Am I right that I can just evaluate if the cosine of the rejection vectors is -1 instead of taking the arc cosine and evaluating if the angle is equal to pi?
   
3. What is a reasonable tolerance to use when evaluating equivalencies given that I am using double precision on a 64-bit operating system? Can someone with experience tell me if 1x10^-9 is too small or too large?
   
4. please let me know if you want me to post any of the formulas I am using to calculate the different steps described above? Also let me know if more pictures are required.
   
   Thanks for the advice,
   LMH

Answer 1

"One question I am addressing is whether or not points A, B, E, and D are in the same plane."

You can do this much simpler than the method you are suggesting.

Method 1: Construct a matrix out of vBA, vBD, and vBE. Calculate its determinant. If it is zero (within tolerance) the vectors are coplanar and point E lies in plane ABD.

Maybe easier to understand, method 2:

$$ \vec {BE} = s \vec {BA} + t \vec {BD} $$

Dot this with the basis vectors.

$$ ( \vec {BA} \cdot \vec {BE} ) = s ( \vec {BA} \cdot \vec {BA} ) + t ( \vec {BA} \cdot \vec {BD} ) $$

$$ ( \vec {BD} \cdot \vec {BE} ) = s ( \vec {BD} \cdot \vec {BA} ) + t ( \vec {BD} \cdot \vec {BD} ) $$

( Note: $ ( \vec {BD} \cdot \vec {BA} ) = ( \vec {BA} \cdot \vec {BD} ) $. )

Solve the two equations two unknowns for s and t. Plug them into the plane equation to get the projection of $\vec {BE}$

$$ \vec {P} = s \vec {BA} + t \vec {BD} $$

Subtract the projection from the vector to get the normal displacement from the plane:

$$ \vec {D} = \vec {BE} - \vec {P} $$

The length of $ \vec {D} $ is how far the point E is from the ABD plane.

Hope this helps.

Ced

P.S. I just saw amd's comment. This is basically that.

---

Followup:

On precision. It doesn't matter how big your system's cpu registers are, what matters is the size of the number being used. A standard "double" is an eight byte floating point which is good to about 15 significant decimal digits. So 10^(-9) is huge in comparison. What really matters more is the nature of your values. If they are real world values measured to say, three significant digits, then that is the tolerance you should be working with.

Answer 2

I have coded up the second method as suggested by Cedron Dawg above and the results I am getting are not quite what I expected. I will first go through an example where point B is not in the plane of points A, D, and E.

For the following original data,

Table 1: original data
id   x          y          z
B   -0.645224  -0.591212  -0.122323
A   -1.088960  -2.853970   1.018110
D   -0.734450  -3.048430   0.922655
E   -0.518990   0.442143   0.046227
C   -0.827525   5.351198  -6.807853

we convert to vectors by taking the vector from point B to each other point, this puts point B at the origin

Table 2: vector version of data
id   x         y         z
B    0.00000   0.00000   0.00000
A   -0.44374  -2.26276   1.14043
D   -0.08923  -2.45722   1.04498
E    0.12623   1.03336   0.16855
C   -0.18230   5.94241  -6.68553

The method outlined by Cedron Dawg is based on,

$$\vec{BA} = s\vec{BA} + t\vec{BD}$$ $$equation:1$$

We create the projection of vector vBE (pBE) on the plane of vectors vBA and vBD. Since pBE will lie on the plane. If the distance between pBE and vBE is zero within floating point tolerance, then all 4 points B, A, D, E lie in the same plane. I assume that the selection of vectors here is arbitrary and we could have also found the projection of vBA on the plane of vBE,vBD, etc.

To create the projection of BE on the plane vBA,vBD, we solve the following system for the constants s and t.

$$ ( \vec{BA} \cdot \vec {BE} ) = s ( \vec {BA} \cdot \vec {BA} ) + t ( \vec {BA} \cdot \vec {BD} ) $$ $$ ( \vec {BD} \cdot \vec {BE} ) = s ( \vec {BD} \cdot \vec {BA} ) + t ( \vec {BD} \cdot \vec {BD} ) $$ $$equation:2$$

This is a system in the form of, $$ ax + by = c $$ $$ dx + ey = f $$ where x = s and y = t from equation 2.

The generalized form of the solution to this system is, $$ x = \frac{bf-ce}{bd-ae} $$ $$ y = \frac{cd-af}{bd-ae} $$ $$equation:3$$

To solve, we calculate the dot products from equation 2 based on the data in Table 2 and assign to the variables for equation 3. The data below is given in double precision.

$$ a = \vec{BA}\cdot\vec{BA} = 6.617562831748998953 $$ $$ b = \vec{BA}\cdot\vec{BD} = 6.791409871054000291 $$ $$ c = \vec{BA}\cdot\vec{BE} = -2.20202699520729972 $$ $$ d = \vec{BD}\cdot\vec{BA} = 6.791409871054000291 $$ $$ e = \vec{BD}\cdot\vec{BD} = 7.137860599084000057 $$ $$ f = \vec{BD}\cdot\vec{BE} = -2.374310923871800316 $$

we plug the above values into equation 3 and derive the solutions for the constants s and t, $$ s = x = 0.3661506250656668127 $$ $$ t = y = -0.6810149660052441378 $$

The projection of vBE on the plane vBA, vBD is thus defined as, $$\vec{pBE} = s\vec{BA} + t\vec{BD}$$ $$\vec{pBE} = 0.3661506250656668127*\vec{BA} + -0.6810149660052441378*\vec{BD}$$

plugging in the vector values from Table 2, $$\vec{pBE} = 0.3661506250656668127*(-0.443736, -2.26276, 1.14043) + -0.6810149660052441378*(-0.089226, -2.45722, 1.04498) $$ $$\vec{pBE} = (-0.162474213764138, -0.828510988373586, 0.417569157343637) + (0.0607642413567839, 1.67340359476741, -0.71164701917616) $$ $$\vec{pBE} = (-0.101709972407354, 0.844892606393824, -0.294077861832523) $$ In the precision if the input data, $$\vec{pBE} = (-0.10170, 0.84489, -0.29408)$$

The displacement, $$ \vec{D} = \vec{BE}-\vec{pBE} $$ $$ (0.126234, 1.03336, 0.16855)-(-0.10170, 0.84489, -0.29408) $$ $$ \vec{D} = (0.227934, 0.18847, 0.46263) $$ Is not zero, which means that vBE and pBE are not at the same coordinates, so the 4 points B,A,D,E are not in the same plane as expected.

This was a relatively straightforward solution, but the coordinates for the projection of vBE on plane vBA,vBD are not in the location I would expect them. Please refer to the next figure.

This figure shows pBE with the other 5 points (point C was not involved in this calculation). If you rotate this image, you can observe that point pBE appears to be in the plane vBA,vBD as expected. I would have expected to see pBE more or less underneath vBE in this figure. I have oriented the figure so that the plane vBA,vBD is more or less on the x-axis. Using the shadow analogy, with a light source on the normal vector to the plane (more or less on the y-axis), the shadow of vBE would be cast somewhere below vBE on the plane, thus creating the projection. I think that pBE should be well to the right of where it is now.

I am wondering if I have calculated this incorrectly or if I am misunderstanding the problem.

LMH