Re: [Algorithms] 3d Lines Intersection

[<ro...@do...>]

Albert_J  wrote:

>Hi,
>
>Is there a good algo for finding a point of intersection between two *3D lines*?
>The lines is specified by their point (x,y,z) and their direction (dx,dy,dz)
>
>I've searched the web, but all of them seems to explain about 2D lines intersection only.

The problem is that in 3D almost all pairs of lines do not intersect
(as opposed to 2D where every pair of non parallel lines intersects).

Even if the logic of the construction of the two lines is such that
you know they must intersect, any computational intersection test will
be subject to rounding errors which, in all probability, will yield a
negative result.

The best you can do computationally  is to find the two points, one on
each line, which are nearest.  You may or may not want to call the
lines "intersecting" depending on the distance between these two
points and your error tolerance.

So if your two lines are specified by  vectors X1, DX1, X2, DX2,
namely
          X1= (x1, y1, z1), 
          DX1 = (dx1, dy1, dz1) 
         X2 = (x2, y2, z2), 
        DX2 = (dx2, dy2, dz2)

write down the vector parametric representations of the general point
one each of the two lines

       X1 + s DX1
       X2  + t DX2.

where s and t are INDEPENDENT real parameters.

You want to find the values of s and t that make the distance between
these two points as small as possible.  Write down the expression for
the distance between the two points, or rather the square of the
distance ( more convenient to avoid square roots), in terms of these
two parameters.   Using elementary vector geometry this squared
distance is 

  distance squared = 
 | (X1 + sDX1) - (X2 +  tDX2) |^2
 = |(X1 - X2) + s DX1 - t DX2)|^2
 = ((X1 - X2) + s DX1 - t DX2) . ( (X1 - X2) + s DX1 - t DX2))
where . means dot product.

Apply the elementary algebraic properties of dot product and collect
terms that are like in s and t and you get
 s^2 DX1.DX1 + t^2 DX2.DX2 - 2st DX1.DX2 + 
 s DX1. (X1 - X2) -  t DX2. (X1 - X2) + (X1 - X2).(X1 - X2)

This distance-squared function is of the form
as^2 + bt^2 + c st + ds + et + f
where a, b, c, d, e, f are known scalars computed from your line data
and s and t are independent variables.

This is a general quadratic function of two independent variables s
and t.  We seek to find the values of s and t that make it as small as
possible.  The best tool for this minization problem is the elementary
differential calculus of two variables, which tells us that the
minimum value occurs where the first partial derivatives with respect
to each of the variables is zero.  So differentiating first with
respect to s and then with respect to t gives two equations
 
	2as + ct + d = 0
                2bt  + cs + e = 0

This is a system of two linear equations in two unknowns s and t,
which you know how to solve from high school algebra.  It will have a
unique solution if the lines are not parallel.

Plug the resulting values of s and t back into the distance-squared
formula and see if that is close enough to zero so that you want to
call it an intersection.  In any case, if you plug the values of s and
t back into the parametric expressions for the two lines you get the
two closest points, one on each line.

P.S. You can avoid calculus and still find the two nearest points by
using the fact from Euclidean solid geometry that the closest two
points are the two feet of the mutual perpendicular between the two
lines, the mutual perpendicular being unique if the lines are not
parallel.  Write down the expression for the vector joining a general
point of each line, as a function of s and t, and write down the
equations that say that the dot product of this vector with each of
DX1 and DX2 is zero (this is the mutual perpendicularity condition). 
 
i.e. 	DX1. ((X1 - X2) + s DX1 + t DX2) = 0
                DX2. ((X1 - X2) + s DX1 + t DX2)) = 0

 The resulting two equations is a 2x2 linear system equivalent to the
one above.  The reason I chose the calculus approach is because I
didn't just want to find the two closest points, but also to test for
(near) intersection, which means I would have had to know the distance
function anyway, so I chose to start with the distance function, and
followed my nose from there.