Для изучения этой статьи вам потребуется базовое понимание методов обнаружения столкновений в узкой фазе, а также знание смежных с данной проблемой геометрических концепций, в частности, суммы Минковского.  

Несколько лет назад я посмотрел отличную презентацию от Дирка. В ней он описывал теорему о разделяющей оси, пролегающей между выпуклыми многогранниками (видео, слайды). Примерно на 18 минуте (слайд 29) он заводит речь о наложении гауссовых отображений выпуклых многогранников — как они помогают найти грани разности Минковского для этих многогранников.

Суть в том, что все грани разности Минковского являются либо:

1.    Гранью формы A

2.    Гранью формы B

3.    Векторным произведением ребра от A и ребра от B

В данном контексте наиболее важно, что в случае №3 можно обойтись без дорогостоящего полного вычисления опорной функции, поскольку известна вершина, расположенная на конкретной грани. Переварив этот факт, я задумался — почему бы не пропускать вычисление опорной функции для всех точек в гауссовом отображении? Оказывается, так в самом деле можно поступать, и в результате получается своего рода проверка с использованием разделяющей оси. В таком случае требуется всего один раз полностью вычислять опорную функцию, не утруждая себя FaceCountA+FaceCountB. Прежде, чем вдаваться в детали, давайте немного обрисуем контекст.

TLDR: Репозиторий на Github.

Формализуем обнаружение столкновений

Пусть A и B — это два ограниченных выпуклых множества, принадлежащих R³. Далее возможны два сценария:

1.    Пересечений нет — нас интересует точка, ближайшая между двумя множествами.

2.    Пересечение — требуется найти минимальный перенос, применив который можно разделить фигуры.

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

Пересечений нет

Допустим, A и B — непересекающиеся множества. Тогда берётся пара ближайших точек (по одной из каждого множества), такая, для которых расстояние между этими множествами является минимальным. Следующие выражения эквивалентны:

1) 
   minx,y |x-y|^2 
   x ∈ A
   y ∈ B

2) 
 = minx,y d   
   d ≥ |x-y|^2
   x ∈ A
   y ∈ B

3)
 = mind,z d       
   d ≥ |z|^2
   z ∈ A-

Обратите внимание: в последней формулировке проявляется разность Минковского. Поскольку все эти выражения выпуклые, выпуклые итерационные решатели должны быстро сходиться к глобальному решению.  Один из наиболее популярных таких решателей — это GJK, но, в принципе, сработает любой алгоритм для выпуклой оптимизации. Кроме того, итерационные методы можно запускать с применением тёплого старта – например, GJK использует найденные ранее решения. Существует великое множество выпуклых итерационных решателей, в качестве примеров можно привести методы внутренней точки, методы градиентного спуска, метод эллипсоидов, т.д.

Разговор об использовании иных решателей (кроме GJK) в играх — это тема для отдельного поста. В оставшейся части этого я сосредоточусь на случае с имеющимся пересечением.

Случай с пересечением

Допустим, A и B пересекаются. Чтобы размежевать эти множества, необходимо найти такие направление v и расстояние s, чтобы полученное переносом множество (B+s*v) не пересекалось с A.

5) 
   mins,v s
   |v| = 1
   (y+s*v) • v ≥ x • v, ∀x ∈ A, ∀y ∈ B
    
6) 
 = mins,v s
   |v| = 1
   (y-x) • v ≥ -s, ∀x ∈ A, ∀y ∈ B

7) 
 = mins,v s
   |v| = 1
   z • v ≤ s, ∀z ∈ A-B

8) 
 = mins,v s
   |v| = 1
   supp(A-B, v) ≤ s

9) 
 = min|v|=1 supp(A-B, v)

Данная задача не является выпуклой, поскольку имеется квадратичное ограничение равенства |v|=1. Применительно к интересующим нас опорным функциям эту задачу можно классифицировать как квадратично ограниченную квадратичную программу (QCQP). С геометрической точки зрения мы в данном случае минимизируем функцию на поверхности сферы, словно ищем самую глубокую низменность на Земле.

Надеюсь, из приведённого графика понятно, что при применении итерационного решателя не гарантируется нахождение глобального минимума. Однако, функция f(v), ограниченная пределами сферы, становится тем предсказуемее, чем ближе начало координат оказывается к границе оболочки. Ниже показано поперечное сечение опорной функции единичного куба по мере растущего сдвига от начала координат. Обратите внимание: опорные функции постепенно начинают выглядеть выпукло, удаляясь от отметки 0 радиан.

Вот мы и сформулировали язык, на котором можно обсуждать данную задачу оптимизации. В оставшейся части статьи выведем именно такое решение и, в частности, вариант алгоритма для создания SAT-решателя.

Решение задачи пересечений

Выше нам удавалось сформулировать стоящие перед нами задачи, не прибегая к каким-либо специализированным геометрическим реалиям, касающимся разностей Минковского, опорных функций, гауссовых отображений и т.д. Естественно, разность Минковского и гауссово отображение попали в окончательную формулировку задачи, но я сознательно не собираюсь подробно на них останавливаться. В этой статье я хочу показать, как данная проблема может быть смоделирована как просто очередная задача на оптимизацию без чрезмерного применения зубодробительной геометрии. Такая стратегия не раз пригодилась мне при решении продвинутых математических задач в ходе разработки игр. Данную философию я почерпнул из отличной лекции Криса Геккера Алгоритм Лемке: молоток из вашего математического арсенала

Начинаем с наиболее абстрактной и общей формулировки задачи, после чего настойчиво оптимизируем наше решение именно для того конкретного случая, с которым имеем дело.

Здесь мы оптимизируем функцию на единичной сфере — выражение (9). Итак, какими же свойствами нашей целевой функции наиболее удобно воспользоваться? Объективная функция имеет вид  

f(v) = maxz∈C (z • v) где C = A-B

Вот какими вопросами о ней можно задаться:

• Поддаётся ли f дифференцированию?

• Насколько быстро вычисляется f?

• Если мы только что вычислили f(v), то можем ли быстро вычислить f(v+e) для небольшого e?

• Т.д.

В зависимости от того, как мы ответим на эти вопросы, оптимальный решатель всегда будет выглядеть по-своему. В данном случае это зависит от формы множества C. Например, если C — это шар, то наша задача минимизации поддаётся аналитическому решению. В самом начале этой статьи я затронул тему выпуклых многогранников, так что теперь мы, наконец-то, подходим к сути.

Задача пересечения выпуклых многогранников

В оставшейся части данной статьи будем исходить из того, что C — это выпуклый многогранник, определяемый как выпуклая оболочка из линейно независимых вершин c₁, c₂, … c_n.

Пример в 2D

Приступая к делу, построим график опорной функции для C, определяемого следующими вершинами

c1 = [-1.2, 1.8]
c2 = [ 2,   2.8]
c3 = [ 1,  -0.7]
c4 = [-1,  -0.2]

Ниже приведены графики многоугольника, график опорной функции в R² и декартов график опорной функции на круге от 0 до 2π.

Из двух последних графиков без труда узнаём несколько полезных свойств f(v):

1.    Функция дифференцируема везде кроме конечного множества тех точек, в которых расположены изломы — разрывы в первой производной.

2.    Каждая дифференцируемая область соответствует одной из вершин нашего многоугольника и в пределах этого региона функция является вогнутой. Назовём эти области сферы вершинными. Помним, что C = A - B, так что область каждой вершины на самом деле соответствует наложенным областям двух вершин, одной от A и одной от B.

3.    Из определения опорной функции следует, что изломы возникают, например, когда
c₁ • v = c₂ • v. Это подразумевает, что они возникают в точке на сфере, удовлетворяющей  
(c₁ - c₂) • v = 0. Таким образом, v непременно должна быть перпендикулярна ребру многоугольника. Иными словами, изломы возникают на нормалях граней C. Это также интуитивно понятно: вообразите, что происходит, если опорное направление является заметающей от одной вершины до другой. В середине эта линия будет перпендикулярна грани многоугольника, и проекция вдоль этого направления будет одинакова для обеих вершин на этой грани.

Если вы стремитесь к максимальной строгости, то эти свойства можно доказать, исходя из определения опорной функции. Математический анализ подсказывает, что минимум функции f(θ), ограниченной областью вершины, должен находиться на границе этой области. В таком случае, глобально минимум должен находиться в одной из точек излома — и так мы получаем глобальный решатель. Мы просто вычисляем f во всех точках излома и выбираем минимальное значение из всех полученных. В данном примере перечислить изломы просто, так как мы указали все вершины, из которых можно произвести нормали граней. На практике мы знаем вершины A и B, но C в явном виде не строим. Однако, это не проблема, поскольку supp(A-B, v) = supp(A, v) - supp(B,v) и точки излома для суммы функций представляют собой объединение точек излома от отдельных функций. Поэтому на практике мы будем вычислять f во всех нормалях A и нормалях B, а затем выбирать минимальное значение. Надеюсь, вы уже догадываетесь, что мы открыли регулярный метод проверки столкновений при помощи разделяющей оси в 2D. При этом, мы не слишком вдавались в геометрию, а пользовались лишь базовыми фактами из математического анализа. Этот случай для 2D очень похож на случай для 3D, который мы и рассмотрим далее.

Общий случай в 3D

В R² разрывы производной, взятой от опорной функции, являются точками и расположены вдоль линий, проходящих через начало координат и пересечения этих линий с окружностью. Неудивительно, что в R³  разрывы опорной функции – это плоскости, проходящие через начало координат, а их пересечения со сферой являются большими кругами.

Если построить график разрывов, взятых одновременно от A и от B, то получится примерно такая картина, какую я привёл в начале этой статьи.

В случае для 3D наша целевая функция похожа по свойствам на аналогичную функцию для 2D:

1.    Вершинная область для c_i ограничена K_i плоскостей. Плоскость, разделяющая две соседние вершинные области, соответствующие c₁ и c₂, определяется как
(c₁ - c₂) • v = 0.

2.    На граничной кривой для вершинной области есть K_i разрывов, которые расположены на нормалях граней. Именно те рассуждения, которыми мы пользовались в 2D, можно обобщить до 3D. Поскольку нормали граней являются пересечениями геодезических кривых, имеем, что нормали граней ортогональны векторному произведению любых двух рёбер на грани.

3.    Минимум функции f в вершинной области находится где-то на нормалях граней. Это можно доказать средствами математического анализа, примерно как и выше. Минимум для участка в 2D находится на границе этого участка. Минимум для участка в 1D находится на его границе.

4.    Можно найти глобальный минимум, перечислив все нормали граней (изломы линий, ограничивающих вершинную область).

Из последнего свойства выводится проверка с применением разделяющей оси, описанная в презентации Дирка. А именно: мы перебираем все нормали граней в разности Минковского и вычисляем проекцию. Дирк указывает одну важную оптимизацию: новые нормали граней, образующиеся в результате векторных произведений, не требуют полного вычисления опорной функции, так как уже известна вершина, лежащая на границе разности Минковского. Можно развить эту идею и обойтись без вычисления опорной функции для всех нормалей граней (кроме одной), поскольку вполне ясно, что мы можем определить соседние вершинные области для нормали каждой грани. 

Оптимизированная SAT-проверка

Теперь к главной идее. Обходя дуги на сфере по соответствующему графовому алгоритму, ведём учёт опорных точек для каждой из областей. Если подробнее:

1.    Выбираем точку w на одном из концов дуги на рис. 7 и вычисляем ту вершинную область, в которой оказываемся в A и -B. То есть, вычисляем supp(A,w) и supp(-B,w) учитываем вершины a, b.

2.    Прослеживая дуги при движении по сфере следим, какие дуги при этом пересекаем. Всякий раз, пересекая дугу, попадаем в новую вершинную область.

3.    Каждое пересечение при таком обходе или его конечная точка w — это нормаль грани, и опорная функция — это просто w • (a - b).

Рассмотрим конкретный пример, воспользовавшись приведённой ниже схемой. Обозначим нормали граней греческими буквами α_i, β_i и γ_i, соответствующими нормалям граней A, -B и новым нормалям граней A-B, соответственно. Аналогично, обозначим латинскими буквами обозначим вершинные области A и -B. Таким образом, если G — это пересечение областей, обозначенных a₁ и b₂ , то

10)
    supp(A-B, v) = (a1 - b2) • v 
    для всех v в G

Вот несколько первых итераций данного примера

1.    Начинаем с произвольной точки α₁  и определяем вершинную область b₁, при этом полностью вычисляя опорную функцию B. Также мы «в подарок» узнаём, что находимся в a₁. Запишем:
f(α₁) = α₁ • (a₁ - b₁).

2.    Начинаем обход дуги от α₁ до α₂ и находим, что пересекаем дугу между β₁ и β₃ . В таком случае, точка пересечения —γ₄. Запишем: f(γ₄) = γ₄ • (a₁ - b₁).

3.    Когда на предыдущем шаге мы пересекали плоскость, мы перешли из области b₁ в область b₂. Продолжаем движение от γ₄ к γ₃ и запишем f(γ₃) = γ₄ • (a₁ - b₂).

4.    Продолжаем движение от γ₃ к α₂ . Поскольку мы пересекли другую плоскость, определяемую дугой между β₂  и β₃, обновим нашу вершинную область до b₃ . Запишем: f(α₂) = γ₄ • (a₁ - b₃).

5.    Продолжаем этот процесс до тех пор, пока не обойдём все дуги A.

Вот, в принципе, и всё. Здесь есть много важных деталей реализации — например, какими структурами данных пользоваться, а также как реализовать прослеживание дуг и проверку пересечений, но эта статья и так получилась слишком длинной. Если вас интересуют все эти тонкости, то можете посмотреть код в репозитории, ссылку на который я дал в самом начале статьи. Широкими мазками я уже обрисовал всё, что нужно.   

Воспользуйтесь каким-либо вариантом структуры данных, именуемой полуребро. Вам потребуется вся топологическая информация, которую только сможете добыть. Для каждой вершины я сохраняю двумерный массив всех её связных вершин. Для каждого ребра также сохраняю каждую из граней, расположенных по бокам от него.

struct Edge 
{ 
    uint8_t vert0, vert1;
    uint8_t face0, face1; 
};

Топологически отсортируем рёбра оболочки, ориентируясь на связность дуг, так, чтобы при обходе было точно известно, что вы успеваете начать обход каждой новой дуги, а обход этой дуги также завершаете именно на ней.

Убедитесь, что ваш построитель оболочек даёт качественные формы. Как и при обычной проверке с применением разделяющей оси, данный алгоритм чувствителен к качеству оболочки. Правда, я до сих пор не заметил никаких проблем с числами, так что он ничем не хуже обычного SAT.

Заключение

В этой статье мы смогли под необычным углом рассмотреть проверку с применением разделяющей оси — не как строго геометрический алгоритм, а как задачу оптимизации на сфере. Анализируя свойства опорной функции, мы обнаружили, что её минимум должен находиться в одном из «изломов», образующихся на пересечениях больших кругов с гауссовым отображением.

В данном случае важнее всего понять, что данное отображение — это просто граф, который можно обходить шаг за шагом. Можно не вычислять опорную функцию заново на каждом узле (нормали грани), а выполнить одно исходное вычисление, и затем, обходя граф, дёшево обновляя каждую из вершин при пересечении границ (дуг) и попадании в новые области. В сущности, это именно такая оптимизация опорной функции, которая во многих физических движках именуется «поиск восхождением к вершине» (hill climbing).  

Если хотите посмотреть, как всё это работает — я собрал для вас демо в Visual Studio для Windows, где можно сравнить два метода и визуализировать гауссовы отображения. В ходе моей проверки способ с обходом графа получался в 5-10 раз быстрее, чем обычный SAT для выпуклых оболочек с большим количеством граней. При сравнении оболочки с треугольником скорость возрастает примерно в 1,2 раза, так как здесь меньше векторных произведении рёбер. Кроме того, в той версии алгоритма, что реализована в демо, обнаружились некоторые проблемы при вычислении треугольников. Не уверен, с чем именно они связаны — с генерацией контактов или с обнаружением столкновений. Эта проблема определённо исправима, и в какой-то момент я получил рабочую модель, но эта история теряется в песках времени. Надеюсь, когда-нибудь найду время, чтобы снова её наладить.

Этот метод обхода я обнаружил независимо, но, как и в случае с большинством задач вычислительной геометрии, вероятно, кто-то решил эту задачу ещё в 1970-е. Если у этого алгоритма есть название, или вы знаете какие-то игры, в которых он используется — напишите об этом в комментариях!

Github