Избегаем тригонометрии

Автор оригинала: Inigo Quilez
  • Перевод

Вступление


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

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

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

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


Пусть у нас будет функция, считающая матрицу поворота вектора вокруг нормализованного вектора $\vec{v}$ на угол $a$. В любом 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 );
}

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

$\vec{a} \cdot \vec{b} = a_xb_x + a_yb_y = ab \cos{\angle \vec{a} \vec{b} }$

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

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, также говорит нам о том, что там, где есть скалярное произведение, там возможно векторное произведение, сообщающее отсутствующую часть информации (перпендикулярную часть, синус).

$||\vec{a}\times\vec{b}|| = ab\sin{\angle\vec{a}\vec{b}}$

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


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

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 и, таким образом, вся функция ($\vec{d}$ и $\vec{z}$ параллельны) переходит в одну (когда $\vec{d}$ и $\vec{z}$ совпадают в этом случае нет четкого поворота). Финальный код выглядит как-то так:

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-движке, тем не только быстрее и яснее он становится, но, прежде всего, математически более элегантнее (правильнее!).

Только зарегистрированные пользователи могут участвовать в опросе. Войдите, пожалуйста.

Альтернативное голосование

Поддержать автора
Поделиться публикацией
AdBlock похитил этот баннер, но баннеры не зубы — отрастут

Подробнее
Реклама

Комментарии 17

    +8
    Для перехода из (x, y, z) к (x', y', z'), где z' = d
    нужна не матрица поворота, а матрица преобразования координат, которая составлена из ортов x, y, z, заданных в координатах (x', y', z'). Которую, в свою очередь, можно найти обратив матрицу, составленную из ортов x', y', z', заданных в x, y, z.

    Возможно, я сноб, но мне почему-то казалось, что это очевидно. Матрица поворота здесь не нужна.
      –2
      Это будут те же самые арифметические действия, только «в профиль». Вам вычисление компонентов векторного произведения AzBx — AxBz не напоминает вычисление детерминанта по минору?
        0

        Тут вопрос не в том какие это будут действия, а в том как к ним идти

          0
          Вы как учитель школьный — толкаете единственно верный ход мыслей?

          P.S. Есть разные инструменты, подходы и интерпретации. Собственно, разный тип мышления выбирает тот интсрумент, который ему понятен.
            +12
            Ну, если человек не прогуливал на первом курсе линейную алгебру, то он знает оба пути, и легко видит, что названный мной — более прямой.
        0
        Пытаясь расшифровать ваше послание, надо вот такую матрицу вычислить?
            | x.x’ y.x’ z.x’ | -1
        M = | x.y’ y.y’ z.y’ |
            | x.z’ y.z’ z.z’ |
        
        +11
        Которую, в свою очередь, можно найти обратив матрицу

        В большинстве случаев (когда нет масштабирования и оси ортогональны) достаточно транспонировать, что ещё дешевле.
          0
          Тем более. Спасибо за дополнение.
          –1
          del
            +5
            Насколько устойчива схема, когда вектора d и z почти противоположны? Ведь c будет стремиться к -1, а k, соответственно, к бесконечности.

            Кстати, схема расчёта поворота ещё сильнее упрощается в плане подготовки преобразований на CPU перед отдачей в GPU-конвейер и размещения в памяти, если вычислять не матрицу поворота, а кватернион:

            vec4 quaternionByTwoVectors(vec3 d, vec3 z) {
                const vec3 vec = cross(d, z);
                const float c = dot(d, z);
                const float k = sqrt((1 + c) * 2);
            
                return vec4(c.x / k, c.y / k, c.z / k, k / 2.0f);
            }


            ну и дальше в вертексном шейдере, при трансформации просто перемножить этот кватернион с координатой каждой вершины (что чуть дороже, чем матричное умножение)

            vec3 rotateVectorByQuaternion(vec4 q, vec3 v) {
              return v + 2.0 * cross(q.xyz, cross(q.xyz, v) + q.w * v);
            }


            P.S.: может, кто подизассемблировать итоговый код, сколько mul/add операций выходит в rotateVectorByQuaternion против перемножения mat3x3 на vec3?
              –2
              И в итоге вместо понятной rotationAxisAngle с понятными любому аргументами мы имеем rotationAxisCosSin с косинусами и синусами вместо аргументов. Да, практически любому ясно как это работает. Но читабельность кода превращается в г…
              Поэтому если кусок кода не критичной по скорости — читабельность в приоритете, оставляем ось и угол.
              Если кусок кода критичен по скорости — реализуем всю математику внутри функции не разбивая её на кучу других.
              А вообще: преждевременная оптимизация — зло. Именно поэтому ни в одном игровом движке мат либа не содержит таких извращений, а дает простые и понятные методы.
                0
                А почему если есть косинус не вычислить синус используя основное тригонометрическое тождество? Или корень более тяжелая операция чем векторное произведение с последующей нормализацией?
                  +1

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

                    0
                    4 операции умножения в быстром алгоритме вычисления обратного корня (1/sqrt(x)): en.wikipedia.org/wiki/Fast_inverse_square_root

                    Но пишут, что sse или avx инструкции (rsqrtss и vrsqrtss соответственно) работают быстрее (за 3-5 тактов)
                      0

                      Так то обратный корень, а для синуса нужен прямой.


                      Вот как раз при нормализации можно и обратный использовать, ежели точность устраивает (а она может и не устроить).

                    0
                    Как минимум, точность убивается в ноль. Допустим, у нас есть угол в районе 0, скажем, 10^-4. Его косинус находится где-то в районе 1, квадрат соответственно тоже. Машинная точность float в районе 1 где-то 10^-8. вычитаем из 1 квадрат косинуса, погрешность остается не меньше 10^-8. берем корень и наша погрешность внезапно становится 10^-4, что соразмеримо с синусом угла, который хотим измерить. В итоге функция измерения синуса через косинус имеет явное квантование в районе 0, что часто неприменимо. А если считать на 16битных числах, то там вообще погрешность 10^-2, что уже вообще использовать нельзя почти ни для чего
                      0
                      Дык, а знак-то к этому корню откуда выкопать?

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

                    Самое читаемое