Check which side of a plane points are on
Asked Answered
L

4

26

I'm trying to take an array of 3D points and a plane and divide the points up into 2 arrays based on which side of the plane they are on. Before I get to heavily into debugging I wanted to post what I'm planning on doing to make sure my understanding of how to do this will work.

Basically I have the plane with 3 points and I use (pseudo code):

var v1 = new vector(plane.b.x-plane.a.x, plane.b.y-plane.a.y, plane.b.z-plane.a.z);
var v2 = new vector(plane.c.x-plane.a.x, plane.c.y-plane.a.y, plane.c.z-plane.a.z);

I take the cross product of these two vectors to get the normal vector.

Then I loop through my array of points and turn them into vectors and calculate the dot product against the normal.

Then i use the dot product to determine the side that the point is on.

Does this sound like it would work?

Laughton answered 28/3, 2013 at 17:23 Comment(2)
It does sound like it would work. However, I will also point out that you can take the "vectorizing" out of the inner loop by multiplying the point plane.a by the normal vector, to get a constant offset. This eliminates 2 subtractions, and is essentially the same as @Ali's solution: his d is - dot(plane_normal, plane.a)Peggie
I needed to relearn this while high af, and the top google results are hard to follow seeing spirits, so I wanted to add a simplified answer: (a, b, c) is the plane's normal, (x, y, z) is the point, there is no "d" it's fake and can't hurt you, move both so the plane rests at (0, 0, 0), no need to normalize anything... it's basically abovePlane = a*x + b*y + c*z > 0.0;Encephalograph
M
42

Let a*x+b*y+c*z+d=0 be the equation determining your plane.

Substitute the [x,y,z] coordinates of a point into the left hand side of the equation (I mean the a*x+b*y+c*z+d) and look at the sign of the result.

The points having the same sign are on the same side of the plane.

Honestly, I did not examine the details of what you wrote. I guess you agree that what I propose is simpler.

Metallophone answered 28/3, 2013 at 20:0 Comment(7)
Putting it in the form dot( (a,b,c,d), (x,y,z,1) ) > 0 where positive dot product is in front of the plane and negative is behind could be useful/faster.Fid
@LucasW It is the dot form: if you expand dot( (a,b,c,d), (x,y,z,1) ) you get exactly a*x+b*y+c*z+d. :)Metallophone
Could you explain this method in the language of linear algebra. What does it mean when the dot product is positive or negative?Dipteran
@dagang The plane splits the space into two half spaces. The normal vector of the plane points into one of these half spaces, lets call this half space H. The dot product is positive if the point is in H, negative otherwise. Zero if it is exactly on the plane.Metallophone
@Metallophone I realize this is old, but maybe you will reply. Can this same test be used for higher order equations of a plane that one encounters when fitting data? For example a 2nd order (in x & y) polynomial: ax^2+by^2+cz^2+dxy+exz+fyz+gx+hy+jz+k=0Iridium
@Iridium With great case, yes, but please keep in mind that higher-order things are more much complicated. How would you define for example "the same side" of a hyperboloid that has two sheets? Things get complicated with higher-order stuff, and therefore great care is necessary. Sorry, I cannot be more specific than this. :-(Metallophone
@Metallophone Thank you for your reply. You are correct and my thinking was too simple for this case. A higher order polynomial is no longer planar and is more difficult to think in terms of one side or another. Thank you for your help :)Iridium
W
3

Following the 'put points into the plane's equation and check the sign' approach given previously. The equation can be easily obtained using SymPy. I used it to find location of points (saved as numpy arrays) in a list of points.

from sympy import Point3D, Plane
plane=Plane(Point3D(point1), Point3D(point2), Point3D(point3))
for point in pointList:
        if plane.equation(x=point[0], y=point[1],z=point[2]) > 0:
            print "point is on side A"
        else:
            print "point is on side B"

I haven't tested its speed compared to other methods mentioned above but is definitely the easiest method.

Wartburg answered 30/6, 2019 at 15:3 Comment(0)
I
1

Your approach sounds good. However, when you say "and turn them into vectors", it might not be good (depending on the meaning of your sentence).

You should "turn your points into vector" by computing the difference in terms of coordinates between the current point and one of the points in the plane (for example, one of the 3 points defining the plane). As you wrote it, it sounds like you might have misunderstood that ; but apart from that, it's ok!

Irrigate answered 28/3, 2013 at 20:14 Comment(0)
K
0

take into account the normal vector of the plane

example: for the point A=[-243.815437431962, -41.7407630281635, 10.0] equation= -2663.1860000000006*Z +21305.488000000005=0 RESULt POSITIVE

but if equation= 2663.1860000000006*Z -21305.488000000005=0

RESULT NEGATIVE

Kiyokokiyoshi answered 4/4, 2022 at 1:22 Comment(0)

© 2022 - 2024 — McMap. All rights reserved.