How to find out the intersection of two coplanar lines in C

Question

I have two 3D lines which lie on the same plane. line1 is defined by a point (x1, y1, z1) and its direction vector (a1, b1, c1) while line2 is defined by a point (x2, y2, z2) and its direction vector (a2, b2, c2). Then the parametric equations for both lines are

x = x1 + a1*t;         x = x2 + a2*s;
y = y1 + b1*t;         y = y2 + b2*s;
z = z1 + c1*t;         z = z2 + c2*s;

If both direction vectors are nonzeros, we can find out the location of intersection node easily by equating the right-hand-side of the equations above and solving t and s from either two of the three. However, it's possible that a1 b1 c1 a2 b2 c2 are not all nonzero so that I can't solve those equations in the same way. My current thought is to deal with this issue case by case, like

case1: a1 = 0, others are nonzero
case2: a2 = 0, others are nonzero
case3: b1 = 0, others are nonzero
...

However, there are so many cases in total and the implementation would become messy. Is there any good ways to tackle this problem? Any reference? Thanks a lot!

Answer 1

You can solve this as a linear system:

| 1 0 0 -a1   0 | | x |   | x1 |
| 0 1 0 -b1   0 | | y |   | y1 |
| 0 0 1 -c1   0 | | z | = | z1 |
| 1 0 0   0 -a2 | | s |   | x2 |
| 0 1 0   0 -b2 | | t |   | y2 |
| 0 0 1   0 -c2 |         | z2 |

x y z is the intersection point, and s t are the coefficients of the vectors. This solves the same equation that @francis wrote, with the advantage that it also obtains the solution that minimizes the error in case your data are not perfect.

This equation is usually expressed as Ax=b, and can be solved by doing x = A^(-1) * b, where A^(-1) is the pseudo-inverse of A. All the linear algebra libraries implement some function to solve systems like this, so don't worry.

Answer 2

It is much more practical to see this as a vector equation. A dot . is a scalar product and A,n,B,m are vectors describing the lines. Point A is on the first line of direction n. Directions are normalized : n.n=1 and m.m=1. The point of intersection C is such that :

 C=A+nt=B+ms

where t and s are scalar parameters to be computed.

Therefore (.n) :

A.n+ t=B.n+m.n s
t= (B-A).n+m.n s

And (.m):

 A.m+n.m t=B.m+ s
 A.m+n.m (B-A).n+(m.n)^2 s=B.m+ s
 n.m(B-A).n+(A-B).m=(1-(m.n)^2).s

Since n.n=m.m=1 and n and m are not aligned, (m.n)^2<1 :

 s=[n.m(B-A).n+(A-B).m]/[1-(m.n)^2]
 t= (B-A).n+m.n s

Answer 3

It might be vital to remember that calculations are never exact, and small deviations in your constants and calculations can make your lines not exactly intersect.

Therefore, let's solve a more general problem - find the values of t and s for which the distance between the corresponding points in the lines is minimal. This is clearly a task for calculus, and it's easy (because linear functions are the easiest ones in calculus).

So the points are

[xyz1]+[abc1]*t
and
[xyz2]+[abc2]*s

(here [xyz1] is a 3-vector [x1, y1, z1] and so on)

The (square of) the distance between them:

([abc1]*t - [abc2]*s + [xyz1]-[xyz2])^2

(here ^2 is a scalar product of a 3-vector with itself)

Let's find a derivative of this with respect to t:

[abc1] * ([abc1]*t - [abc2]*s + [xyz1]-[xyz2])    (multiplied by 2, but this doesn't matter)

(here the first * is a scalar product, and the other *s are regular multiplications between a vector and a number)

The derivative should be equal to zero at the minimum point:

[abc1] * ([abc1]*t - [abc2]*s + [xyz1]-[xyz2]) = 0

Let's use the derivative with respect to s too - we want it to be zero too.

 [abc1]*[abc1]*t - [abc1]*[abc2]*s = -[abc1]*([xyz1]-[xyz2])
-[abc2]*[abc1]*t + [abc2]*[abc2]*s =  [abc2]*([xyz1]-[xyz2])

From here, let's find t and s.

Then, let's find the two points that correspond to these t and s. If all calculations were ideal, these points would coincide. However, at this point you are practically guaranteed to get some small deviations, so take and of these points as your result (intersection of the two lines).

It might be better to take the average of these points, to make the result symmetrical.