给定两条线段，找出两条线段之间的距离为d的两个点。
这类似于“两条线段之间的最短距离问题”，不同之处在于我们求解的是线段上相隔给定距离d的两个点。
每条线段由两个 3 维点组成。
我通过谷歌搜索找到的数学让我既害怕又困惑。我是一名程序员，但我很难理解解决此类问题背后的证明和分析。
输入： 2 条线段和距离d
输出：每个线段上的 2 个点彼此相距 d，如果不存在两个点，则 无
【问题讨论】：
这是一个非迭代的解决方案。我担心数学可能会激怒你，虽然这里没有什么复杂的。
首先，最简单的方法是使用整个距离的平方。
一条三维线由点 P 和 Q 描述，另一条由点 R 和 S 描述，那么说明问题的一种方式是我们要找到标量 m 和 n 使得点 a 之间的距离平方沿着第一条线的分数 m，沿着第二条线的分数 n 是给定的 dsq。
我们必须将 m 和 n 限制在 0 和 1（含）之间，这样我们的点才能真正在线段上。
如果有m和n，那么点是
X = P + m*(Q-P)
Y = R + n*(S-R)
假设我们首先找到 Dsq 的最小值和最大值。这将告诉我们是否有解决方案：如果 dsq 的给定值小于最小值或大于最大值，则没有解决方案，我们可以停止。
令出现最小值的点的 m 和 n 值为 m_min 和 n_min，最大值为 m_max 和 n_max。如果我们引入一个新变量 z（在 [0,1] 中），那么我们可以考虑 m,n 值的“线”：
m(z) = m_min + z*(m_max-m_min)
n(z) = n_min + z*(n_max-n_min)
当 z 为 0 时，这些是最小 Dsq 的值，而当 z=1 时，它们是最大 Dsq 的值。所以当我们将 z 从 0 增加到 1 时，Dsq 的值必须通过我们想要的值！也就是说，我们只需要搜索使 Dsq 成为我们想要的值的 z 值。
使问题（相对）简单的原因是 X 和 Y 之间的距离平方是 m 和 n 中的二阶多项式。具体来说，一些繁琐的代数表明，如果 Dsq 是 X 和 Y 之间的距离的平方，
Dsq = a + 2*b*m + 2*c*m + d*m*m + 2*e*m*n + f*m*m
where, in terms of dot products
a = (P-R).(P-R)
b = (P-R).(Q-P)
c =-(P-R).(S-R)
d = (Q-P).(Q-P)
e =-(Q-P).(S-R)
f = (S-R).(S-R)
最大值和最小值必须出现在任一角（(m,n)=(0,0) 或 (0,1) 或 (1,0) 或 (1,1)）或沿其中一个角边缘（对于某些 n，在 (0,n) 或 (1,n) 处，或对于某些 m，在 (m,0) 或 (m,1) 处）或在 Dsq 的导数（相对于到 m 和 n) 都是 0)。 请注意，例如，在边上说 (0,n)，我们得到了 Dsq 在 n 中的二次方，因此很容易找到其中的最大值。
此外，当我们沿着最小值和最大值之间的“线”看时，如果我们将 m(z) 和 n(z) 代入 Dsq 的公式中，经过更繁琐的代数，我们会得到z 的二次方，因此很容易找到 z 的值，它会给出所需的 Dsq 值。
嗯，这篇文章已经很长了，所以这里是实现这些想法的 C 代码。我为这些点尝试了一百万个随机值，当距离总是在最大值和最小值之间时，它总是找到合适的 3d 点。在我的（相当普通的）Linux 桌面上，这需要几秒钟。
// 3d vectors
static void v3_sub( double* P, double* Q, double* D)
{ D[0] = P[0]-Q[0];
D[1] = P[1]-Q[1];
D[2] = P[2]-Q[2];
}
static double v3_dot( double* P, double* Q)
{ return P[0]*Q[0] + P[1]*Q[1] + P[2]*Q[2];
}
// quadratic in one variable
// return *x so X -> r[0] + 2*r[1]*X + r[2]*X*X has minumum at *x
static int quad_min( const double*r, double* x)
{ if ( r[2] <= 0.0)
{ return 0;
}
*x = -r[1]/r[2];
return 1;
}
// return x so r[0] + 2*r[1]*x + r[2]*x*x == d, and whether 0<=x<=1
static int solve_quad( const double* r, double d, double* x)
{
double ap = r[0] - d;
if ( r[1] > 0.0)
{
double root1 = -(r[1] + sqrt( r[1]*r[1] - ap*r[2])); // < 0
*x = ap/root1;
}
else
{
double root1 = (-r[1] + sqrt( r[1]*r[1] - ap*r[2])); // >= 0
if ( root1 < r[2])
{ *x = root1/r[2];
}
else
{ *x = ap/root1;
}
}
return 0.0 <= *x && *x <= 1.0;
}
// quadratic in 2 variables
typedef struct
{ double a,b,c,d,e,f;
} quad2T;
static double eval_quad2( const quad2T* q, double m, double n)
{
return q->a
+ 2.0*(m*q->b + n*q->c)
+ m*m*q->d + 2.0*m*n*q->e + n*n*q->f
;
}
// eval coeffs of quad2 so that quad2(m,n) = distsq( P+m*(Q-P), R+n*(S-R))
static quad2T set_quad2( double* P, double* Q, double* R, double* S)
{
double D[3]; v3_sub( P, R, D);
double U[3]; v3_sub( Q, P, U);
double V[3]; v3_sub( S, R, V);
quad2T q;
// expansion of lengthSq( D+m*U-n*V)
q.a = v3_dot( D, D);
q.b = v3_dot( D, U);
q.c = -v3_dot( D, V);
q.d = v3_dot( U, U);
q.e = -v3_dot( U, V);
q.f = v3_dot( V, V);
return q;
}
// if gradient of q is 0 in [0,1]x[0,1], return (m,n) where it is zero
// gradient of q is 2*( q->b + m*q->d + n*q->e, q->c + m*q->e + n*q->f)
// so must solve ( q->d q->e ) * (m) = -(q->b)
// ( q->e q->f ) (n) (q->c)
static int dq_zero( const quad2T* q, double* m, double* n)
{
double det = q->d*q->f - q->e*q->e;
if ( det <= 0.0)
{ // note matrix be semi-positive definite, so negative determinant is rounding error
return 0;
}
*m = -( q->f*q->b - q->e*q->c)/det;
*n = -(-q->e*q->b + q->d*q->c)/det;
return 0.0 <= *m && *m <= 1.0
&& 0.0 <= *n && *n <= 1.0
;
}
// fill *n with minimising value, if any in [0,1], of n -> q(m0,n)
static int m_edge_min( const quad2T* q, double m0, double* n)
{
double r[3]; // coeffs of poly in n when m == m0
r[0] = q->a + 2.0*m0*q->b + m0*m0*q->d;
r[1] = q->c + m0*q->e;
r[2] = q->f;
return ( quad_min( r, n)
&& *n > 0.0 && *n < 1.0
);
}
// fill *m with minimising value, if any in [0,1], of m -> q(m,n0)
static int n_edge_min( const quad2T* q, double* m, double n0)
{
double r[3]; // coeffs of poly in m when n == n0
r[0] = q->a + 2.0*n0*q->c + n0*n0*q->f;
r[1] = q->b + n0*q->e;
r[2] = q->d;
return ( quad_min( r, m)
&& *m > 0.0 && *m < 1.0
);
}
// candidates for min, man
typedef struct
{ double m,n; // steps along lines
double d; // distance squared between points
} candT;
static int find_cands( const quad2T* q, candT* c)
{
int nc = 0;
double x, y;
// the corners
c[nc++] = (candT){ 0.0,0.0, eval_quad2( q, 0.0, 0.0)};
c[nc++] = (candT){ 0.0,1.0, eval_quad2( q, 0.0, 1.0)};
c[nc++] = (candT){ 1.0,1.0, eval_quad2( q, 1.0, 1.0)};
c[nc++] = (candT){ 1.0,0.0, eval_quad2( q, 1.0, 0.0)};
// the edges
if ( m_edge_min( q, 0.0, &x))
{ c[nc++] = (candT){ 0.0,x, eval_quad2( q, 0.0, x)};
}
if ( m_edge_min( q, 1.0, &x))
{ c[nc++] = (candT){ 1.0,x, eval_quad2( q, 1.0, x)};
}
if ( n_edge_min( q, &x, 0.0))
{ c[nc++] = (candT){ x, 0.0, eval_quad2( q, x, 0.0)};
}
if ( n_edge_min( q, &x, 1.0))
{ c[nc++] = (candT){ x, 1.0, eval_quad2( q, x, 1.0)};
}
// where the derivatives are 0
if ( dq_zero( q, &x, &y))
{ c[nc++] = (candT){ x, y, eval_quad2( q, x, y)};
}
return nc;
}
// fill in r so that
// r[0] + 2*r[1]*z + r[2]*z*z = q( minm+z*(maxm-minm), minn+x*(maxn-minn))
static void form_quad
( const quad2T* q
, double minm, double maxm
, double minn, double maxn
, double* r
)
{
double a = minm;
double c = maxm-minm;
double b = minn;
double d = maxn-minn;
r[0] = q->a + 2.0*q->b*a + 2.0*q->c*b + q->d*a*a + 2.0*q->e*a*b + q->f*b*b;
r[1] = q->b*c + q->c*d + q->d*a*c + q->e*(a*d+b*c) + q->f*b*d;
r[2] = q->d*c*c + 2.0*q->e*c*d + q->f*d*d;
}
static int find_points
( double* P, double* Q, double* R, double* S, double dsq, double* X, double* Y
)
{
double m, n;
quad2T q = set_quad2( P, Q, R, S);
candT c[9];
int nc = find_cands( &q, c); // find candidates for max and min
// find indices of max and min
int imin = 0;
int imax = 0;
for( int i=1; i<nc; ++i)
{ if ( c[i].d < c[imin].d)
{ imin = i;
}
else if ( c[i].d > c[imax].d)
{ imax = i;
}
}
// check if solution is possible -- should allow some slack here!
if ( c[imax].d < dsq || c[imin].d > dsq)
{ return 0;
}
// find solution
double r[3];
form_quad( &q, c[imin].m, c[imax].m, c[imin].n, c[imax].n, r);
double z;
if ( solve_quad( r, dsq, &z))
{ // fill in distances along
m = c[imin].m + z*(c[imax].m - c[imin].m);
n = c[imin].n + z*(c[imax].n - c[imin].n);
// compute points
for( int i=0; i<3; ++i)
{ X[i] = P[i] + m*(Q[i]-P[i]);
Y[i] = R[i] + n*(S[i]-R[i]);
}
return 1;
}
return 0;
}
【讨论】：
这可以使用初等代数来解决，只需求解二次多项式。看下面的推导：
给定由点 P1 和 P2 定义的线段 P 和由点 Q1 和 Q2 定义的线段 Q，我们可以将射线 P(t) 定义为：
P(t) = P1 + t V
其中 t 是正标量，V 是单位向量：
V = (P2 - P1) / |P2 - P1|
而线段Q(t)为：
Q(t) = Q1 + t (Q2 - Q1)
其中 t 是区间 [0 1] 中的正标量。
直线Q(t)中任意一点到直线P的最短距离由点在直线P上的投影给出。投影沿直线P的法向量。
Q(t)
|
|
P1 ------------x------------ P2
所以我们在 P 行中寻找一个点 x，使得向量 (x - Q(t)) 的长度等于 d：
|x - Q(t)|^2 = d^2
可以使用射线 P(t) 计算点 x 因为 t = (Q(t) - P1) • V：
x = P((Q(t) - P1) • V)
x = P1 + ((Q(t) - P1) • V) V
x = P1 - (P1 • V) V + (Q(t) • V) V
x = P1 - (P1 • V) V + (Q1• V) V + t ((Q2 - Q1) • V ) V
其中 (•) 是点积。
x = C + t D
在哪里
C = P1 - (P1 • V) V + (Q1• V) V
D = ((Q2 - Q1) • V ) V
现在等式如下所示：
|C + t D - Q1 - t (Q2 - Q1)|^2 = d^2
|C - Q1 + t (D - Q2 + Q1)|^2 = d^2
分组术语：
|t A + B|^2 = d^2
在哪里
A = (D - Q2 + Q1)
B = C - Q1
以我们拥有的正方形：
(t A + B) • (t A + B) = d^2
t^2 (A • A) + 2 t (A • B) + (B • B - d^2) = 0
这是一个简单的二次多项式。求解 t 我们可以得到两个值，如果两者都是复数，则没有真正的答案。如果两者都是真实的，那么我们可能有两个答案，可能由于对称性，我们必须在区间 [0 1] 中选择一个 t。
一旦我们有了 t，我们就可以使用 Q(t) 和线 P 中的对应点 x 计算线段 Q 中的点
x = P((Q(t) - P1) • V)
如果参数 (Q(t) - P1) • V 在区间 [0 L] 内，其中 L 是向量 (P2 - P1) 的长度，则 x 位于线段 P 的两端，否则x 在外面，然后找不到答案。
【讨论】：
假设 A 线的 2 个端点与 B 线上最近的端点的距离不同，我将使用蛮力法。我会选择 A 线的中心点作为 C 线的一端，并以“插入距离”的步骤在 B 线上滑动 C 线的另一端，直到我处于距离“d”的“插入距离”之内。
如果我最接近“d”的值太大，我会将 A 上 C 线的新端点设置为 A 中心点的一半 b/t 以及 A 线的最近端点到 B 线上最近的端点. 如果我最接近“d”的位置太小，我会将 A 上的新端点移动到 A 的中心点和 A 线的最远端点到 B 线最近的端点的一半处。
如果在达到最大迭代次数之前未找到可接受的值，则对“插入步骤”迭代重复此过程并返回给我最接近“d”距离的端点。然后，您可以确定您的算法是否允许足够的步数，或者是否具有接近“d”的过于严格的值。
如果 A 线的 2 个端点与 B 线最近的端点距离相同，则使用 B 线的最远端点。如果这两个端点相同，则初始步骤发生的方向是任意的。
此外，您可以使用相同的算法跳到越来越小的中点（在正确的方向上），而不是简单地在 B 线上滑动第二个端点，以节省计算次数。
【讨论】：
这是一个迭代算法，它需要一些数学知识，但不需要对数学优化有深入的了解。它很强大，但可能不是特别快。
在高层次上，这个算法就像二分搜索（技术上，三元搜索）。在每一对迭代中，我们切掉每个段剩余部分的恒定部分，如果存在一个有效的解决方案，请注意保留一个有效的解决方案。我们可以从数学上证明，随着迭代次数的增加，在极限情况下，两个段都缩小为点，如果存在这些点，这些点就是有效的解决方案。在实践中，我们会在一些迭代后停止（例如，一百次，或者当段足够短时）并返回每个段上的任意点。
此算法使用两个数学成分。第一个是计算点和线段之间距离的公式。第二个事实是，当我们沿着一段扫过一个点时，到另一段的距离会先减小后增大。
如果我有时间，我会展开这个描述。
from __future__ import division
def squared_distance_between_points(p, q):
"""Returns the squared distance between the point p and the point q."""
px, py, pz = p
qx, qy, qz = q
return (px - qx)**2 + (py - qy)**2 + (pz - qz)**2
def squared_distance_between_point_and_segment(p, q, r):
"""Returns the squared distance between the point p and the segment qr."""
px, py, pz = p
qx, qy, qz = q
rx, ry, rz = r
# Translate the points to move q to the origin (p' = p - q and r' = r - q).
px -= qx
py -= qy
pz -= qz
rx -= qx
ry -= qy
rz -= qz
# Project p' onto the line 0r'.
# The point on this line closest to p' is cr'.
c = (px * rx + py * ry + pz * rz) / (rx * rx + ry * ry + rz * rz)
# Derive c' by clamping c. The point on the segment 0r' closest to p is c'r'.
c = min(max(c, 0), 1)
# Compute the distance between p' and c'r'.
return squared_distance_between_points((px, py, pz),
(c * rx, c * ry, c * rz))
def trisect(p, q):
"""Returns the point one-third of the way from the point p to the point q."""
px, py, pz = p
qx, qy, qz = q
return ((2 * px + qx) / 3, (2 * py + qy) / 3, (2 * pz + qz) / 3)
def find_points_on_segments_at_distance(p, r, s, u, d, iterations=100):
"""Returns a point q on the segment pr and a point t on the segment su
such that the distance between q and t is approximately d.
If this is not possible (or barely possible), returns None."""
d **= 2
feasible = False
for i in range(2 * int(iterations)):
q1 = trisect(p, r)
d1 = squared_distance_between_point_and_segment(q1, s, u)
q2 = trisect(r, p)
d2 = squared_distance_between_point_and_segment(q2, s, u)
if d <= min(d1, d2):
# Use convexity to cut off one third of the search space.
if d1 <= d2:
r = q2
else:
p = q1
elif d <= max(d1, d2):
# There is certainly a solution in the middle third.
feasible = True
p = q1
r = q2
elif (d <= squared_distance_between_points(p, s)
or d <= squared_distance_between_points(p, u)):
# There is certainly a solution in the first third.
feasible = True
r = q1
elif (d <= squared_distance_between_points(r, s)
or d <= squared_distance_between_points(r, u)):
# There is certainly a solution in the last third.
feasible = True
p = q2
else:
# Definitely infeasible.
return None
# Swap the segments.
p, r, s, u = s, u, p, r
if not feasible:
return None
return p, r
【讨论】：