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

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

В этой части туториала мы рассмотрим следующие темы:

• Простое распознавание коллизий
• Генерирование простого многообразия
• Разрешение импульсов силы

Вот небольшое демо:


Примечание: хоть этот туториал написан на C++, вы сможете использовать те же техники и концепции почти в любой среде разработки игр.

---

Необходимые знания

Для понимания этой статьи требуется неплохое знание математики и геометрии, и в гораздо меньшей степени — непосредственно программирования. В частности, вам потребуется следующее:

• Базовое понимание основ векторной математики
• Умение выполнять алгебраические вычисления

---

Распознавание коллизий

В Интернете достаточно статей и туториалов о распознавании коллизий, поэтому я не буду подробно рассматривать эту тему.

Ограничивающий прямоугольник, выровненный по координатным осям

Ограничивающий прямоугольник, выровненный по координатным осям (Axis Aligned Bounding Box, AABB) — это прямоугольник, четыре оси которого выровнены относительно системы координат, в которой он находится. Это значит, что прямоугольник не может вращаться и всегда находится под углом в 90 градусов (обычно выровнен относительно экрана). Обычно его называют «ограничивающим прямоугольником», потому что AABB используются для ограничения других, более сложных форм.


Пример AABB.

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

struct AABB
{
  Vec2 min;
  Vec2 max;
};

Эта форма позволяет задать AABB двумя точками. Точка min обозначает нижние границы по осям x и y, а max обозначает верхние границы — иными словами, они обозначают верхний левый и нижний правый углы. Чтобы определить, пересекаются ли два AABB, необходимо базовое понимание теоремы о разделяющей оси (Separating Axis Theorem, SAT).

Вот быстрая проверка, взятая с сайта Real-Time Collision Detection Кристера Эриксона, в которой используется SAT:

bool AABBvsAABB( AABB a, AABB b )
{
  // Выходим без пересечения, потому что найдена разделяющая ось
  if(a.max.x < b.min.x or a.min.x > b.max.x) return false
  if(a.max.y < b.min.y or a.min.y > b.max.y) return false

  // Разделяющая ось не найдена, поэтому существует по крайней мере одна пересекающая ось
  return true
}

Окружности

Окружность задаётся радиусом и точкой. Вот как может выглядеть структура окружности:

struct Circle
{
  float radius
  Vec position
};

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

Важная оптимизация, позволяющая избавиться от оператора квадратного корня:

float Distance( Vec2 a, Vec2 b )
{
  return sqrt( (a.x - b.x)^2 + (a.y - b.y)^2 )
}

bool CirclevsCircleUnoptimized( Circle a, Circle b )
{
  float r = a.radius + b.radius
  return r < Distance( a.position, b.position )
}

bool CirclevsCircleOptimized( Circle a, Circle b )
{
  float r = a.radius + b.radius
  r *= r
  return r < (a.x + b.x)^2 + (a.y + b.y)^2
}

В общем случае умножение — это гораздо менее затратная операция, чем получение квадратного корня значения.

---

Разрешение импульсов силы

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

В общем случае объект в физическом движке имеет три основные степени свободы (в двух измерениях): движение в плоскости xy и вращение. В этой статье мы намеренно ограничили вращение и используем только AABB с окружностями, поэтому единственная степень свободы, которую нам нужно будет рассматривать — это движение в плоскости xy.

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


Простой пример того, чего можно достичь с помощью разрешения импульсов силы.

Чтобы достигнуть такого эффекта и при этом следовать интуитивному пониманию того, как должны вести себя объекты, мы используем твёрдые тела и немного математики. Твёрдое тело — это просто форма, задаваемая пользователем (то есть вами, разработчиком), которая явно определяется как недеформируемая. И AABB, и окружности в этой статье недеформируемы, и всегда будут являться либо AABB, либо окружностью. Все сжатия и растяжения запрещены.

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

Объекты столкнулись — что дальше?

Предположим, что мы обнаружили столкновение двух тел. Как их разделить? Будем считать, что распознавание коллизий предоставляет нам две важные характеристики:

• Нормаль коллизии
• Глубина проникновения

Чтобы применить импульс силы к обоим объектам и оттолкнуть их друг от друга, нам нужно знать, в каком направлении и насколько их отталкивать. Нормаль коллизии — это направление, в котором будет приложен импульс силы. Глубина проникновения (вместе с некоторыми другими параметрами) определяет, насколько большим будет используемый импульс силы. Это значит, что единственное значение, которое нам нужно вычислить — это величина импульса силы.

Теперь давайте подробно рассмотрим, как же вычислить величину импульса силы. Начнём с двух объектов, для которых обнаружено пересечение:

Уравнение 1

$$

Заметьте, что для создания вектора из положения A в положение B необходимо выполнить: endpoint - startpoint. $$ — это относительная скорость из A в B. Это уравнение можно выразить относительно нормали коллизии $$, то есть мы хотим узнать относительную скорость из A в B вдоль направления нормали коллизии:

Уравнение 2

$$

Теперь мы используем скалярное произведение. Скалярное произведение — это просто сумма покомпонентных произведений:

Уравнение 3

$$

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

Чтобы выбрать нужную упругость (обозначаемую как $$, «эпсилон»), отвечающую интуитивно ожидаемым результатам, нам следует использовать наименьшую задействованную упругость:

// Два заданных объекта A и B
e = min( A.restitution, B.restitution )

Получив $$, мы можем подставить его в уравнение вычисления величины импульса силы.

Ньютоновский закон восстановления гласит следующее:

Уравнение 4

$$

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

Уравнение 5

$$

Заметьте, что здесь появилось отрицательное значение. Notice how we introduced a negative sign here. По ньютоновскому закону восстановления $$, результирующий вектор после отталкивания, действительно направляется в обратную сторону от V. Так как же представить противоположные направления в нашем уравнении? Ввести знак «минус».

Теперь нам нужно выразить эти скорости под воздействием импульса силы. Вот простое уравнение для изменения вектора на скаляр импульса силы $$ в определённом направлении $$:

Уравнение 6

$$

Надеюсь, это уравнение вам понятно, потому что оно очень важно. У нас есть единичный вектор $$, обозначающий направление. Также у нас есть скаляр $$, обозначающий длину вектора $$. При суммировании отмасштабированного вектора $$ с $$ мы получаем $$. Это просто сложение двух векторов, и мы можем использовать это небольшое уравнение для приложения импульса силы одного вектора к другому.

Здесь нам ещё предстоит проделать небольшую работу. Формально импульс силы определяется как изменение импульса. Импульс — это масса * скорость. Зная это, мы можем выразить импульс в соответствии с формальным определением так:

Уравнение 7

$$

Три точки в форме треугольника ($$) можно прочитать как «следовательно». Это обозначение используется для того, чтобы из предшествующего ему вывести истинность последующего.

Мы неплохо двигаемся! Однако нам нужно выразить импульс силы с помощью $$ относительно двух разных объектов. Во время коллизии объектов A и B объект A отталкивается в противоположном от B направлении:

Уравнение 8

$$

Эти два уравнения отталкивают A от B вдоль единичного вектора направления $$ на скаляр импульса силы (величины $$) $$.

Всё это нужно для объединения уравнений 8 и 5. Конечное уравнение будет выглядеть примерно так:

Уравнение 9

$$

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

Уравнение 10

$$

Ого, довольно много вычислений! Но на этом всё. Важно понимать, что в окончательной форме уравнения 10 слева у нас $$ (величина), а всё справа нам уже известно. Это значит, что мы можем написать пару строк кода для вычисления скаляра импульса силы $$. И этот код гораздо более читаем, чем математическая запись!

void ResolveCollision( Object A, Object B )
{
  // Вычисляем относительную скорость
  Vec2 rv = B.velocity - A.velocity

  // Вычисляем относительную скорость относительно направления нормали
  float velAlongNormal = DotProduct( rv, normal )

  // Не выполняем вычислений, если скорости разделены
  if(velAlongNormal > 0)
    return;

  // Вычисляем упругость
  float e = min( A.restitution, B.restitution)

  // Вычисляем скаляр импульса силы
  float j = -(1 + e) * velAlongNormal
  j /= 1 / A.mass + 1 / B.mass

  // Прикладываем импульс силы
  Vec2 impulse = j * normal
  A.velocity -= 1 / A.mass * impulse
  B.velocity += 1 / B.mass * impulse
}

В этом примере кода нужно заметить два важных аспекта. Во-первых, посмотрите на строку 10, if(VelAlongNormal > 0). Эта проверка очень важна, она гарантирует, что мы разрешаем коллизию, только если объекты движутся друг к другу.


У двух объектов возникла коллизия, но скорость разделит их в следующем кадре. Не разрешаем этот тип коллизии.

Если объекты движутся в противоположные друг от друга стороны, мы ничего не делаем. Благодаря этому мы не будем разрешать коллизии тех объектов, которые на самом деле не сталкиваются. Это важно для создания симуляции, соответствующей интуитивным ожиданиям о том, что должно происходить при взаимодействии объектов.

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

A.inv_mass = 1 / A.mass

Во многих физических движках необработанная масса на самом деле не хранится. Часто физические движки хранят только величину, обратную массе. Просто так бывает, что в большинстве математических расчётов используется масса в виде 1/масса.

И последнее, что нужно заметить, что мы должны с умом распределить наш скаляр импульса силы $$ на два объекта. Мы хотим, чтобы мелкие объекты отлетали от крупных с большей долей $$, а скорости больших объектов изменялись на очень небольшую долю $$.

Для этого можно сделать следующее:

float mass_sum = A.mass + B.mass
float ratio = A.mass / mass_sum
A.velocity -= ratio * impulse

ratio = B.mass / mass_sum
B.velocity += ratio * impulse

Важно осознавать, что этот код аналогичен приведённому выше примеру функции ResolveCollision(). Как объяснялось выше, обратные массы довольно полезны в физическом движке.

Тонущие объекты

Если мы воспользуемся уже написанным кодом, то объекты будут сталкиваться и отлетать друг от друга. Это отлично, но что случится, если один из объектов имеет бесконечную массу? Нам потребуется удобный способ задания в нашей симуляции бесконечной массы.

Я предлагаю использовать в качестве бесконечной массы ноль — однако если мы попробуем вычислить обратную массу объекта с нулевой массой, мы получим деление на ноль. Решить эту проблему при вычислении обратной массы можно следующим образом:

if(A.mass == 0)
  A.inv_mass = 0
else
  A.inv_mass = 1 / A.mass

Значение «ноль» приведёт к верным вычислениям при разрешении импульсов силы. Это нас устраивает. Проблема тонущих объектов возникает, когда какой-нибудь объект начинает «тонуть» в другом из-за гравитации. Иногда объект с низкой упругостью ударяется о стену с бесконечной массой и начинает тонуть.

Такое утопание возникает из-за ошибок вычислений с плавающей запятой. Во время каждого вычисления с плавающей запятой добавляется небольшая ошибка из-за ограничений оборудования. (Подробнее см. [Floating point error IEEE754] в Google.) Со временем эта ошибка накапливается в ошибку позиционирования, что приводит к утоплению объектов друг в друге.

Для исправления этой ошибки позиционирования необходимо её учитывать, поэтому я покажу вам способ, называемый «линейным проецированием». Линейное проецирование на небольшой процент снижает проникновение двух объектов друг в друга. Оно выполняется после приложения импульса силы. Исправление положения выполняется очень просто: перемещаем каждый объект вдоль нормали коллизии $$ на процент глубины проникновения:

void PositionalCorrection( Object A, Object B )
{
  const float percent = 0.2 // обычно от 20% до 80%
  Vec2 correction = penetrationDepth / (A.inv_mass + B.inv_mass)) * percent * n
  A.position -= A.inv_mass * correction
  B.position += B.inv_mass * correction
}

Учтите, что мы масштабируем penetrationDepth на общую массу системы. Это даст нам коррекцию положения, пропорциональную величине массы. Мелкие объекты отталкиваются быстрее, чем тяжёлые.

Однако в этой реализации есть небольшая проблема: если мы всегда разрешаем ошибку позиционирования, то объекты всегда будут дрожать, пока они находятся друг на друге. Чтобы устранить дрожание, нужно задать небольшой допуск. Мы будем выполнять корректировку положения только если проникновение выше определённого произвольного порога, который мы назовём «погружением» («slop»):

void PositionalCorrection( Object A, Object B )
{
  const float percent = 0.2 // обычно от 20% до 80%
  const float slop = 0.01 // обычно от 0.01 до 0.1
  Vec2 correction = max( penetration - k_slop, 0.0f ) / (A.inv_mass + B.inv_mass)) * percent * n
  A.position -= A.inv_mass * correction
  B.position += B.inv_mass * correction
}

Это позволит объектам немного проникать друг в друга без задействования коррекции положения.

---

Генерирование простого многообразия

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

Вот как выглядит объявление стандартного многообразия:

struct Manifold
{
  Object *A;
  Object *B;
  float penetration;
  Vec2 normal;
};

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

Окружность-окружность

Давайте начнём с простейшего алгоритма коллизии: коллизия окружность-окружность. Эта проверка в большей степени тривиальна. Можете ли вы представить, каким будет направление разрешения коллизии? Это вектор от окружности A к окружности B. Его можно получить вычитанием положения B из положения A.

Глубина проникновения связана с радиусами окружностей и расстоянием между ними. Наложение окружностей можно вычислить вычитанием из суммы радиусов расстояния до каждого из объектов.

Вот полный пример алгоритма генерирования многообразия коллизии окружность-окружность:

bool CirclevsCircle( Manifold *m )
{
  // Объявление пары указателей на каждый объект
  Object *A = m->A;
  Object *B = m->B;

  // Вектор от A к B
  Vec2 n = B->pos - A->pos

  float r = A->radius + B->radius
  r *= r

  if(n.LengthSquared( ) > r)
    return false

  // У окружностей распознана коллизия, вычисляем многообразие
  float d = n.Length( ) // вычисляем sqrt

  // Если расстояние между окружностями не равно нулю
  if(d != 0)
  {
    // Расстояние - это разность между радиусом и расстоянием
    m->penetration = r - d

    // Используем d, потому что мы уже вычислили sqrt в Length( )
    // Направлен из A в B, и это единичный вектор
    c->normal = t / d
    return true
  }

  // Окружности имеют одинаковое положение
  else
  {
    // Выбираем случайные (но согласованные) значения
    c->penetration = A->radius
    c->normal = Vec( 1, 0 )
    return true
  }
}

Здесь стоит заметить следующее: мы не выполняем вычислений квадратного корня, пока без этого можно обойтись (если у объектов нет коллизии), и мы проверяем, не находятся ли окружности в одной точке. Если они находятся в одной точке, то расстояние будет равно нулю и нужно избежать деления на ноль при вычислении t / d.

AABB-AABB

Проверка AABB-AABB немного более сложна, чем окружность-окружность. Нормаль коллизии не будет вектором из A в B, а будет нормалью к ребру. AABB — это прямоугольник с четырьмя рёбрами. Каждое ребро имеет нормаль. Эта нормаль обозначает единичный вектор, перпендикулярный к ребру.

Исследуем общее уравнение прямой в 2D:

$$

В уравнении выше a и b — это вектор нормали к прямой, а вектор (a, b) считается нормализованным (длина вектора равна нулю). Нормаль коллизии (направление разрешения коллизии) будет направлена в сторону одной из нормалей рёбер.

Знаете ли вы, что обозначает c в общем уравнении прямой? c — это расстояния до начала координат. Как мы увидим в следующей части статьи, это очень полезно для проверки того, на какой стороне от прямой находится точка.

Всё, что теперь нужно — определить, какое из рёбер одного объекта сталкивается с другим объектом, после чего мы получим нормаль. Однако иногда могут пересекаться несколько рёбер двух AABB, например, при пересечении двух углов. Это значит, что нам нужно определить ось наименьшего проникновения.


Две оси проникновения; горизонтальная ось X — ось наименьшего проникновения, поэтому эту коллизию нужно разрешать вдоль оси X.

Вот полный алгоритм генерирования многообразия AABB-AABB и распознавания коллизий:

bool AABBvsAABB( Manifold *m )
{
  // Задание пары указателей для каждого из объектов
  Object *A = m->A
  Object *B = m->B
 
  // Вектор из A в B
  Vec2 n = B->pos - A->pos
 
  AABB abox = A->aabb
  AABB bbox = B->aabb
 
  // Вычисление половины ширины вдоль оси x для каждого объекта
  float a_extent = (abox.max.x - abox.min.x) / 2
  float b_extent = (bbox.max.x - bbox.min.x) / 2
 
  // Вычисление наложения по оси x
  float x_overlap = a_extent + b_extent - abs( n.x )
 
  // Проверка SAT по оси x
  if(x_overlap > 0)
  {
    // Вычисление половины ширины вдоль оси y для каждого объекта
    float a_extent = (abox.max.y - abox.min.y) / 2
    float b_extent = (bbox.max.y - bbox.min.y) / 2
 
    // Вычисление наложения по оси y
    float y_overlap = a_extent + b_extent - abs( n.y )
 
    // Проверка SAT по оси y
    if(y_overlap > 0)
    {
      // Определяем, по какой из осей проникновение наименьшее
      if(x_overlap > y_overlap)
      {
        // Указываем в направлении B, зная, что n указывает в направлении от A к B
        if(n.x < 0)
          m->normal = Vec2( -1, 0 )
        else
          m->normal = Vec2( 0, 0 )
        m->penetration = x_overlap
        return true
      }
      else
      {
        // Указываем в направлении B, зная, что n указывает в направлении от A к B
        if(n.y < 0)
          m->normal = Vec2( 0, -1 )
        else
          m->normal = Vec2( 0, 1 )
        m->penetration = y_overlap
        return true
      }
    }
  }
}

Окружность-AABB

Последняя проверка, которую я рассмотрю — проверка окружность-AABB. Идея здесь заключается в том, чтобы вычислить ближайшую к окружности точку AABB; после этого проверка упрощается до чего-то вроде проверки окружность-окружность. После вычисления ближайшей точки и распознавания коллизий нормаль будет направлением к ближайшей точке из центра окружности. Глубина проникновения — это разность между расстояниями до ближайшей к окружности точки и радиусом окружности.


Схема пересечения AABB-окружность.

Есть один хитрый особый случай — если центр окружности находится внутри AABB, то нужно центр окружности отсечь до ближайшего ребра AABB, а нормаль отразить.

bool AABBvsCircle( Manifold *m )
{
  // Задание пары указателей для каждого из объектов
  Object *A = m->A
  Object *B = m->B

  // Вектор от A к B
  Vec2 n = B->pos - A->pos

  // Ближайшая к центру B точка A
  Vec2 closest = n

  // Вычисление половины ширины вдоль каждой оси
  float x_extent = (A->aabb.max.x - A->aabb.min.x) / 2
  float y_extent = (A->aabb.max.y - A->aabb.min.y) / 2

  // Ограничиваем точку ребром AABB
  closest.x = Clamp( -x_extent, x_extent, closest.x )
  closest.y = Clamp( -y_extent, y_extent, closest.y )

  bool inside = false

  // Окружность внутри AABB, поэтому нам нужно ограничить центр окружности
  // до ближайшего ребра
  if(n == closest)
  {
    inside = true

    // Находим ближайшую ось
    if(abs( n.x ) > abs( n.y ))
    {
      // Отсекаем до ближайшей ширины
      if(closest.x > 0)
        closest.x = x_extent
      else
        closest.x = -x_extent
    }

    // ось y короче
    else
    {
      // Отсекаем до ближайшей ширины
      if(closest.y > 0)
        closest.y = y_extent
      else
        closest.y = -y_extent
    }
  }

  Vec2 normal = n - closest
  real d = normal.LengthSquared( )
  real r = B->radius

  // Если радиус меньше, чем расстояние до ближайшей точки и
  // Окружность не находится внутри AABB
  if(d > r * r && !inside)
    return false

  // Избегаем sqrt, пока он нам не понадобится
  d = sqrt( d )

  // Если окружность была внутри AABB, то нормаль коллизии нужно отобразить
  // в точку снаружи
  if(inside)
  {
    m->normal = -n
    m->penetration = r - d
  }
  else
  {
    m->normal = n
    m->penetration = r - d
  }

  return true
}

---

Заключение

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

• Сортировка и отсечение контактных пар
• Широкая фаза
• Расслоение
• Интеграция
• Такты
• Пересечение полупространств
• Модульность (материалы, масса и силы)