Check which side of a plane points are on

Question

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?

Solution 1

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.

Solution 2

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.

Solution 3

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!

Comments

[Lucas W]

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.

[Ali]

@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. :)

[Dagang]

Could you explain this method in the language of linear algebra. What does it mean when the dot product is positive or negative?

[Ali]

@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.

[Saladsamurai]

@Ali 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=0

[Ali]

@Saladsamurai 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. :-(

[Saladsamurai]

@Ali 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 :)