РАЗНОЕ
7.7. Кватернионы

Кватернион, он же гиперкомплексное число, представляет собой набор четырех чисел. Иногда будет удобно представлять себе кватернион как 4D-вектор, иногда как набор четырех чисел, иногда как число и 3D-вектор, а иногда и как гиперкомплексное число с тремя мнимыми единицами i, j, k; таким образом, имеем следующие представления:

q = [x1,x2,x3,x4] = [scalar,(vector)] = [x1,(x2,x3,x4)] = x1+x2*i+x3*j+x4*k.

Сложить или вычесть два кватерниона, а также умножить кватернион на число можно, как обычно, покомпонентно; с умножением ситуация более сложная. Умножение кватернионов должно в результате дать тоже кватернион, то есть конструкцию, содержащую лишь слагаемые вида r и r*l, где r - действительное число, а l - одна из мнимых единиц. Поэтому надо как-то определить операцию умножения для любых двух мнимых единиц. Определяется она так, что умножение получается некоммутативным, т.е. от перестановки мест множителей произведение меняется, и x*y != y*x. Поэтому умножение двух кватернионов приходится выполнять не по привычным правилам арифметики, а по следующим аксиомам:

a*(b*c) = (a*b)*c,      (ассоциативность)
(a+b)*c = a*c+b*c,      (транзитивность)
a*(b+c) = a*b+a*c,      (транзитивность)
a*1 = 1*a = a,          (существование единицы)
a*0 = 0*a = 0,          (существование нуля)
i*i = j*j = k*k = -1,   (свойство мнимых единиц)
i*j = -j*i = k.         (связь между мнимыми единицами i, j, k)

Из этих правил, кстати, следует, что

j*k = -k*j = i,
k*i = -i*k = j,

и получается такая вот таблица умножения комплексных единиц (умножение действительных чисел между собой и на комплексные единицы действует по обычным правилам, так что все свойства кватернионов определяются, в общем, этой таблицей):

второй множитель первый множитель i j k i -1 k -j j -k -1 i k j -i -1

Кроме того, из этих правил можно вывести правило для умножения кватернионов, заданных в форме [scalar,vector]:

q1 = [s1,v1],
q2 = [s2,v2],
q1*q2 = [s1*s2 - v1*v2, s1*v2 + s2*v1 + v1xv2].

Здесь v1*v2 - скалярное произведение векторов v1, v2; v1xv2 - векторное, все остальные произведения обычные (либо число на число, либо число на вектор).

Нужны же кватернионы для представления и интерполяции поворотов. Поворот относительно оси (x,y,z) (иными словами, поворот вокруг вектора (x,y,z), проведенного из начала координат) на угол angle представляется кватернионом q, лежащим на единичной 4D-сфере (то есть, 4D-вектором длины 1):

s = cos(angle/2),
v = (x,y,z) * sin(angle/2) / |(x,y,z)|,
q = [s,v].

Что интересно, в такой форме поворот, соответствующий комбинации поворотов q1 и q2, просто равен их произведению. В случае с 3D Studio это позволяет быстро и просто перевести сохраненные в CHUNK_TRACKROTATE относительные повороты в абсолютные: просто читаем эти самые повороты (а записаны они как раз в форме [angle,(x,y,z)], причем длина вектора (x,y,z) уже приведена к единичной), переводим их в кватернионную форму, получаем набор кватернионов q0, q1, ..., q(n-1), qn. Здесь q0 и так задает абсолютный поворот, а вот все остальные придется переводить (умножение здесь, конечно, кватернионное):

absolute_q0 = q0,
absolute_q1 = q1*absolute_q0,
absolute_q2 = q2*absolute_q1,
...
absolute_qn = qn*absolute_q(n-1).

Получаем набор кватернионов, задающих абсолютные повороты, или абсолютную ориентацию объекта в какие-то моменты времени. Для того же, чтобы получить поворот-ориентацию в любой момент времени, придется как-то интерполировать повороты между этими заданными ключевыми значениями.
Все кватернионы, задающие повороты, должны лежать на единичной 4D-сфере, поэтому простейший метод (линейная интерполяция) несколько усложнится: мы вынуждены интерполировать не по прямой между двумя векторами, а по дуге на этой 4D-сфере, являющейся сечением сферы плоскостью, проходящей через центр сферы и наши два вектора, то есть две точки на сфере. Все это называется сферической линейной интерполяцией (spherical linear interpolation, если скоращенно, slerp) и определяется следующим образом:

slerp(q1,q2,t) = (q1*sin((1-t)*a) + q2*sin(t*a)) / sin(a),

где t - локальное время (см.п.7.6), a - угол между векторами q1, q2;

0 <= t <= 1,
cos(a) = (q1,q2)/(|q1|*|q2|) = (q1,q2).

То есть q1, q2 здесь уже рассматриваем как 4D-вектора. Приведенную формулу нетрудно вывести (для лучшего понимания): нам нужна такая точка q, которая лежит на единичной сфере, лежит в одной плоскости с q1 и q2 и центром (то есть нулем), причем угол между векторами q и q1 меняется линейно и, таким образом, равен t*a. Раз точка лежит в одной плоскости с 0, q1, q2, то вектор q равен линейной комбинации векторов q1, q2:

q = k1*q1 + k2*q2,

где k1, k2 - какие-то (пока неизвестные) коэффициенты. q лежит на сфере, значит, длина q равна 1, отсюда имеем:

|q| = (q,q) = 1,
(k1*q1+k2*q2, k1*q1+k2*q2) = 1,
k1*k1*(q1,q1) + k2*k2*(q2,q2) + 2*k1*k2*(q1,q2) = 1,
k1*k1 + k2*k2 + 2*k1*k2*(q1,q2) = 1.

Угол между q и q1 равен t*a, отсюда:

cos(q,q1) = cos(t*a),
(q,q1) = cos(t*a),
k1*(q1,q1) + k2*(q1,q2) = cos(t*a),
k1 + k2*(q1,q2) = cos(t*a).

Получили систему уравнений для k1, k2:

k1 + k2*(q1,q2) = cos(t*a),
k1*k1 + k2*k2 + 2*k1*k2*(q1,q2) = 1,

или

k1 + k2*cos(a) = cos(t*a),
k1*k1 + k2*k2 + 2*k1*k2*cos(a) = 1.

Отсюда k1 = (cos(t*a) - k2*cos(a)), и получаем квадратное уравнение:

cos(t*a)^2 - 2*k2*cos(a)*cos(t*a) + k2^2*cos(a)^2 + k2^2 +
2*k2*cos(a)*cos(t*a) - 2*k2^2*cos(a)^2 = 1,

cos(t*a)^2 + k2^2*(1 - cos(a)^2) = 1,

k2^2 * sin(a)^2 = sin(t*a)^2,

k2 = sin(t*a) / sin(a),

k1 = cos(t*a) - sin(t*a)*cos(a) / sin(a) =
   = (cos(t*a)*sin(a) - sin(t*a)*cos(a)) / sin(a) =
   = sin(a - t*a) / sin(a) = sin((1 - t)*a) / sin(a).

Если a - очень маленький угол, настолько, что могут возникнуть ошибки при делении на sin(a), можно использовать обычную линейную интерполяцию (так как при маленьких значениях a sin(a) ~= a, sin(t*a) ~= t*a, и так далее).

Итак, мы умеем задавать повороты кватернионами, мы умеем их интерполировать линейно по множеству их возможных значений, то есть, поверхности сферы. Но хочется ведь интерполировать сплайнами Кочанека-Бартельса (далее везде, где используется термин "сплайны", подразумеваются именно такие сплайны), так как ориентация объекта должна меняться плавно, а не рывками, и желательно по совршенно той же траектории, что и в 3D Studio. Причем строить сплайны надо на поверхности четырехмерной сферы, иначе результаты интерполяции не будут соответствовать поворотам; кватернион-поворот должен обязательно лежать на единичной 4D-сфере. Естественное, возникает вопрос - как все это сделать?

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

g(t) = v1*(1-t)^3 + c1*3*t*(1-t)^2 + c2*3*t^2*(1-t) + v2*t^3,

и считаем ее значение в произвольно взятой точке t, используя только линейную интерполяцию (linear interpolation, lerp): пусть

lerp(a, b, t) = a*(1-t) + b*t,

тогда g(t) можно посчитать вот так:

tmp1 = lerp(v1, c1, t),
tmp2 = lerp(c1, c2, t),
tmp3 = lerp(c2, v2, t),
tmp4 = lerp(tmp1, tmp2, t),
tmp5 = lerp(tmp2, tmp3, t),
g(t) = lerp(tmp4, tmp5, t).

Нам же надо интерполяцию сплайнами по поверхности сферы, это можно получить, всего-навсего заменив в приведенных выше формулах линейную интерполяцию lerp на наш сферический вариант, slerp. Далее, сравнивая g(t) с полученной в п.7.6. интерполяционной функцией f(t) можно заметить, что, если

p1 = v1              <=>   v1 = p1,
r1 = (c1 - v1) * 3   <=>   c1 = (p1 + r1) / 3,
r2 = (v2 - c2) * 3   <=>   c2 = (p2 - r2) / 3,
p2 = v2              <=>   v2 = p2,

то g(t) совпадает с f(t). Длинный и нудный вывод для v1, v2, c1, c2 через функции lerp()/slerp() делать, пожалуй, смысле нет, так что ограничимся конечными результатами. А именно, для каждой точки-ключа cur имеем

g1 = slerp(cur, prev, -(1+bias)/3.0);
g2 = slerp(cur, next, (1-bias)/3.0);
g3 = slerp(g1, g2, 0.5 + 0.5*continuity);
g4 = slerp(g1, g2, 0.5 - 0.5*continuity);
cur.ra = slerp(cur, g3, (tension-1));
cur.rb = slerp(cur, g3, -(tension-1));

Для начальной и конечной точки, соответственно, имеем следующее:

q0.rb = slerp(p0, p1, (1-tension)*(1+continuity*bias)/3.0);
qn.ra = slerp(pn, p(n-1), (1-tension)*(1-continuity*bias)/3.0);

При интерполяции между какими-то точками a и b просто полагаем

v1 = a,
c1 = a.rb,
c2 = b.ra,
v2 = b,

и считаем g(t) по приведенному выше алгоритму. Здесь мы до сих пор не учли параметры ease to и ease from, но это дело одной строки кода - посчитать на самом деле надо не g(t), а g(ease(t)). Впрочем, обычно ease(t) = t. Что это за функция ease(), откуда она берется, для чего нужна и как рассчитывается, написано в п.7.6.

Таким образом, получаем кватернион, соответствующий повороту, задающий ориентацию объекта. Осталось выяснить, как из этого кватерниона получить что-нибудь более привычное - скажем, матрицу поворота. С одной стороны, вспомнив, как делается перевод в кватернионную форму для поворота на угол angle относительно оси (x,y,z), можно написать, что

q = [s,v],
angle = 2 * arccos(angle),
(x,y,z) = v / sin(angle/2),

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

q = [w,(x,y,z)],

    [ 1-2*(y*y+z*z)   2*(x*y-w*z)     2*(x*z+w*y)   ]
A = [ 2*(x*y+w*z)     1-2*(x*x-z*z)   2*(y*z-w*x)   ]
    [ 2*(x*z-w*y)     2*(y*z+w*x)     1-2*(x*x-y*y) ].

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

• все повороты надо заранее перевести в кватернионную форму;

• расчеты производных и интерполяция делаются по формулам, данным именно в этом пункте, а не по "обычным" для сплайнов;

• из кватернионной формы результирующий поворот обычно надо перевести в привычную матричную.

Все остальное совпадает с планом действий при "обычной" сплайновой интерполяциии и изложено в пункте 7.6.