Closest points on 2d segments, passing through third 2d segment
computational-geometry, geometry, intersection, math, mathematical-optimization
Solution
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:
- Intersection point
- Edge of the box `[0,1]X[0,1]`
- 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.
Problem
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.) - There's obviously a case where the unconstrained `Px` and `Qx` already happen to satisfy the `C` intersection condition. - 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. - 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. - 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.