Равномерное перемещение объекта вдоль кривой



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

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

В нашем случае сплайн — это функция, отображающая набор входных параметров (контрольные точки) и значение аргумента $t$ (обычно принимающего значения от 0 до 1) в точку на плоскости или в пространстве. Получившаяся кривая — это множество значений функции-сплайна при $t\in[0,1]$.

В качестве примера можно рассмотреть кубическую кривую Безье, которая задаётся следующим уравнением:
image


image
Кубическая кривая Безье

На рисунке изображены две кубических кривые Безье, задаваемых четырьмя точками (через две из них кривая проходит, оставшиеся две задают кривизну).

image
Анимация отображения параметра t в точку кривой

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


Кусочный сплайн

С заданием и параметризацией маршрута разобрались, теперь переходим к основному вопросу. Чтобы наш условный самолёт двигался вдоль маршрута с постоянной скоростью, нам нужно в любой момент времени уметь вычислять точку на кривой в зависимости от пройденного вдоль этой кривой расстояния $s(len)$, при этом имея лишь возможность вычислить положение точки на кривой по значению параметра $t$ (функцию $y(t)$). Именно на этом этапе начинаются сложности.

Первой мыслью было сделать линейное отображение $len\in[0, Length] \Rightarrow t\in[0,1]$ и вычислить $y(t)$ от получившегося значения — легко, вычислительно дёшево, в общем, то, что нужно.

image

Проблема данного способа очевидна сразу же — на самом деле пройденная дистанция $s$ зависит от $t$ нелинейно, и чтобы в этом убедиться достаточно расставить по маршруту равномерно распределенные по значению $t$ точки:


«Равномерно» распределенные по пути точки

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


Визуализация неправильной параметризации кривой

Обратившись за советом в поисковик, на stackoverflow и youtube, я обнаружил второй способ вычисления $s(len) = g(t)$, а именно представление кривой в виде кусочно-линейной функции (расчета набора равноудаленных друг от друга по кривой точек):


Представление кривой в виде кусочно-линейного сплайна

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

Пример кода
Источник

public Vector2[] CalculateEvenlySpacedPoints(float spacing, float resolution = 1)
  {
    List<Vector2> evenlySpacedPoints = new List<Vector2>();
    evenlySpacedPoints.Add(points[0]);
    Vector2 previousPoint = points[0];
    float dstSinceLastEvenPoint = 0;

    for (int segmentIndex = 0; segmentIndex < NumSegments; segmentIndex++)
    {
      Vector2[] p = GetPointsInSegment(segmentIndex);
      float controlNetLength = Vector2.Distance(p[0], p[1]) + Vector2.Distance(p[1], p[2]) + Vector2.Distance(p[2], p[3]);
      float estimatedCurveLength = Vector2.Distance(p[0], p[3]) + controlNetLength / 2f;
      int divisions = Mathf.CeilToInt(estimatedCurveLength * resolution * 10);
      float t = 0;
      while (t <= 1)
      {
        t += 1f/divisions;
        Vector2 pointOnCurve = Bezier.EvaluateCubic(p[0], p[1], p[2], p[3], t);
        dstSinceLastEvenPoint += Vector2.Distance(previousPoint, pointOnCurve);

        while (dstSinceLastEvenPoint >= spacing)
        {
          float overshootDst = dstSinceLastEvenPoint - spacing;
          Vector2 newEvenlySpacedPoint = pointOnCurve + (previousPoint - pointOnCurve).normalized * overshootDst;
          evenlySpacedPoints.Add(newEvenlySpacedPoint);
          dstSinceLastEvenPoint = overshootDst;
          previousPoint = newEvenlySpacedPoint;
        }

        previousPoint = pointOnCurve;
      }
    }

    return evenlySpacedPoints.ToArray();
  }

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

Вследствие чего я продолжил поиски и наткнулся на отличную статью «Moving Along a Curve with Specified Speed», на основе которой строятся дальнейшие рассуждения.

Значение скорости можно вычислить как $\sigma(t) = |V(t)|=|\frac{dX}{dt}|$, где $X(t)$ — функция сплайна. Чтобы решить поставленную задачу, нужно найти функцию $Y(t)=X(s)$, где $s$ — длина дуги от начальной точки до искомой, и это выражение устанавливает соотношение между $t$ и $s$.

Используя замену переменной дифференцирования, мы можем записать

$\frac{dY}{dt}=\frac{dX}{ds}\frac{ds}{dt}.$

Учитывая, что скорость — величина неотрицательная, получаем

$|\frac{dY}{dt}|=|\frac{dX}{ds}||\frac{ds}{dt}|=\frac{ds}{dt},$

так как $|\frac{dX}{ds}| = 1$ по условию неизменности модуля вектора скорости (вообще говоря, $|\frac{dX}{ds}|$ равен не единице, а константе, но для простоты выкладок примем эту константу равной единице).

Теперь получим соотношение между $t$ и $s$ в явном виде:

$s=g(t)=\int_{0}^{t}|\frac{dY(\tau)}{dt}|d\tau,$

откуда полная длина кривой $L$, очевидно, вычисляется как $g(1)$. С помощью полученной формулы можно, имея значение аргумента $t$, вычислить длину дуги до точки, в которую это значение $t$ отображается.

Нас же интересует обратная операция, то есть вычисление значения аргумента $t$ по заданной длине дуги $s$:

$t=g^{-1}(s).$

Так как не существует общего аналитического способа нахождения обратной функции, придется искать решение численное. Обозначим $F(t)=g(t)-s.$ При заданном $s$, требуется найти такое значение $t$, при котором $F(t)=0$. Таким образом, задача превратилась в задачу поиска корня уравнения, с чем может прекрасно справиться метод Ньютона.

Метод образует последовательность приближений вида

$t_{i+1}=t_i-\frac{F(t_i)}{F'(t_i)},$

где

$F'(t)=\frac{dF}{dt}=\frac{dg}{dt}=|\frac{dY}{dt}|.$

Для вычисления $F(t)$ требуется вычислить $g(t)$, вычисление которого, в свою очередь, требует численного интегрирования.

Выбор $t_0=\frac{s}{L}$ в качестве начального приближения теперь выглядит разумным (вспоминаем первый подход к решению задачи).

Далее выполняем итерации методом Ньютона до тех пор, пока погрешность решения не станет приемлемой, либо количество проделанных итераций — непозволительно большим.

При использовании стандартного метода Ньютона потенциально может возникнуть проблема. Функция $F(t)$ — неубывающая, так как её производная $F'(t)=|dY/dt|$ неотрицательна. Если вторая производная $F''(t)>0$, то функцию называют выпуклой и метод Ньютона на ней гарантированно сходится к корню. Однако, в нашем случае, $F''(t)$ может оказаться отрицательной, из-за чего метод Ньютона может сойтись за пределами диапазона определения $t\in[0,1]$. Автор статьи предлагает использовать модификацию метода Ньютона, которая, в случае если очередной результат итерации методом Ньютона не попадает в текущий ограничивающий корень интервал, вместо него выбирает середину интервала (метод дихотомии). Вне зависимости от результата вычисления на текущей итерации, диапазон сужается, а значит, метод сходится к корню.

Остается выбрать метод численного интегрирования — я воспользовался квадратурами метода Гаусса, так как он обеспечивает достаточно точный результат при небольших затратах.

Код функции, вычисляющей t(s)
public float ArcLength(float t) => Integrate(x => TangentMagnitude(x), 0, t);

private float Parameter(float length)
{
  float t = 0 + length / ArcLength(1);
  float lowerBound = 0f; 
  float upperBound = 1f;

  for (int i = 0; i < 100; ++i)
  {
    float f = ArcLength(t) - length;

    if (Mathf.Abs(f) < 0.01f)
      break;

    float derivative = TangentMagnitude(t);
    float candidateT = t - f / derivative;

    if (f > 0)
    {
      upperBound = t;
      if (candidateT <= 0)
        t = (upperBound + lowerBound) / 2;
      else
        t = candidateT;
    }
    else
    {
      lowerBound = t;
      if (candidateT >= 1)
        t = (upperBound + lowerBound) / 2;
      else
        t = candidateT;
    }
  }
  return t;
}


Код функции численного интегрирования
private static readonly (float, float)[] CubicQuadrature =
{(-0.7745966F, 0.5555556F), (0, 0.8888889F), (0.7745966F, 0.5555556F)};

public static float Integrate(Func<float, float> f, in float lowerBound, in float uppedBound)
{
  var sum = 0f;
  foreach (var (arg, weight) in CubicQuadrature)
  {
    var t = Mathf.Lerp(lowerBound, uppedBound, Mathf.InverseLerp(-1, 1, arg));
    sum += weight * f(t);
  }

  return sum * (upperBound - lowerBound) / 2;
}

Далее можно настроить погрешность метода Ньютона, выбрать более точный метод численного интегрирования, но задача, по сути, решена — был построен алгоритм вычисления аргумента $t$ сплайна для заданной длины дуги. Для объединения нескольких участков кривых в один и написания обобщенной процедуры вычисления $s(t)$ можно хранить длины всех участков и предварительно вычислять индекс участка, на котором требуется произвести итеративную процедуру модифицированным методом Ньютона.


Равномерно распределенные вдоль пути точки


Теперь самолет движется равномерно

Таким образом, нами были рассмотрены несколько способов параметризации сплайна относительно пройденного расстояния, в использованной мной в качестве источника статье в качестве варианта автор предлагает численно решать дифференциальное уравнение, но, по его собственной ремарке, предпочитает модифицированный метод Ньютона:
I prefer the Newton’s-method approach because generally t can be computed faster than with differential equation solvers.

В качестве заключения приведу таблицу времени расчета положения точки на представленной на скриншотах кривой в одном потоке на процессоре i5-9600K:
Кол-во вычислений Среднее время, мс
100 0,62
1000 6,24
10000 69,01
100000 672,81
Считаю, что такая скорость вычислений вполне позволяет применять их в различных играх и симуляциях. Буду рад уточнениям и критике по существу в комментариях.
Ads
AdBlock has stolen the banner, but banners are not teeth — they will be back

More

Comments 18

    +3

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

      +1
      Спасибо за статью.

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

        +1

        Что-то меня интуитивно тянет к дифференциальной геометрии и формулам Френе. Кмк красивое решение должно быть там. Как результат естественной параметризации.

          0
          ответ на одни вопросы порождает другие:
          крайняя правая точка кривой с учётом её тангенса создаёт нереальную для полёта самолёта кривую. Вот о таких моментах и их автоматическом исправлении путём добавления новых точке-автокоррекции тангенса, других уловках было бы интересно почитать а если ещё с учётом, чтоб это решалось в рамках навмеша было бы вообще замечательно(применять поиск пути и трейсить каждый кадр — решение очевидное и общеизвестное но явно не удобное и чересчур дорогое, интересно получить решение именно путём правки траектории и тангенса)
            0
            решение кажется неоптимальным, даже гугл дает быстрым поиском более быстрое решение: wiki.unity3d.com/index.php/Hermite_Spline_Controller?_ga=2.106834496.230616163.1595920433-1836999005.1595920433
              0
              Не обнаружил по приведённой ссылке процедуры, решающей ту же задачу. Единственный похожий метод public Vector3 GetHermiteAtTime(float timeParam) выполняет простую линейную интерполяцию. Не там смотрю?
              +1

              Надеюсь эту статью прочтут в Anywayanyday.


              А то их скринсэйвер жутко бесит угловатой траекторией.

              image

                0
                А нельзя реализовать управляемое моделирование движения самолета, задав ему точки, через которые он должен пролететь? Тогда траекторию будет задать чуть сложнее, зато двигаться он будет реалистичнее.
                  0
                  Можно, конечно, но тогда хорошо бы ему еще задать интерполяцию вращения по всем осям, и задавать точки маршрута с учетом совокупного воздействия подъемной силы и гравитации.
                  0
                  Алгоритм выглядит весьма неоптимальным — много раз (столько, сколько итераций в методе Ньютона) считается длина дуги от начальной точки, вместо того, чтобы считать длину дуги на интервале. И похоже это число итераций метода Ньютона весьма велико — возможно всегда равно 100 и условие
                  if (Mathf.Abs(f) < 0.01f)
                        break;
                  

                  никогда (или очень редко срабатывает). Иначе не знаю как объяснить такую медленную скорость работы — почти 1 секунда для 100000 точек.
                    0
                    Количество итераций методом Ньютона в среднем составляет от 1 до 3. И при численном интегрировании не имеет разницы в количестве вычислений, считается ли длина дуги от начала или на участке. Разве что точность будет чуть выше.
                      0
                      Имелось ввиду, что алгоритм можно переписать так, чтобы не считать длину дуги от начальной точки на каждой итерации, а использовать значение длины, вычисленное на предыдущей итерации, и считать длину участка дуги между предыдущей точкой и текущей.
                    +1

                    Я тоже решал такую задачу. Сделал так: разбил кривую на несколько частей, посчитал длину каждой части численно и построил интерполяционный полином, который выдавал t(L), где L — необходимая длина от начала сплайна, t аргумент функции сплайна. У меня кривая задавалась кубическим b-сплайном, полином, который выдавал t(L) был четвертой степени, но так как t(0)=0, то хватило четырех коэффициентов. Таким образом сплайн и полином L(t) имели по 4 коэффициента, их удобно было в шейдеры передавать. Плюс моего метода ещё и в том, что коэффициенты t(L) можно посчитать заранее, а не в рантайме.

                      0

                      Вроде бы спираль Эйлера (спираль Корню, клотоида) позволяет вычислять точки на себе в зависимости от пройденного расстояния вдоль кривой, а не на основе абстрактного t, как у Безье.


                      Когда-то пытался сделать сплайны на основе спиралей Эйлера как раз для простоты равномерного движения — пускать паровозики с вагончиками, но сходу что-то не пошло и сделал просто на дугах окружностей.

                        0
                        Интересная тема, не сталкивался раньше. Но, как мне кажется, это будет потеря вариативности в угоду меньшей вычислительной сложности.
                        +1
                        Нас же интересует обратная операция, то есть вычисление значения аргумента по заданной длине дуги
                        Таким образом, задача превратилась в задачу поиска корня уравнения, с чем может прекрасно справиться метод Ньютона.

                        Такую задачу ещё можно решать бинарным поиском.
                        Выходит, конечно же, медленнее, но для интервала длины 1 число итераций в худшем (любом) случае O(log(1/precision)) и нет проблем со схожденим, как в методе Ньютона. Для дебага — самое то.

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

                          Спасибо, но обозначения немного путанные мне кажется. Что такое X и Y не понял, и s и t. Под интегралом тоже путанно, там наверное должно быть dY/d\tau. И упомянуть что у нас переменные это вектора в 3d тоже не помешает.

                          Only users with full accounts can post comments. Log in, please.