10

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!

fduff
  • 3,671
  • 2
  • 30
  • 39
Jin Yan
  • 199
  • 12
  • 1
    why dont you represent the lines by their projection onto the plane? Otherwise it looks like you're trying to solve line intersection in threespace. (note that the cross product between the two vectors should give you a normal vector for the plane, use this to do the projection) – Steve Cox Mar 20 '14 at 20:43
  • Sounds like a case for matrices... – EOF Mar 20 '14 at 20:58
  • Wait a second, which information are given to find the coordinates of the intersection point? I mean, out of `x1`, `y1`, ..., `c2` and so on – Utkan Gezer Mar 20 '14 at 21:08
  • @SteveCox What do you mean by projecting the lines onto the plane? Which plane are you referring to? I know the normal vector of the plane on which the two lines are residing. – Jin Yan Mar 21 '14 at 00:23
  • @ThoAppelsin That's what I only have... – Jin Yan Mar 21 '14 at 00:24

3 Answers3

5

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
francis
  • 9,525
  • 2
  • 25
  • 41
3

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.

ChronoTrigger
  • 8,459
  • 1
  • 36
  • 57
  • Great answer! Very simple and definitely solved my problem. Now I'm sure that the two lines intersect with each other, what if they are parallel? Need to check whether the matrix A here is rank-deficient or not? – Jin Yan Mar 21 '14 at 00:43
  • If the lines are parallel (both vectors `(abc)_1` and `(abc)_2` are the same up to scale), matrix `A` will be rank-deficient (it will have rank 4 instead of 5). Again, if your data is noisy and the vectors are not exactly the same, `A` will be rank 5 but your solution will not be reliable. To prevent those cases, check the condition number of `A`: get its singular values and divide the largest one by the smallest. The higher, the worse. If lines are perfectly parallel, the smallest singular value will be 0. – ChronoTrigger Mar 21 '14 at 01:15
  • Of course, you could also check the angle between the vectors with the dot product, which may be even easier. – ChronoTrigger Mar 21 '14 at 01:18
  • I see. Dot product seems to be easier. So before forming this linear system, I'd better first check whether the two lines are parallel to each other to guarantee the algorithm is stable. – Jin Yan Mar 21 '14 at 01:36
2

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.

anatolyg
  • 26,506
  • 9
  • 60
  • 134