Нахождение расстояния между двумя отрезками в пространстве

Question

Необходимо найти кратчайшее расстояние между двумя отрезками AB и CD (координаты концов отрезков заданы в виде X,Y,Z) которые заданы в пространстве.

Я реализовал эту задачу программно перебором всех точек и нахождение самого минимального расстояния. Как это можно формульно реализовать, чтобы подставив координаты получить нужное решение? Ниже моя реализация

double DIST(double Ax, double Ay, double Az, double Bx, double By, double Bz, double Cx, double Cy, double Cz, double Dx, double Dy, double Dz){
double min = 100000;
for (double t=0; t<=1; t+=0.0005)
    for (double s = 0; s <= 1; s += 0.0005)
    {
        double EF = sqrt(pow(t * Ax + (1 - t) * Bx - s * Cx - (1 - s) * Dx, 2) + pow(t * Ay + (1 - t) * By - s * Cy - (1 - s) * Dy, 2) + pow(t * Az + (1 - t) * Bz - s * Cz - (1 - s) * Dz, 2));
        if (min > EF)
            min = EF;
    }
return min;}

Answer 1

Мне кажется, что проще всего -

составляем параметрическое уравнение для первого отрезка типа x(t) = x0 + t(x1-x0), y(t) = y0 + t(y1-y0) ... и такие же для второго x(u), y(u) ... t меняется от 0 до 1, u тоже.
Записываем квадрат расстояния между точками как функцию d(u,t) = (x(t)-x(u))² + ...
Ищем ее экстремум дифференцированием. Если он в пределах [0,1] и для t, и для u - ура, найдено. Если нет - надо смотреть предельные расстояния для концов отрезков.

Словом, типичная задача поиска минимума функции двух переменных.

Answer 2

Реализация на Delphi нахождения расстояния между скрещивающимися прямыми в пространстве (wiki) (С версия).

В Вашем случае (отрезки) нужно посмотреть, куда попадают концы общего перпендикуляра. Если параметры s и t выходят за пределы диапазона 0..1, то кратчайшее расстояние будет отсчитываться от соответствующего конца (концов) отрезка.

function LineLineDistance(const L0, L1: TLine3D; var s, t: Double): Double;
var
  u: TPoint3D;
  a, b, c, d, e, det, invdet:Double;
begin
  u := VecDiff(L1.Base, L0.Base);
  a := Dot(L0.Direction, L0.Direction);
  b := Dot(L0.Direction, L1.Direction);
  c := Dot(L1.Direction, L1.Direction);
  d := Dot(L0.Direction, u);
  e := Dot(L1.Direction, u);
  det := a * c - b * b;
  if det < eps then   
    Result := -1
  else begin
    invdet := 1 / det;
    s := invdet * (b * e - c * d);
    t := invdet * (a * e - b * d);

    Result := Distance(PointAtParam(L0, s), PointAtParam(L1, t));
  end;
end;

БЛОГ НА HUSL

Нахождение расстояния между двумя отрезками в пространстве

Ответы (2 шт):