Вступление

Мне кажется, что нам надо использовать меньше тригонометрии в компьютерной графике. Хорошее понимание проекций, отражений и векторных операций (как в истинном значении скалярного (dot) и векторного (cross) произведений векторов) обычно приходит с растущим чувством беспокойства при использовании тригонометрии. Точнее, я считаю, что тригонометрия хороша для ввода данных в алгоритм (для понятия углов это интуитивно понятный способ измерения ориентации), я чувствую, что что-то не так, когда вижу тригонометрию, находящуюся в глубинах какого-нибудь алгоритма 3D-рендеринга. На самом деле, я думаю, что где-то умирает котенок, когда туда закрадывается тригонометрия. И я не так беспокоюсь о скорости или точности, но с концептуальной элегантностью я считаю… Сейчас объясню.

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

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

Неправильный вариант поворота пространства или объекта

Пусть у нас будет функция, считающая матрицу поворота вектора вокруг нормализованного вектора $$ на угол $$. В любом 3D-движке или математической библиотеке реального времени будет одна такая функция, которая скорее всего слепо скопирована с другого движка, википедии или туториала по OpenGL… (да, к этому моменту вы должны признать, и в зависимости от вашего настроения, возможно переживать из-за этого).

Функция будет выглядеть примерно так:

mat3x3 rotationAxisAngle( const vec3 & v, float a )
{
    const float si = sinf( a );
    const float co = cosf( a );
    const float ic = 1.0f - co;

    return mat3x3( v.x*v.x*ic + co,     v.y*v.x*ic - si*v.z,  v.z*v.x*ic + si*v.y,
                   v.x*v.y*ic + si*v.z, v.y*v.y*ic + co,      v.z*v.y*ic - si*v.x,
                   v.x*v.z*ic - si*v.y, v.y*v.z*ic + si*v.x,  v.z*v.z*ic + co );
}

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

$$

Также, вы знаете, что иногда acosf() может вернуть странные значения, если скалярное произведение выходит за диапазон [-1; 1], и вы решаете изменить его значение так, чтобы оно попало в этот диапазон (прим. пер. to clamp) (в этот момент вы можете даже осмелиться винить точность вашего компьютера, потому что длина нормализованного вектора не равна в точности 1). К этому моменту погиб один котенок. Но пока вы не знаете об этом, вы продолжаете писать свой код. Далее вы вычисляете ось поворота, и вы знаете, что это — векторное произведение вектора $$ вашего объекта и выбранного направления $$, все точки в вашем объекте будут вращаться в плоскостях, параллельных той, которая определяется этими двумя векторами, так, на всякий случай… (котенка возродили и снова убили). В итоге, код выглядит как-то так:

const vec3   axi = normalize( cross( z, d ) );
const float  ang = acosf( clamp( dot( z, d ), -1.0f, 1.0f ) );
const mat3x3 rot = rotationAxisAngle( axi, ang );

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

const vec3   axi = normalize( cross( z, d ) );
const float  ang = acosf( clamp( dot( z, d ), -1.0f, 1.0f ) );
const float  co = cosf( ang );
const float  si = sinf( ang );
const float  ic = 1.0f - co;
const mat3x3 rot = mat3x3(
    axi.x*axi.x*ic + co,       axi.y*axi.x*ic - si*axi.z, axi.z*axi.x*ic + si*axi.y,
    axi.x*axi.y*ic + si*axi.z, axi.y*axi.y*ic + co,       axi.z*axi.y*ic - si*axi.x,
    axi.x*axi.z*ic - si*axi.y, axi.y*axi.z*ic + si*axi.x, axi.z*axi.z*ic + co);

Как вы уже могли заметить, мы выполняем довольно неточный и дорогой вызов acos, чтобы его сразу же отменить вычислением косинуса возвращаемого значения. И появляется первый вопрос: «а почему бы не пропустить цепь вызовов acos()--->cos() и сохранить процессорное время?» Более того, не говорит ли нам это о том, что мы делаем что-то неправильное и сильно запутанное, и что к нам приходит какой-то простой математический принцип, проявляющийся через упрощение этого выражения?

Вы можете поспорить, что упрощение нельзя сделать, так как вам нужен будет угол для вычисления синуса. Однако, это не так. Если вы знакомы с векторным произведением векторов, то вы знаете, что так же как скалярное произведение содержит косинус, векторное — содержит синус. Большинство графических программистов понимают зачем нужно скалярное произведение векторов, но не все понимают, зачем нужно векторное произведение (и используют его только для того, чтобы считать нормали и оси поворота). В основном, математический принцип, помогающий нам избавиться от пары cos/acos, также говорит нам о том, что там, где есть скалярное произведение, там возможно векторное произведение, сообщающее отсутствующую часть информации (перпендикулярную часть, синус).

$$

Правильный вариант поворота пространства или объекта

Теперь мы можем извлечь синус угла между $$ и $$ просто посмотрев на длину их векторного произведения… — помним, что $$ и $$ нормализованы! А это значит, что мы можем (мы должны!!) переписать функцию таким образом:

const vec3   axi = cross( z, d );
const float  si  = length( axi );
const float  co  = dot( z, d );
const mat3x3 rot = rotationAxisCosSin( axi/si, co, si );

и удостовериться, что наша новая функция построения матрицы поворота, rotationAxisCosSin(), нигде не вычисляет синусы и косинусы, а принимает их в качестве аргументов:

mat3x3 rotationAxisCosSin( const vec3 & v, const float co, const float si )
{
    const float ic = 1.0f - co;

    return mat3x3( v.x*v.x*ic + co,       v.y*v.x*ic - si*v.z,    v.z*v.x*ic + si*v.y,
                   v.x*v.y*ic + si*v.z,   v.y*v.y*ic + co,        v.z*v.y*ic - si*v.x,
                   v.x*v.z*ic - si*v.y,   v.y*v.z*ic + si*v.x,    v.z*v.z*ic + co );
}

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

mat3x3 rotationAlign( const vec3 & d, const vec3 & z )
{
    const vec3  v = cross( z, d );
    const float c = dot( z, d );
    const float k = (1.0f-c)/(1.0f-c*c);

    return mat3x3( v.x*v.x*k + c,     v.y*v.x*k - v.z,    v.z*v.x*k + v.y,
                   v.x*v.y*k + v.z,   v.y*v.y*k + c,      v.z*v.y*k - v.x,
                   v.x*v.z*K - v.y,   v.y*v.z*k + v.x,    v.z*v.z*k + c    );
}

Позже, Золтан Врана (Zoltan Vrana) подметил, что k можно упростить до k = 1/(1+c), что не только математически элегантнее выглядит, но также перемещает две особенности в k и, таким образом, вся функция ($$ и $$ параллельны) переходит в одну (когда $$ и $$ совпадают в этом случае нет четкого поворота). Финальный код выглядит как-то так:

mat3x3 rotationAlign( const vec3 & d, const vec3 & z )
{
    const vec3  v = cross( z, d );
    const float c = dot( z, d );
    const float k = 1.0f/(1.0f+c);

    return mat3x3( v.x*v.x*k + c,     v.y*v.x*k - v.z,    v.z*v.x*k + v.y,
                   v.x*v.y*k + v.z,   v.y*v.y*k + c,      v.z*v.y*k - v.x,
                   v.x*v.z*K - v.y,   v.y*v.z*k + v.x,    v.z*v.z*k + c    );
}

Мы не только избавились от трех тригонометрических функций и избавились от уродливого clamp'а (и нормализации!), но и концептуально упростили нашу 3D математику. Никаких трансцендентных функций, здесь используются только векторы. Векторы создают матрицы, изменяющие другие векторы. И это важно, потому что чем меньше тригонометрии в вашем 3D-движке, тем не только быстрее и яснее он становится, но, прежде всего, математически более элегантнее (правильнее!).