Mif 15

Задача. Вычислить расстояние между двумя отрезками [latex]AB[/latex] и [latex]CD[/latex], заданных координатами вершин в четырехмерном пространстве.

Тесты

(1)

[latex]A_0[/latex] [latex]A_1 [/latex] [latex]A_2[/latex] [latex]A_3[/latex] [latex]B_0[/latex] [latex]B_1[/latex] [latex]B_2[/latex] [latex]B_3[/latex]
1  1 0 0 0 2 0 0 0
2  1 0 0 0 3 0 0 0
3  -1 1 0 6 -2 -1 1 6
4  1 1 1 1 2 2 2 2
5  1 1 1 2 8 5 7 10
6  1 1 0 0 1 0 1 1
7  3 4 7 8 9 5 6 2
8  1 2 2 3 1 4 8 9

(2)

[latex]C_0[/latex] [latex]C_1[/latex] [latex]C_2[/latex] [latex]C_3[/latex] [latex]D_0[/latex] [latex]D_1[/latex] [latex]D_2[/latex] [latex]D_3[/latex] [latex]r[/latex]
1  1 1 0 0 2 1 0 0 1.000000
2  2  4  0  0 2 4 3 0 4.000000
3  -1 -1  0 6 2 1 1 6 1.154701
4  -9 -9 -9 -9 -5 -5 -5 -5 12.000000
5  8 0  0 0 2 1 1 1 1.414214
6  2  3 2 0 1 4 9 1 3.000000
7  7 4 11 15 15 5 12 9 9.000000
8  5 7 2 23 4 8 8 21 13.000000

Код программы 

Алгоритм и его обоснование

Расстояние между отрезками в четырехмерном пространстве находится по-разному, в зависимости от взаимного расположения этих отрезков. Тут мы можем выделить два основных случая:

  1. Отрезки лежат на параллельных прямых или на одной прямой.
  2. Отрезки лежат на пересекающихся либо на скрещивающихся прямых.

Чтобы выяснить, с каким случаем мы имеем дело, рассмотрим общую картину взаимного расположения отрезков и опишем ее математически:
PICTURE1
По условию нам заданы 4 точки: [latex]A[/latex], [latex]B[/latex], [latex]C[/latex] и [latex]D[/latex] — концы двух отрезков. Для удобства представления уравнений и точек, связанных с ними, обозначим их [latex]P_0[/latex], [latex]P_1[/latex], [latex]Q_0[/latex] и [latex]Q_1[/latex] соответственно. Через эти пары точек мы можем провести 2 прямые [latex]p[/latex] и [latex]q[/latex], параметрические уравнения которых имеют вид:

[latex]

\begin{matrix}
\vec p = P_0 + \vec u \cdot s \\
\vec q = Q_0 + \vec v \cdot t
\end{matrix}

[/latex],

где векторы:

[latex]

\begin{matrix}

\vec u = P_1 — P_0\\

\vec v = Q_1 — Q_0

\end{matrix}

[/latex],

а   [latex]s[/latex]   и   [latex]t[/latex]   — параметры. При [latex]s=0[/latex]   или   [latex] t=0[/latex]   мы получаем начальную точку соответствующего отрезка, а при [latex]s=1[/latex]   или [latex]t=1[/latex]   — конечную. При произвольном значении параметра мы получаем произвольную точку на прямой.

Рассмотрим вектор [latex]\vec w = Q — P[/latex] , соединяющий 2 произвольные точки на этих прямых. Легко показать, что вектор [latex]\vec w[/latex]   соединяет 2 ближайшие точки  [latex]Q_c[/latex]   и   [latex]P_c[/latex]   при условии:

[latex]\vec w \perp p[/latex] и [latex]\vec w \perp q[/latex].

Этому условию соответствует система из двух уравнений:

[latex] \begin{cases}
\vec u \cdot \vec w = 0\\
\vec v \cdot \vec w = 0
\end{cases}

[/latex]

Распишем ее для   [latex]\vec w = Q_0 — P_0 + \vec v \cdot t — \vec u \cdot s = \vec w_0 + \vec v \cdot t — \vec u \cdot s[/latex] :
[latex]

\begin{cases}
\vec u \cdot ( \vec w_0 + \vec v \cdot t — \vec u \cdot s ) = 0\\
\vec v \cdot ( \vec w_0 + \vec v \cdot t — \vec u \cdot s ) = 0
\end{cases}

[/latex]

Введем вспомогательные скалярные переменные:
[latex]

\begin{matrix}

a&=&\vec u \cdot \vec u\\
b&=&\vec u \cdot \vec v\\
c&=&\vec v \cdot \vec v\\
d&=&\vec u \cdot \vec w_0\\
e&=&\vec v \cdot \vec w_0

\end{matrix}

[/latex]

Теперь наша система будет выглядеть так:

[latex]

\begin{cases}
d — a \cdot s + b \cdot t = 0 \\
e — b \cdot s + c \cdot t = 0
\end{cases}
[/latex]

Перепишем систему в удобном для нас виде:

[latex] \begin{cases}
a \cdot s — b \cdot t = d \\
b \cdot s — c \cdot t = e
\end{cases}
[/latex]

Решение этой системы мы можем получить, например, методом Крамера.

Главный определитель системы:   [latex]D = b^2 — a \cdot c[/latex]

Два вспомогательных определителя:
[latex] \begin{matrix}
D_1 = b \cdot e — c \cdot d\\
D_2 = a \cdot e — b \cdot d\\
\end{matrix}
[/latex] Если [latex]D \neq 0[/latex],   то существует единственное решение:

[latex]

\begin{cases}
s_c = \frac{D_1}{D} \\
t_c = \frac{D_2}{D}
\end{cases}
[/latex]

Если же мы получаем   [latex]D = 0[/latex],   легко показать, что отрезки параллельны. То есть мы имеем дело со случаем 1.

Тогда:

а) Если хотя бы одна точка одного отрезка проецируется на другой отрезок, то расстояние между отрезками равняется расстоянию между прямыми.

Найдем проекцию точки   [latex]P_0[/latex]   на линию   [latex]q[/latex]. Для этого сначала найдем вектор, который является проекцией вектора   [latex]\vec w_0[/latex]   на линию   [latex]q[/latex].

[latex]\vec w_q=(\vec w_0 \cdot \vec v) \cdot \frac{\vec v}{v^2}[/latex].
Конец полученного вектора находится в точке   [latex]Q_0[/latex],   а начало в новой точке   [latex]P_{0q}=Q_0-\vec w_q[/latex]. Соединим точки   [latex]P_0[/latex]   и   [latex]P_{0q}[/latex] вектором [latex]\vec w_p = P_{0q} — P_0[/latex]. Длина полученного вектора и будет искомым расстоянием:   [latex]r = \left| P_0 P_{0q} \right|[/latex].

RESULT

Для проверки условия а) необходимо получить проекции остальных исходных точек на отрезки:
[latex] \begin{matrix}
P_{1q} = P_{0q} + \vec u\\
Q_{0p} = P_0 + \vec w_q\\
Q_{1p} = Q_{0p} + \vec v
\end{matrix}
[/latex]

Если точка   [latex]P_{0q}[/latex]   лежит на прямой   [latex]q[/latex],    задаваемой уравнением:
[latex]\vec q = Q_0 + \vec v \cdot t[/latex],
то определить, принадлежит ли точка [latex]P_{0q}[/latex]   отрезку [latex]Q_0 Q_1[/latex]   можно, решив уравнение:
[latex]P_{0q} = Q_0 + \vec v \cdot t[/latex].
[latex] 0 = \vec w_q + \vec v \cdot t[/latex].
Домножив обе части скалярно на вектор   [latex]\vec v[/latex],   мы получим уравнение: [latex] 0 = e + c \cdot t[/latex], отсюда   [latex]t = \frac{-e}{c}[/latex].

Если [latex]t \in \left[0,1\right][/latex], то точка [latex]P_{0q}[/latex] лежит на отрезке [latex]Q_0 Q_1[/latex]. Если же нет, переходим к аналогичной проверке следующих точек:

[latex]P_{1q}:[/latex] [latex]P_{1q} = Q_0 + \vec v \cdot t[/latex].
[latex] 0 = \vec w_q — \vec u + \vec v \cdot t[/latex].
Опять домножив обе части скалярно на вектор   [latex]\vec v[/latex],   мы получим уравнение:

[latex] 0 = e — b + c \cdot t[/latex],
отсюда   [latex]t = \frac{-e+b}{c}[/latex].
[latex]Q_{0p}:[/latex] [latex]Q_{0p} = P_0 + \vec u \cdot s[/latex].
[latex]0 =-\vec w_q + \vec u \cdot s[/latex].
Опять домножив обе части скалярно на вектор   [latex]\vec u[/latex],   мы получим уравнение:

[latex] 0 = -\frac{e \cdot b}{c} + a \cdot s[/latex],
отсюда   [latex]s = \frac{e \cdot b}{c \cdot a}[/latex].
[latex]Q_{1p}:[/latex] [latex]Q_{1p} = P_0 + \vec u \cdot s[/latex].
[latex] 0 = -\vec w_q — \vec v + \vec u \cdot s[/latex].
Опять домножив обе части скалярно на вектор   [latex]\vec u[/latex],   мы получим уравнение:

[latex] 0 = -\frac{e \cdot b}{c} — b + a \cdot s[/latex],
отсюда   [latex]s = \frac{(e — c) \cdot b}{c \cdot a}[/latex].
б) В противном случае, расстояние между отрезками равняется минимальному расстоянию между их концами. Здесь задача предельно упрощается. Мы находим длины отрезков, попарно соединяющих 4 исходные точки, и выбираем наименьший из них.

Если же исходные отрезки лежат на пересекающихся либо на скрещивающихся прямых, мы также рассматриваем 2 случая:

а) Оба конца кратчайшего отрезка, соединяющего прямые, лежат на соответствующих исходных отрезках:
[latex]P_c \in P_0 P_1[/latex]   и   [latex]Q_c \in Q_0 Q_1[/latex].

В этом случае пара параметров   [latex](s_c, \; t_c)[/latex]   принадлежит области:   [latex](s,t):\left[0,1\right]\times \left[0,1\right].[/latex]

То есть, решение тривиально: ответом будет дина вектора   [latex]\vec w_c[/latex]

б) Хотя бы один из концов кратчайшего отрезка, соединяющего прямые, не лежит на исходном отрезке, то есть:
[latex]P_c \not\in P_0 P_1[/latex] или [latex]Q_c \not\in Q_0 Q_1[/latex],
что соответствует значениям параметров   [latex]s_c \not\in \left[0,1\right][/latex]   или   [latex]t_c \not\in \left[0,1\right][/latex].

В этом случае минимальное расстояние между отрезками определяется на границе области:   [latex](s,t):\left[0,1\right]\times \left[0,1\right][/latex]   (см. рисунок ниже):

elipsoid

Здесь решением является длина кратчайшего отрезка.

Длину отрезка, соединяющего 2 прямые, можно оценивать по квадрату длины вектора   [latex]\vec v[/latex]: [latex]w^2=(\vec w)^2=(\vec w_0 — \vec u \cdot s + \vec v \cdot t)^2[/latex].

В частности, минимум   [latex]w^2[/latex] достигается в точке   [latex](s_c,t_c)[/latex].
Однако в случае б) мы должны найти минимум расстояния на границе нашей области, то есть решить задачу нахождения минимума при ограничениях (решить задачу условной минимизации). В нашем случае ограничения имеют очень простой вид — оси координат, и две линии, параллельные им. Поэтому мы можем решить на четырех границах 4 упрощенные задачи минимизации, а затем выбрать наименьшее решение.

Замечание: В пространстве параметров функция [latex]w^2(s,t)[/latex] представляет из себя эллиптический параболоид. Однако для простоты мы выше изобразили его линии уровня в виде окружностей. Типичный вид эллиптического параболоида и его линий уровня представлен на рисунках ниже:
3dellipticparabolloid
2dmap_ellipticparabolloid

Рассмотрим поочередно все 4 ограничения и решим задачу для них:

(1) Пусть [latex]t=t_1=0[/latex].

Тогда: [latex]{w^2\mid_{t_1=0}} = (\vec w_0-\vec u \cdot s_1)^2[/latex].

Для определения экстремума приравняем производную к нулю:[latex] \begin{array}{r}
\frac{d}{ds_1}{(\vec w_0-\vec u \cdot s_1)^2}=0\\
2 \cdot (\vec w_0-\vec u \cdot s_1) \cdot (- \vec u)=0\\
-d +a \cdot s_1=0\\
s_1=\frac{d}{a}
\end{array}
[/latex]

Легко показать, что при [latex]s_1>1[/latex] мы должны присвоить ему значение [latex]s_1=1[/latex], а если [latex]s<0[/latex] — значение [latex]s_1=0[/latex], так как мы не должны выходить за границы исходных отрезков.

Подставим полученное значение [latex]s[/latex] в уравнение прямой [latex]p[/latex] для точки [latex]P_c[/latex]:[latex]P_c = P_0 + \vec u \cdot s.[/latex]

А точка [latex]Q_c[/latex] совпадает с точкой [latex]Q_0[/latex]. Тогда первый минимум равен: [latex]r_1 = Q_0 P_c[/latex].

Аналогично найдем три остальных минимума [latex]r_2, r_3, r_4[/latex], приравняв [latex]s[/latex] к нулю, а затем [latex]t[/latex] и [latex]s[/latex] к единице. Наименьший из них и есть искомое расстояние [latex]r[/latex].

Код программы

4 thoughts on “Mif 15

  1. Молодец.
    Большая серьёзная работа. Признаться, я предполагал что будут записаны параметрические уравнения обоих отрезков, далее вычисляется квадрат расстояния между точками отрезков как функция обоих параметров. Останется только отыскать нули производной и найти минимум внутри единичного квадрата.

    Давайте пока я внимательно всё перепроверю в выкладках, Вы поправите отступы в программе?

    • Да, спасибо, исправила. Я использовала стиль GNU с небольшим дополнением, для удобства чтения и анализа кода: если одна операция зависит от другой по данным, он отделяется отступом в один пробел.

  2. — В 53-й строке else закомментарен сознательно?

    Что касается отступов.
    Программирование обычно коллективная работа. Если программист-одиночка пишет код для себя, то самовыражение через отступы вполне уместно. Но это редкость. Обычно отступы несут существенную коммуникативную нагрузку. Использовать «более правильные» отступы в коде, который будут читать люди так же нелепо, как использовать «более правильные» слова в языке на котором общаешься.

    Что касается данной модификации отступов.
    Она имеет существенные недостатки. Например, отступ в 1 пробел делается после условного оператора или цикла даже если последующий оператор не зависит по данным от условного, а только может зависеть. Т.е. только во время выполнения программы станет достоверно известно, нужен ли был дополнительный отступ в коде. Приходится его ставить всегда на всякий случай.
    Назначение отступов (индент табуляции) — демонстрации иерархической структуры программы. Модификация (1 пробел) — демонстрация зависимости по данным. Т.е. один механизм для разных вещей. При нескольких отступах второго типа они неотличимы от отступов первого типа.

    Признаться, мне нравится это дополнение.
    Я склонен проявить малодушие и разрешить Вам пользоваться такими отступами. Хоть и обоснованно предполагаю, что Вам это может повредить в дальнейшем при работе в команде.
    Но Вы ведь умный человек, как-то выкрутитесь.
    Или нет.

Добавить комментарий