Closest points on 2d segments, passing through third 2d segment

Question

This is easier with a diagram. CaRMetal, attack!

I have two 2D line segments, P and Q. I want to find the point Px on P, and Qx on Q, such that dist(Px,Qx) is minimized. So far so good; this is a pretty straightforward task.

Now comes the wrinkle. I want to constrain Px and Qx such that the line PxQx containing them must intersect a third line segment C. (Feel free to assume that none of the original line segments intersect, BTW.)

1. There's obviously a case where the unconstrained Px and Qx already happen to satisfy the C intersection condition.
2. There's also an infeasible case, where C doesn't even contain the convex hull of P and Q. Those are trivial cases to check for.
3. In other cases, the line PxQx must contain either Ca or Cb. This does not seem to reduce cleanly to a system of linear equations.
4. And in a final set of cases, PxQx contains not only an endpoint of C, but an endpoint of P or Q as well. These seem simple to check.

It's the cases in (3) that I'm worried about, because I don't see how to get a nice closed form for it without invoking unpleasantly high-degree polynomials. I could just hurl an iterative constrained optimizer at the whole thing, of course, but I'd like to maximize performance, and high accuracy in near-degenerate cases may be important.

Answer 1

I'd like to write down a short little formula here and say, 'voila', but that's not going to happen and here's why. This problem is an extension of the find the shortest distance between two line segments which is well described by Dan Sunday in Distance between Segments and Rays. Using the labels and notation in this diagram we can parametrize P and Q by P(t) = P_0 + t(P_1 - P_0) and Q(s) = Q_0 + s(Q_1 - Q_0) where subtraction of points is done coordinate wise, i.e. Q_1 - Q_0 = (m1-m0,n1-n0). With this parametrization the problem of finding the shortest distance between the line segments P and Q is simply minimize the distance^2,

f(s,t) =  (a0 + (a1 - a0)*t - m0 - (m1 - m0)*s)^2 + (b0 + (b1-b0)*t - n0 - (n1-n0)*s)^2 

over the region in s,t space bounded by 0<=s<=1 and 0<=t<=1. (This transform avoids dealing with square roots while preserving the location of the minimum) Note that unless the segments intersect, a minimum occurs on the boundary.

We however have one more constraint---we are only considering (s,t) pairs where the line connecting P(t) and Q(s) passes through C. Fix one point Cp = (c,d) on C. Then if the line associated to a pair (s,t) passes through Cp, the vectors P(t)->Cp and Q(S)->Cp must be parallel. As this is a 2D question, we set the z direction to 0 and use the cross product (which must be zero for parallel vectors) to get that any such (s,t) must satisfy the relationship

g(s,t) = (a0 + (a1 - a0)*t - c)*(n0 + (n1 - n0)*s -d) - (b0 +(b1 - b0)*t - d)*(m0 + (m1 - m0)*s -d) = 0

So the solutions of g(s,t) = 0 are all pairs (s,t) where the line connecting their associated points go through Cp. If we parametrize C by C(r) = C_0 + r(C_1 - C_0), then we can consider that associated to each r is the set of solutions to g(s,t) = 0 where we let Cp = C(r). Plotting g(s,t) on [0,1]X[0,1] we can see what these curves look like. For the particular case in our diagram the contours increase as r increases with the contours for C(0) and C(1) forming boundaries of the feasible region. The underlying colors in the plot are the values of f(s,t). More generally, the feasible region (i.e. those values of (s,t) that are possible solutions to the problem) with our extra constraint is the points in the box [0,1]X[0,1] that are also between the contours for C(0) and C(1). It's certainly possible that neither contour intersects the box [0,1]X[0,1] at which point either then entire box is feasible or there are no solutions. As before, unless P and Q intersect in the feasible region, the minimum will occur on the boundary.

Hence this is constrained optimization problem. The boundary of the feasible region falls into one of three types:

1. Intersection point
2. Edge of the box [0,1]X[0,1]
3. Contour g(s,t) = 0 for either C(0) or C(1)

The first two are rather straight forward to check. The third however is where it gets interesting. The good news is that we are trying to minimize a quadratic function subject to a quadratic constraint. The bad news is that this means that we will need to solve a cubic to do so. The closed form for the minimum value of minimize f(s,t) subject to g(s,t) = 0 is the solution to the system of equations (reference Any calculus textbook or Wikipedia, Lagrange Multiplier)

• (partial g)/(partial s) (partial f)/partial t) = (partial g)/(partial t) (partial f)/partial s)
• g(s,t) = 0

also known as

• -2*((b0 - b1)*((n0 - n1)*s - (b0 - b1)*t + b0 - n0) + (a0 - a1)*((m0 - m1)*s - (a0 - a1)*t + a0 - m0))*((n0 - n1)*((a0 - a1)*t - a0 + c) - (m0 - m1)*((b0 - b1)*t - b0 + d)) + 2*((b0 - b1)*((m0 - m1)*s + d - m0) - (a0 - a1)*((n0 - n1)*s + d - n0))*((n0 - n1)*((n0 - n1)*s - (b0 - b1)*t + b0 - n0) + (m0 - m1)*((m0 - m1)*s - (a0 - a1)*t + a0 - m0)) = 0
• (a0 + (a1 - a0)*t - c)*(n0 + (n1 - n0)*s -d) - (b0 +(b1 - b0)*t - d)*(m0 + (m1 - m0)*s -d) = 0

There's a closed form for cubic equations and so a solution exists subject to the appropriate boundary assumptions. (In particular you need the coefficient of s and t in g to be non-zero.) The result is lengthy.

Alternately, we are minimizing a parabaloid and as such it is quadratic on the nice regular quadratic contour that is g(s,t) = 0. This is very amenable to binary search which has the added bonus of converging quickly and not requiring any square roots.

Answer 2

Px = P(t1) = Pa · (1 - t1) + Pb · t1
Qx = Q(t2) = Qa · (1 - t2) + Qb · t2

Minimize f( t1, t2 ) = |Px - Qx|²

Using the equations from "How do you detect where two line segments intersect?", we have:

t(t1, t2) = (Ca - Px) ⨯ (Cb - Ca) / (Qx - Px) ⨯ (Cb - Ca)
u(t1, t2) = (Ca - Px) ⨯ (Qx - Px) / (Qx - Px) ⨯ (Cb - Ca)

(⨯ = cross product)

Both these values must be between 0 and 1 to for the segments to intersect. t is the position along the line Px Qx, and u is the position along the C line. If you expand the formulas, they will be quotients of two linear or quadratic functions.

Since you only will be comparing t and u to zero and one, they can be simplified slightly:

t(t1, t2) = 0
(Ca - Px) ⨯ (Cb - Ca) / (Qx - Px) ⨯ (Cb - Ca) = 0
(Ca - Px) ⨯ (Cb - Ca) = 0, Qx ≠ Px

and

t(t1, t2) = 1
(Ca - Px) ⨯ (Cb - Ca) / (Qx - Px) ⨯ (Cb - Ca) = 1
(Ca - Px) ⨯ (Cb - Ca) - (Qx - Px) ⨯ (Cb - Ca) = 0
(Ca - Qx) ⨯ (Cb - Ca) = 0, Qx ≠ Px

and

u(t1, t2) = 0
(Ca - Px) ⨯ (Qx - Px) / (Qx - Px) ⨯ (Cb - Ca) = 0
(Ca - Px) ⨯ (Qx - Px) = 0, Qx ≠ Px

and

u(t1, t2) = 1
(Ca - Px) ⨯ (Qx - Px) / (Qx - Px) ⨯ (Cb - Ca) = 1
(Ca - Px) ⨯ (Qx - Px) - (Qx - Px) ⨯ (Cb - Ca) = 0
(Ca - Px) ⨯ (Qx - Px) + (Cb - Ca) ⨯ (Qx - Px) = 0
(Cb - Px) ⨯ (Qx - Px) = 0, Qx ≠ Px

We have four constraints (or eight depending on how you count):

• 0 ≤ t1 ≤ 1
• 0 ≤ t2 ≤ 1
• 0 ≤ t(t1, t2) ≤ 1
• 0 ≤ u(t1, t2) ≤ 1

You have four constrained variables, and two boundaries for each. This makes a total of eight cases you must consider, each forming a curve or line-segment in (t1, t2)-space.

The constraints form a region in (t1, t2)-space of allowed values. You must trace along the outer edges of this region and find the point that minimizes the distance between Px and Qx. As long as the segments P and Q does not intersect, the minimum will always lie on the outer border. Though in some cases, there won't be any valid solutions.

---

To find a minimum (t1, t2), you need to evaluate all candidate points:

• Interior minimum; The point where the lines P and Q cross. (1 point)
• Boundary minimum; The minimum along any of the boundary curves. (8 points)
• Intersection minimum; The points where any two boundary curves cross. (24 points)

For each point, you need to examine if it falls within all other constraints, and pick the point that generates the minimum distance between Px and Qx.

To find the interior minimum point, solve Px = Py for t1, t2. If this point falls within the other boundaries, the lines cross, and pass trough the C line. (Very unlikely)

To find the boundary minimum points, you need to look at the slope along the curves. This can be found by solving

{ ∇f(t1, t2) ⨯ ∇G(t1, t2) = 0, G(t1, t2) = k }

for t1, t2, where ∇ is the Nabla operator (Vector of the first order derivatives of the function) and G(t1, t2) = k is a boundary condition, such as t(t1, t2) = 1.

To find intersection minimum points, you need to equate two curves, and solve for t1, t2.

One way to organize this, would be to calculate the polynomial coefficients of each constraint curve, and write a function to check for intersections.

(A₁ + A₂·t1 + A₃·t1² + A₄·t2 + A₅·t1·t2 + A₆·t2²) / (A₇ + A₈·t1 + A₉·t2)

Answer 3

As for your case 3 - it looks like you can find the segment with minimal length analytically using some vector algebra and differentiation. Let's forget for the moment, that your original segments have boundary points, and express them as vectors, depending on parameter:

P(x) = P0 + p * x
Q(y) = Q0 + q * y

where P0 and Q0 are any points, lying on segments P and Q, p and q - unit vectors, and x and y - floating parameters. Also let's denote a boundary point of segment C as C0.

Any vector with the beginning on the line P and the end at the point C0 can be expressed as:

R(x) = C0 - (P0 + p * x)

Also, any vector with the beginning on the line P and the end on the line Q can be expressed as:

S(x, y) = (Q0 + q * y) - (P0 + p * x)

We want the vector S(x, y) to pass through the point C0 - it's enough to require that vectors R(x) and S(x, y) are collinear. Using this condition we can find the y variable as a function of x (and all the other parameters - P0, p, Q0, q, C0), so we'll get the vector S as a function of x only:

S(x) = S(x, y(x))

So, now we can express the norm of the S as a function of x. We can find a first derivative of this function, set it to zero and find the x value, which gives this function a minimum.

I did these calculations in a particular case when the line P lies on the X axis and got a rather complex expression for the vector length. Then I used an online symbolic differentiation to find its first derivative, which gave me a cubic equation. This equation could be solved analytically, but it's long and complex process, so a numerical solution might be another option.