The algorithm to find the point of intersection of two 3D line segment

Most 3D lines do not intersect. A reliable method is to find the shortest line between two 3D lines. If the shortest line has a length of zero (or distance less than whatever tolerance you specify) then you know that the two original lines intersect.

A method for finding the shortest line between two 3D lines, written by Paul Bourke is summarized / paraphrased as follows:

In what follows a line will be defined by two points lying on it, a
point on line “a” defined by points P1 and P2 has an equation

Pa = P1 + mua (P2 - P1)

similarly a point on a second line “b” defined by points P4 and P4
will be written as

Pb = P3 + mub (P4 - P3)

The values of mua and mub range from negative to positive infinity.
The line segments between P1 P2 and P3 P4 have their corresponding mu
between 0 and 1.

There are two approaches to finding the shortest line segment between
lines “a” and “b”.

Approach one:

The first is to write down the length of the line
segment joining the two lines and then find the minimum. That is,
minimise the following

|| Pb - Pa ||^2

Substituting the equations of the lines gives

|| P1 - P3 + mua (P2 - P1) - mub (P4 - P3) ||^2

The above can then be expanded out in the (x,y,z) components.

There are conditions to be met at the minimum, the derivative with
respect to mua and mub must be zero. …the above function only has
one minima and no other minima or maxima. These two equations can then
be solved for mua and mub, the actual intersection points found by
substituting the values of mu into the original equations of the line.

Approach two:

An alternative approach but one that gives the exact same equations is
to realise that the shortest line segment between the two lines will
be perpendicular to the two lines. This allows us to write two
equations for the dot product as

(Pa - Pb) dot (P2 - P1) = 0
(Pa - Pb) dot (P4 - P3) = 0

Expanding these given the equation of the lines

( P1 - P3 + mua (P2 - P1) - mub (P4 - P3) ) dot (P2 - P1) = 0
( P1 - P3 + mua (P2 - P1) - mub (P4 - P3) ) dot (P4 - P3) = 0

Expanding these in terms of the coordinates (x,y,z) …
the result is as follows

d1321 + mua d2121 - mub d4321 = 0
d1343 + mua d4321 - mub d4343 = 0

where

dmnop = (xm - xn)(xo - xp) + (ym - yn)(yo - yp) + (zm - zn)(zo - zp)

Note that dmnop = dopmn

Finally, solving for mua gives

mua = ( d1343 d4321 - d1321 d4343 ) / ( d2121 d4343 - d4321 d4321 )

and back-substituting gives mub

mub = ( d1343 + mua d4321 ) / d4343

This method was found on Paul Bourke’s website which is an excellent geometry resource. The site has been reorganized, so scroll down to find the topic.