РАЗНОЕ
7.6. Сплайны Кочанека-Бартельса

Пусть у нас имеется набор моментов времени и заданные значения какой-то функции (например, координат положения объекта) в эти моменты. Для того, чтобы получить значение функции в какой-то промежуточный момент времени, между заданными моментами ее придется интерполировать. Самый простой случай, конечно, линейная интерполяция (точнее говоря, кусочно-линейная), когда между любыми двумя заданными точками p1, p2 значение интерполируется по следующей формуле:

f(t) = p1.value + (p2.value - p1.value) * t.

Здесь t - "локальное" время, t=0 в точке p1 и t=1 в точке p2, pi.value - значение интерполируемой величины в точке pi. Перевод "локального" времени t в "глобальное" время T и наоборот тривиален:

t = (T - p1.T) / (p2.T - p1.T),
T = p1.T + t * (p2.T - p1.T).

Результат получится такого вида:

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

Интерполяция сплайнами - такой метод, который дает в результате гладкую функцию. Достигается это следующим образом: между каждыми двумя точками функция интерполируется кубической параболой (то есть, полиномом третьей степени, вида P(x) = a*x^3+b*x^2+c*x+d), причем функция подбирается так, чтобы на концах отрезка (как раз в точках задания значений) совпадали с какими-то заранее заданными числами и сами значения функции, и значения производной. Совпадение самих значений функции требуется для непрерывности получаемой приближающей функции, а совпадение значений производной - для непрерывности производной, то есть для "гладкости" приближающей функции. Таким образом, результат интерполяции у нас получается составленным не из прямых, как при кусочно-линейной интерполяции, а из "сшитых" кусков разных кубических парабол. Причем эта кривая получается еще и гладкой, так как мы "сшиваем" не только сами значения функции, но и значения ее производной.

Пусть на концах отрезка "локального" времени [0,1] заданы значения функции p1, p2 и значения производной r1, r2. Приближающая функция - кубическая, то есть

f(t) = a*t^3 + b*t^2 + c*t + d,
f'(t) = 3*a*t^2 + 2*b*t + c,

откуда и получаем:

f(0) = p1,
f(1) = p2,
f'(0) = r1,
f'(1) = r2,

d = p1,
a + b + c + d = p2,
c = r1,
3*a + 2*b + c = r2.

Решив эту систему уравнений, получаем:

a = 2*p1 - 2*p2 + r1 + r2,
b = -3*p1 + 3*p2 - 2*r1 -r2,
c = r1,
d = p1,

отсюда

f(t) = (2*p1 - 2*p2 + r1 + r2) * t^3 +
(-3*p1 + 3*p2 - 2*r1 - r2) * t^2 + r1 * t + p1,

или, что равносильно,

f(t) = p1 * (2*t^3 - 3*t^2 + 1) + r1 * (t^3 - 2*t^2 + t) +
p2 * (-2*t^3 + 3*t^2) + r2 * (t^3 - t^2).

Подставляем сюда значения интерполируемой величины в начале и конце отрезка (p1 и p2), желаемые значения производных в начале и конце отрезка (r1 и r2), локальное время для отрезка интерполяции t, и получаем интерполированное значение функции в любой точке отрезка - кусок кубической параболы, который проходит через точки (0,p1), (1,p2) и имеет в точках t=0 и t=1 производные r1 и r2 соответственно.

Осталось выяснить, откуда брать значения производных r1, r2. Они считаются через заданные (они же ключевые) значения нашей функции в самой точке-ключе, следующей и предыдущей точке, а также параметры сплайна tension, continuity и bias в данной точке. Все это (значения функции и набор параметров сплайна в каком-то наборе точек) задается извне; то есть, например, читается из 3DS-файла.

Введем следующие обозначения. Пусть cur - текущая точка, next - предыдущая, prev - следующая, r - производная в текушей точке cur, которую мы и должны как-то посчитать. Пусть

g1 = cur.value - prev.value,
g2 = next.value - cur.value,
g3 = g2 - g1.

Тогда по умолчанию (это когда параметры tension, continuity, bias равны 0) производная считается как

r = g1 + 0.5 * g3.

Параметры tension, continuity и bias соответственно изменяют вес g3, а также значения r, g1, g2; с учетом этих параметров формула для производной выглядит следующим образом:

g1 = (cur.value - prev.value) * (1 + bias),
g2 = (next.value - cur.value) * (1 - bias),
g3 = g2 - g1,
ra = (1 - tension) * (g1 + 0.5 * g3 * (1 + continuity)),
rb = (1 - tension) * (g1 + 0.5 * g3 * (1 - continuity)).

Здесь уже появляются два разных значения производной ra и rb. ra - это то значение производной, которое используется, когда точка является началом отрезка сплайн-интерполяции; rb - когда концом. То есть, при интерполяции между какими-то точками p1, p2 используются значения r1 = p1.ra, r2 = p2.rb. При continuity = 0 они (ra и rb) равны, при continuity = 1 ra удваивается, а rb становится равным нулю. Таким образом, параметр continuity контролирует "изломанность", негладкость сплайна.

Геометрический смысл всех этих параметров нагляднее всего можно показать в случае, когда value - 2D-вектор. Тогда cur, prev, next - какие-то точки на плоскости, сплайн - проведенная через них кривая, r - вектор-градиент кривой в точке cur, он же является касательной (точнее, направляющим вектором прямой, являющейся касательной) в этой точке. Параметр tension непосредственно изменяет длину градиента; bias меняет веса g1, g2 в g3; continuity заставляет конец градиента "ездить" вдоль g3. Все это соотвествующим образом влияет на вид самой кривой. Для лучшего понимания посмотрите на приведенную картинку, представьте себе гладкую линию (кривую), проведенную через точки prev, cur, next и учтите, что r - это желаемое нами положение касательной прямой к этой кривой в точке cur.

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

ra = next.value - cur.value,
rb = (1.5 * (next.value - cur.value) - 0.5 * next.ra) * (1 - tension).

Для конечной точки

ra = (1.5 * (cur.value - prev.value) - 0.5 * prev.rb) * (1 - tension),
rb = cur.value - prev.value.

ra для начальной точки и rb для конечной на самом деле не используются, они приведены здесь лишь для того, чтобы формулы давали правильный результат в случае, если есть всего две точки p0, p1; тогда

p0.rb = p1.ra = (p1.value - p0.value) * (1 - tension).

Осталось упомянуть про параметры ease to и ease from. Они контролируют нелинейность движения по полученной кривой в зависимости от времени - или, если угодно, нелинейность самого времени. Не особо вникая в детали, корректно отработать все это можно с помощью следующего куска кода:

// ...

float ease(float t, float a, float b) {
  float k;
  float s = a + b;

  if (s == 0.0) return t;
  if (s > 1.0) {
    a = a / s;
    b = b / s;
  }
  k = 1.0 / (2.0 - a - b);
  if (t < a)
    return ((k / a) * t * t);
  else {
    if (t < 1.0 - b){
      return (k * (2 * t - a));
    } else {
      t = 1.0 - t;
      return (1.0 - (k / b) * t * t);
    }
  }
}

// ...

time = ease(time, beginKey.easeFrom, endKey.easeTo);

// ...

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

• найти те два ключа, между которыми попадает этот момент времени (если ключ всего один, значение просто постоянно);
• перевести глобальное время в локальное для отрезка между этими ключами;
• скорректировать локальное время с помощью функции ease();
• посчитать "выходящую" производную rb для начального ключа и "входящую" производную ra для конечного ключа;
• по значениям величины в ключах, посчитанным производным и локальному времени посчитать интерполированное значение.