Вот неитеративное решение.Боюсь, что математика может вас раздражать, хотя здесь нет ничего сложного.
Во-первых, проще всего работать с квадратами на всем протяжении.
Это одна трехмерная линия, описанная точками P иQ, а другой по точкам R и S, то один из способов постановки проблемы состоит в том, что мы хотим найти скаляры m и n, чтобы расстояние возводилось в квадрат между точкой дроби 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 должно пройти через значение, которое мы хотим!То есть нам просто нужно найти значение z, которое делает Dsq желаемым значением.
Что делает проблему (относительно) прямой, так это то, что квадрат расстояния между 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)) или вдоль одного из ребер (в (0, n) или (1, n) для некоторого n, или (m, 0)или (m, 1) для некоторого m) или в точке посередине, где оба производных Dsq (по m и n) равны 0).Обратите внимание, например, что на ребре скажем (0, n) мы получаем квадратик по n для Dsq, поэтому легко найти максимум этого значения.
Более того, когда мы посмотрим вдоль линииМежду минимальными и максимальными значениями, если мы подставим 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;
}