Sign in to follow this  
jollyjeffers

Finding the line of intersection between two planes

Recommended Posts

Hi all, I've been trying to figure this one out for the last hour or so, but despite searching the forums, wikipedia, google, mathworld and my geometry book I can't crack it [sad]. For example: Wikipedia: Line of intersection between two planes and Mathworld: Plane-plane intersection. They make sense, but I can't quite translate it into a concrete/usable set of equations/algorithm. I tried searching for some source code and work backwards, but didn't get very far there. Situation, I have two planes: a1*x + b1*y + c1*z + d1 = 0 a2*x + b2*y + c2*z + d2 = 0 And I want to find a line in the form p(t) = A + tB. From what I've read I gather than the directional component of the line ('B') is computed as: B = normalize( cross_product( <a1,b1,c1>, <a2,b2,c2> ) ) But then I would need the starting point, A. This is where I get stuck - and the various resources I've got cease to make any (practical) sense. I'm reckoning its gotta be pretty straight-forward and its just that my algebra is rustier than I thought ([headshake]). But equally, the fact that all of my searches turn up very little I'm wondering if its something that isn't really being used much [smile] Any help would be appreciated! Cheers, Jack

Share this post


Link to post
Share on other sites
Wouldn't A just be any point that satisfied a1*x + b1*y + c1*z + d1 = a2*x + b2*y + c2*z + d2 ? Or, is the question you're asking more how to solve for x, y, z?

Either way, I *know* my algebra is wicked rusty, so I'm not even going to take a crack at it (plus, I'm at work, and they frown on that sort of thing =). I'll leave that to someone who actually knows what they're doing.

Share this post


Link to post
Share on other sites
Perhaps I'm missing something, but:

a1*x + b1*y + c1*z + d1 = 0
a2*x + b2*y + c2*z + d2 = 0

Solving for x, y, z will give you a solution with a free variable, which defines the line of intersection. Setting free variable to 0 for example will give you a point A.

EDIT: typo, added example:

a1 = 1
b1 = 0
c1 = -2
d1 = -9

a2 = 0
b2 = 1
c2 = 1
d2 = -3


[1 0 -2 9
0 1 1 3]

Solution:

x1 = 9 + 2x3
x2 = 3 - x3
x3 free


Hope this helps

Share this post


Link to post
Share on other sites
Delphi code:

function PlaneIntersect(P, Q: TPlane): TLine3D;
var Temp: TVector3D;
begin
// Intersection of 2 planes. Direction is perpendicular to both normals
// (i.e. on both planes).
try
Result.Direction := VectorNormalize(VectorProduct(P.Normal, Q.Normal));
Result.Magnitude := -1; // unbounded
// Get the constant for both. C = -(aX + bY + cZ)
PlaneConstant(P); PlaneConstant(Q);
// Get a 3rd (temp) plane for equations
Temp.X := 1; Temp.Y := 1.2; Temp.Z := 1.3; Temp.C := 0;
while (ScalarProduct(Temp, P.Normal) = 0) or (ScalarProduct(Temp, Q.Normal) = 0) or
((P.Normal.Y <> 0) and ((Temp.X / Temp.Y) = (P.Normal.X / P.Normal.Y))) or
((Q.Normal.Y <> 0) and ((Temp.X / Temp.Y) = (Q.Normal.X / Q.Normal.Y)))
do
Temp.X := Temp.X + 1; // Make sure it's not parallel to either
// Solve 3 simultaneous equations (find intersection point)
Result.Point := SolveSimul3D(P.Normal, Q.Normal, Temp);
except
raise EMath3DError.Create('Planes do not intersect (parallel)');
end;
end;




So I just pick a random third plane and find the intersection of the line in question with that. I suspect there has to be a better way too. This code was a quick hack to get something that worked. (And don't ask me how it works, heh, I haven't looked at it for nearly two years.)

Edit: beaten by a cleverer answer.

Share this post


Link to post
Share on other sites
Quote:
Original post by jollyjeffers
Hi all,

I've been trying to figure this one out for the last hour or so, but despite searching the forums, wikipedia, google, mathworld and my geometry book I can't crack it [sad].

For example: Wikipedia: Line of intersection between two planes and Mathworld: Plane-plane intersection.

They make sense, but I can't quite translate it into a concrete/usable set of equations/algorithm. I tried searching for some source code and work backwards, but didn't get very far there.

Situation, I have two planes:

a1*x + b1*y + c1*z + d1 = 0
a2*x + b2*y + c2*z + d2 = 0

And I want to find a line in the form p(t) = A + tB.

From what I've read I gather than the directional component of the line ('B') is computed as:

B = normalize( cross_product( <a1,b1,c1>, <a2,b2,c2> ) )

But then I would need the starting point, A. This is where I get stuck - and the various resources I've got cease to make any (practical) sense.

I'm reckoning its gotta be pretty straight-forward and its just that my algebra is rustier than I thought ([headshake]). But equally, the fact that all of my searches turn up very little I'm wondering if its something that isn't really being used much [smile]

Any help would be appreciated!

Cheers,
Jack
jollyjeffers,

Here is some old code of mine for finding the line of intersection between two planes:

// Code to compute the line of intersection of two planes. On output,
// O is the origin of the line and D is the direction.

D = n1.Cross(n2);
// GetLargestAbsoluteComponent() returns the index of the element
// with the greatest absolute value
unsigned int i = D.GetLargestAbsoluteComponent();
if (Math<T>::Fabs(D[i]) < Math<T>::EPSILON)
{
O.Zero();
D.Zero();
}
else
{
unsigned int j = (i + 1) % 3;
unsigned int k = (i + 2) % 3;

T p[3];

p[i] = (T)0.0;
p[j] = (d1 * n2[k] - d2 * n1[k]) / D[i];
p[k] = (d2 * n1[j] - d1 * n2[j]) / D[i];

O.Set(p);
D.Normalize();
}




This is older code and doesn't reflect my current (and better) practices, but it's been tested and I'm almost positive it works correctly. Ideally this should be put in a function that returns false on failure (where I currently have calls to .Zero()), and true otherwise.

Let me know if you have any questions or problems converting the code. Also FYI this simply implements what you most likely read in the aforementioned references; the line of intersection is the solution space to an under-determined linear system with the 'free' variable chosen so as to make the solution maximally robust.

[Crap, three new replies in the time it took me to post that! Oh well...]

Share this post


Link to post
Share on other sites
Thanks for the quick replies - greatly appreciated!

Quote:
Original post by Falhar
what about picking 3rd random plane and compute intersection point of theese 3 planes?? :)
Yup, seen this suggested and seems to be the way forwards.

Quote:
Have you already tried looking at Soft Surfer's algorithms?!? I think there's an article in the algorithms section that talks about plane-plane intersections.
there is indeed... Thanks for the link, surprised that one didn't turn up in my googling, looks like a pretty good site [smile]

if d1 = d2 = 0, then they cross at the origin

if the cross_product( <a1,b1,c1>, <a2,b2,c2> ) is 0 (I presume it means <0,0,0> as cross product returns a vector..?) then they're parallel and there is no intersection.

To work out the other cases...

Given that B is the result of the norm/cross operation, we take a non-zero component with the largest magnitude. Thus we eliminate one of the coordinates and can directly compute the other two:

(Using Bz = 0 for example - looking for an intersection point of P<X, Y, 0>):

a1*x + b1*y + d1 = 0
a2*x + b2*y + d2 = 0

PX = (b1*d2 - b2*d1) / (a1*b2 - a2*b1)
PY = (a2*d1 - a1*d2) / (a1*b2 - a2*b1)

Right, I think I follow this now - although I still feel very dumb for some reason [headshake]

Cheers,
Jack

Share this post


Link to post
Share on other sites
I would probably use a matrix-based method, since the process is extremely straightforward and can be extended to an arbitrary number of dimensions.

The first step is to put your coefficients into a matrix (with dn on the right side of the plane equations), in this case a 2x4 matrix. Then put it into reduced row echelon form. You'll end up with a matrix that looks like this:

| 1 0 z1 d1 |
| 0 1 z2 d2 |

If the bottom row ends up being all 0's, then the planes are coincident. If only the first three elements are 0 (d2 is non-zero), then the plane normals are dependent but the distance is different so there's no intersection.

The line of intersection is then just [d1 d2 0]T + [-z1 -z2 1]T*t for all real t. Omid Ghavami provided a similar solution using free variables which is basically what t was here. I like this method a lot since all you need is a general MxN matrix class with a rref() implementation, and you can generalize this process to any dimension with any number of free variables. This formal "process" isn't very clear when you're only dealing with planes but I can explain it in a little more depth if you want.

Share this post


Link to post
Share on other sites
Thanks for the additional reply Zipster.

The matrix method does sound pretty neat - especially if you say it can be extended for an arbitrary number of dimensions.

The only bit I didn't get was 'Reduced Echelon Form' - but mathworld points me to Gaussian Elimination as a way of generating this, which seems reasonable enough.

Although, for my purpose the general algebraic method seems to be sufficient. Something that works is perfectly fine - a clever/generalised approach is just a bonus!

Being a graphics programmer I have the luxury of working in 2 or 3 dimensions for 99.99% of the time [smile]

I'll make a note of the technique should it be more appropriate in the future...

Cheers,
Jack

Share this post


Link to post
Share on other sites

Create an account or sign in to comment

You need to be a member in order to leave a comment

Create an account

Sign up for a new account in our community. It's easy!

Register a new account

Sign in

Already have an account? Sign in here.

Sign In Now

Sign in to follow this