How do I get three non-colinear points on a plane?

Question

I'm trying to implement at line-plane intersection algorithm. According to Wikipedia I need three non-colinear points on the plane to do that.

I therefore tried implementing this algorithm in C++, however. Something is definitely wrong, cause it makes no sense that I can choose whatever x and y coordinates and they'll fit in the plane. What if the plane is vertical and along the x-axis? No point with y=1 would then be in the plane.

I realize that this problem has been posted a lot on StackOverflow, and I see lots of solutions where the plane is defined by 3 points. But I only have a normal and a position. And I can't test my line-plane intersection algorithm before I sort out my non-colinear point finder.

The problem right now, is that I'm dividing by normal.z, and that obviously won't work when normal.z is 0.

I'm testing with this plane: Plane* p = new Plane(Color(), Vec3d(0.0,0.0,0.0), Vec3d(0.0,1.0,0.0)); // second parameter : position, third parameter : normal

The current code gives this incorrect answer:

{0 , 0 , 0} // alright, this is the original
{12.8377 , 17.2728 , -inf} // obviously this is not a non-colinear point on the given plane

Here's my code:

std::vector<Vec3d>* Plane::getThreeNonColinearPoints() {
    std::vector<Vec3d>* v = new std::vector<Vec3d>();

    v->push_back(Vec3d(position.x, position.y, position.z)); // original position can serve as one of the three non-colinear points.

    srandom(time(NULL));

    double rx, ry, rz, start;

    rx = Plane::fRand(10.0, 20.0);
    ry = Plane::fRand(10.0, 20.0);
    // Formula from here: http://en.wikipedia.org/wiki/Plane_(geometry)#Definition_with_a_point_and_a_normal_vector
    // nx(x-x0) + ny(y-y0) + nz(z-z0) = 0
    // |-----------------| <- this is "start"
    //I'll try to insert position as x0,y0,z0 and normal as nx,ny,nz, and solve the equation
    start = normal.x * (rx - position.x) + normal.y * (ry - position.y);
    // nz(z-z0) = -start
    start = -start;
    // (z-z0) = start/nz
    start /= normal.z; // division by zero
    // z = start+z0
    start += position.z;
    rz = start;

    v->push_back(Vec3d(rx, ry, rz));

    // TODO one more point

    return v;
}

I realize that I might be trying to solve this totally wrong. If so, please link a concrete implementation of this. I'm sure it must exist, when I see so many line-plane intersection implementations.

Thanks in advance.

Answer 1

Where N=(Nx,Ny,Nz) is the normal, you could project the points N, (Ny,Nz,Nx), (Nz,Nx,Ny) onto the plane: they're guaranteed to be distinct.

Alternatively, if P and Q are on the plane, P+t(Q-P)xN is also on the plane for any t!=0 where x is the cross product.

Alternatively if M!=N is an arbitrary vector, K=MxN and L=KxN are colinear with the plane and any point p on the plane can be written as p=Origin+sK+tL for some s,t.

Answer 2

A plane can be defined with several ways. Typically a point on the plane and a normal vector is used. To get the normal vector from three points (P1, P2, P3 ) take the cross product of the side of the triangle

P1 = {x1, y1, z1};
P2 = {x2, y2, z2};
P3 = {x3, y3, z3};

N = UNIT( CROSS( P2-P1, P3-P1 ) );
Plane P = { P1, N }

The reverse, to go from a point P1 and normal N to three points, you start from any direction G not along the normal N such that DOT(G,N)!=0. The two orthogonal directions along the plane are then

//try G={0,0,1} or {0,1,0} or {1,0,0}
G = {0,0,1};
if( MAG(CROSS(G,N))<TINY ) { G = {0,1,0}; }
if( MAG(CROSS(G,N))<TINY ) { G = {1,0,0}; }
U = UNIT( CROSS(N, G) );  
V = CROSS(U,N);
P2 = P1 + U;
P3 = P1 + V;

A line is defined by a point and a direction. Typically two points (Q1, Q2) define the line

Q1 = {x1, y1, z1};
Q2 = {x2, y2, z2};
E = UNIT( Q2-Q1 );
Line L = { Q1, E }

The intersection of the line and plane are defined by the point on the line r=Q1+t*E that intersects the plane such that DOT(r-P1,N)=0. This is solved for the scalar distance t along the line as

t = DOT(P1-Q1,N)/DOT(E,N);

and the location as

r = Q1+(t*E);

NOTE: The DOT() returns the dot-product of two vector, CROSS() the cross-product, and UNIT() the unit vector (with magnitude=1).

DOT(P,Q) = P[0]*Q[0]+P[1]*Q[1]+P[2]*Q[2];
CROSS(P,Q) = { P[1]*Q[2]-P[2]*Q[1], P[2]*Q[0]-P[0]*Q[2], P[0]*Q[1]-P[1]*Q[0] };
UNIT(P) = {P[0]/sqrt(DOT(P,P)), P[1]/sqrt(DOT(P,P)), P[2]/sqrt(DOT(P,P))};
t*P =  { t*P[0], t*P[1], t*P[2] };
MAG(P) = sqrt(P[0]*P[0]+P[1]*P[1]+P[2]*P[2]);

Answer 3

One approach you may find easy to implement is to see where the plane intersects the coordinate axes. For the plane given by the equationaX + bY + cZ - d = 0 hold two variables at 0 and solve for the third. So the solutions would be (assuming a, b, c, and d are all non-zero):

(d/a, 0, 0)
(0, d/b, 0)
(0, 0, d/c)

You will need to consider the cases where one or more of the coefficients are 0 so you don't get a degenerate or colinear solutions. As an example if exactly one of the coefficients is 0 (say a=0) you instead use

(1, d/b, 0)
(0, d/b, 0)
(0, 0, d/c)

If exactly two of the coefficients are 0 (say a=0 and b=0) you can use:

(1, 0, d/c)
(0, 1, d/c)
(0, 0, d/c)

If d=0, the plane intersects the three axes at the origin, and so you can use:

(1, 0, -a/c)
(0, -c/b, 1)
(-b/a, 1, 0)

You will need to work out simular cases for d and exactly one other coefficient being 0, as well as d and two others being 0. There should be a total of 16 cases, but there are a few things that come to mind which should make that somewhat more manageable.