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.

I found a solution: it's here.

The idea is to make use of vector algebra, to use the dot and cross to simply the question until this stage:

a (V1 X V2) = (P2 - P1) X V2

and calculate the a.

Note that this implementation doesn't need to have any planes or axis as reference.

The original source you mention is only for the 2d case. The implementation for perp is fine. The use of x and y are just variables not an indication of preference for a specific plane.

The same site has this for the 3d case: "http://geomalgorithms.com/a07-_distance.html"

Looks like Eberly authored a response: "https://www.geometrictools.com/Documentation/DistanceLine3Line3.pdf"

Putting this stuff here because google points to geomalgorithms and to this post.

But finding the point of intersection for two 3D line segment is not, I afraid.

I think it is. You can find the point of intersection in exactly the same way as in 2d (or any other dimension). The only difference is, that the resulting system of linear equations is more likely to have no solution (meaning the lines do not intersect).

You can solve the general equations by hand and just use your solution, or solve it programmatically, using e.g. Gaussian elemination.

I found an answer!

in an answer from above, I found these equations:

Eq#1: var s = Vector3.Dot(Vector3.Cross(dc, db), Vector3.Cross(da, db)) / Norm2(Vector3.Cross(da, db));

Eq#2: var t = Vector3.Dot(Vector3.Cross(dc, da), Vector3.Cross(da, db)) / Norm2(Vector3.Cross(da, db));

Then I modified #3rd Equation:

Eq#3:

if ((s >= 0 && s <= 1) && (k >= 0 && k <= 1))
{
   Vector3 res = new Vector3(this.A.x + da.x * s, this.A.y + da.y * s, this.A.z + da.z * s);
}

And while keeping Eq#1 and Eq#2 just the same, I created this equations:

MyEq#1: Vector3f p0 = da.mul(s).add(A<vector>); MyEq#2: Vector3f p1 = db.mul(t).add(C<vector>);

then I took a wild guess at creating these three more equations:

MyEq#3: Vector3f p0z = projUV(da, p0).add(A<vector>); MyEq#4: Vector3f p1z = projUV(db, p1).add(C<vector>);

and finally to get the subtraction of the two magnitudes of the projUV(1, 2) gives you the margin of the error between 0 and 0.001f to find whether the two lines intersect.

MyEq#5: var m = p0z.magnitude() - p1z.magnitude();

Now I mind you, this was done in Java. This explanation is not java convention ready. Just put it to work from the above equations. (Tip: Don't transform to World Space yet so that both projection of UV equations fall exactly where you want them).

And these equations are visually correct in my program.

In addition to Bobs answer:

I find on testing that the intersection() function as written solves half the original problem, which was an algorithm to find the point of intersection of two 3D line segments.

Assuming the lines are coplanar, there are 5 possible outcomes to this question:

1. The line segments are parallel, so they don't intersect, or,

2. The line segments aren't parallel, and the infinite length lines they lie upon do intersect, but the intersection point is not within the bounds of either line segment, or,

3. The lines intersect and the intersection point is within the bounds of line a but not line b, or,

4. The lines intersect and the intersection point is within the bounds of line b but not line a, or,

5. The lines intersect and the intersection point is within the bounds of both line segments.

Bob's intersection() function returns true when the lines intersect and the point of intersection is within the bounds of line a, but returns false if the lines intersect and the point of intersection is within the bounds of only line b.

But if you call intersect() twice, first with lines a then b and then a second time with lines b and a (first and second params swapped), then if both calls return true, then the intersect is contained within both line segments (case 5). If both calls return false, then neither line segment contains the intersect (case 2). If only one of the calls returns true, then the segment passed as the first parameter on that call contains the point of intersection (cases 3 or 4).

Also, if the return from the call to norm2(cross(da,db)) equals 0.0, then the line segments are parallel (case 1).

The other thing I noted in testing is that with fixed precision floating point numbers of the kind code is often implemented with, it can be quite unusual for dot(dc, cross(da,db)) to return 0.0, so returning false when its not the case might not be what you want. You might want to introduce a threshold below which the code continues to execute rather than return false. This indicates the line segments are skew in 3D, but depending on your application you might want to tolerate a small amount of skew.

The final thing I noticed was this statement in Bill's code:

ip = a.first + da * Coord(s,s,s);

That da * Coord(s,s,s) looks to be a vector times vector multiply. When I replaced it with a scalar multiple of da * s, I found it worked fine.

But in any case, many thanks to Bob. He figured out the hard part.