Calculating the shortest distance between two lines (line segments) in 3D

Question

I have two line segments: X1,Y1,Z1 - X2,Y2,Z2 And X3,Y3,Z3 - X4,Y4,Z4

I am trying to find the shortest distance between the two segments.

I have been looking for a solution for hours, but all of them seem to work with lines rather than line segments.

Any ideas how to go about this, or any sources of furmulae?

Evolutionary

Take a look at the paper Robust Computation of Distance Between Line Segments – Triboluminescence

Answer 1

One basic approach is the same as computing the shortest distance between 2 lines, with one exception.

If you look at most algorithms for finding the shortest distance between 2 lines, you'll find that it finds the points on each line that are the closest, then computes the distance from them.

The trick to extend this to segments (or rays), is to see if that point is beyond one of the end points of the line, and if so, use the end point instead of the actual closest point on the infinite line.

For a concrete sample, see:

http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm

More specifically:

http://softsurfer.com/Archive/algorithm_0106/algorithm_0106.htm#dist3D_Segment_to_Segment()

Wickliffe

The SoftSurfer code is fairly good. It has trouble with almost parallel lines. I ended up writing a check for "almost parallel" lines. Once I did this, it worked well. I'm not sure why the check for "almost parallel" lines built into SoftSurfer didn't work out. Just thought the next user would like to know.... – Gelatinoid

I found the same thing as @TimPerry. I used this parallel check: bool parallel = dot(u, v) / (u.length() * v.length()) > 1 - SMALL_NUM; Modify for your vector library. – Miscible

@ReedCopsey I tried dist3D_Segment_to_Segment() with the following line segments: LineSegment lineSegmentA; lineSegmentA.startPoint = PointType(0.,0.,0.); lineSegmentA.endPoint = PointType(1.,0.,0.); LineSegment lineSegmentB; lineSegmentB.startPoint = PointType(.5,0.2,0.); lineSegmentB.endPoint = PointType(.5,1.2,0.); The closest distance should be from the endPoint of lineSegmentB to the center point of lineSegmentA, which should have value "0.2". However, the function returns the distance as "0.538". Any thoughts? – Touchandgo

@DavidDoria You have almost certainly transcribed the function incorrectly. That's the value you'd get if the D < SMALL_NUM route is chosen. It should not be chosen and indeed is not chosen. You need to make sure you transcribe the code correctly. – Ungenerous

Those links are dead. Are there new links you can source? Is there a new source that can be linked? – Curet

This is why you don't post links... – Wuhu

Answer 2

I'll answer this in terms of matlab, but other programming environments can be used. I'll add that this solution is valid to solve the problem in any number of dimensions (>= 3).

Assume that we have two line segments in space, PQ and RS. Here are a few random sets of points.

> P = randn(1,3)
P =
     -0.43256      -1.6656      0.12533

> Q = randn(1,3)
Q =
      0.28768      -1.1465       1.1909

> R = randn(1,3)
R =
       1.1892    -0.037633      0.32729

> S = randn(1,3)
S =
      0.17464     -0.18671      0.72579

The infinite line PQ(t) is easily defined as

PQ(u) = P + u*(Q-P)

Likewise, we have

RS(v) = R + v*(S-R)

See that for each line, when the parameter is at 0 or 1, we get one of the original endpoints on the line returned. Thus, we know that PQ(0) == P, PQ(1) == Q, RS(0) == R, and RS(1) == S.

This way of defining a line parametrically is very useful in many contexts.

Next, imagine we were looking down along line PQ. Can we find the point of smallest distance from the line segment RS to the infinite line PQ? This is most easily done by a projection into the null space of line PQ.

> N = null(P-Q)
N =
     -0.37428     -0.76828
       0.9078     -0.18927
     -0.18927      0.61149

Thus, null(P-Q) is a pair of basis vectors that span the two dimensional subspace orthogonal to the line PQ.

> r = (R-P)*N
r =
      0.83265      -1.4306

> s = (S-P)*N
s =
       1.0016     -0.37923

Essentially what we have done is to project the vector RS into the 2 dimensional subspace (plane) orthogonal to the line PQ. By subtracting off P (a point on line PQ) to get r and s, we ensure that the infinite line passes through the origin in this projection plane.

So really, we have reduced this to finding the minimum distance from the line rs(v) to the origin (0,0) in the projection plane. Recall that the line rs(v) is defined by the parameter v as:

rs(v) = r + v*(s-r)

The normal vector to the line rs(v) will give us what we need. Since we have reduced this to 2 dimensions because the original space was 3-d, we can do it simply. Otherwise, I'd just have used null again. This little trick works in 2-d:

> n = (s - r)*[0 -1;1 0];
> n = n/norm(n);

n is now a vector with unit length. The distance from the infinite line rs(v) to the origin is simple.

> d = dot(n,r)
d =
       1.0491

See that I could also have used s, to get the same distance. The actual distance is abs(d), but as it turns out, d was positive here anyway.

> d = dot(n,s)
d =
       1.0491

Can we determine v from this? Yes. Recall that the origin is a distance of d units from the line that connects points r and s. Therefore we can write dn = r + v(s-r), for some value of the scalar v. Form the dot product of each side of this equation with the vector (s-r), and solve for v.

> v = dot(s-r,d*n-r)/dot(s-r,s-r)
v =
       1.2024

This tells us that the closest approach of the line segment rs to the origin happened outside the end points of the line segment. So really the closest point on rs to the origin was the point rs(1) = s.

Backing out from the projection, this tells us that the closest point on line segment RS to the infinite line PQ was the point S.

There is one more step in the analysis to take. What is the closest point on the line segment PQ? Does this point fall inside the line segment, or does it too fall outside the endpoints?

We project the point S onto the line PQ. (This expression for u is easily enough derived from similar logic as I did before. Note here that I've used \ to do the work.)

> u = (Q-P)'\((S - (S*N)*N') - P)'
u =
      0.95903

See that u lies in the interval [0,1]. We have solved the problem. The point on line PQ is

> P + u*(Q-P)
ans =
      0.25817      -1.1677       1.1473

And, the distance between closest points on the two line segments was

> norm(P + u*(Q-P) - S)
ans =
        1.071

Of course, all of this can be compressed into just a few short lines of code. But it helps to expand it all out to gain understanding of how it works.

Ulcerate

Answer 3

I would parameterize both line segments to use one parameter each, bound between 0 and 1, inclusive. Then find the difference between both line functions and use that as the objective function in a linear optimization problem with the parameters as variables.

So say you have a line from (0,0,0) to (1,0,0) and another from (0,1,0) to (0,0,0) (Yeah, I'm using easy ones). The lines can be parameterized like (1*t,0*t,0*t) where t lies in [0,1] and (0*s,1*s,0*s) where s lies in [0,1], independent of t.

Then you need to minimize ||(1*t,1*s,0)|| where t, s lie in [0,1]. That's a pretty simple problem to solve.

Jackquelinejackrabbit

Given line segments from p1 to p2 and from q1 to q2 you need to compute all of the following distances and take the minimum: (line1, line2), (p1, line2), (p1, q1), (p1, q2), (p2, line2), (p2, q1), (p2, q2), (line1, q1), (line1, q2). (Or maybe you can show mathematically that some can be eliminated.) – Katiakatie

Answer 4

Finding the distance between two finite lines based on finding between two infinite lines and then bound the infinite lines to the finite lines doesn't work always. for example try this points

Q=[5 2 0]
P=[2 2 0]
S=[3 3.25 0]
R=[0 3 0]

Based on infinite approach the algorithm select R and P for distance calculation (distance=2.2361), but somewhere in the middle of R and S has got a closer distance to the P point. Apparently, selecting P and [2 3.166] from R to S line has lower distance of 1.1666. Even this answer could get better by precise calculation and finding orthogonal line from P to R S line.

Shamus

I have tried different more points and methods and this is my finding till now. When you are using project method to find distance between two finite lines you must perform projection in either side. It means that you have to first project RS to PQ first and then vise-versa to get the right answer. In this method two different distance is calculated and the lowest one is what you are seeking. – Shamus

Answer 5

This question is the topic of the article On fast computation of distance between line segments by Vladimir J. Lumelksy 1985. It goes even further by finding not only the minimal Euclidean distance (MinD) but a point on each segment separated by that distance.

The general algorithm is as follows:

1. Compute the global MinD (global means the distance between two infinite lines containing the segments) and coordinates of both points (bases) of the line of minimum distances, see skew lines; if both bases lie inside the segments, then actual MinD is equal to the global MinD; otherwise, continue.
2. Compute distances between the endpoints of both segments (a total of four distances).
3. Compute coordinates of the base points of perpendiculars from the endpoints of one segment onto the other segment; compute the lengths of those perpendiculars whose base points lie inside the corresponding segments (up to four base point, and four distances).
4. Out of the remaining distances, the smallest is the sought actual MinD. Altogether, this represents the computation of six points and of nine distances.

The article then describes and proves how to reduce the amount of tests based on the data received in initial steps of the algorithm and how to handle degenerate cases (e.g. equal endpoints of a segment).

C-language implementation by Eric Larsen can be found here, see SegPoints() function.

Triboluminescence

Answer 6

How about extending the line segments into infinite lines and find the shortest distance between the two lines. Then find the points on each line that are the end points of the shortest distance line segment.

If the point for each line is on the original line segment, then the you have the answer. If a point for each line is not on the original segment, then the point is one of the original line segments' end points.

Appetizer

Answer 7

First, find the closest approach Line Segment bridging between their extended lines. Let's call this LineSeg BR.

If BR.endPt1 falls on LS1 and BR.endPt2 falls on LS2, you're done...just calculate the length of BR.

If the bridge BR intersects LS1 but not LS2, use the shorter of these two distances: smallerOf(dist(BR.endPt1, LS2.endPt1), dist(BR.endPt1, LS2.endPt2))

If the bridge BR intersects LS2 but not LS1, use the shorter of these two distances: smallerOf(dist(BR.endPt2, LS1.endPt1), dist(BR.endPt2, LS1.endPt2))

If none of these conditions hold, the closest distance is the closest pairing of endpoints on opposite Line Segs.

Arsenite