Вычисление касательного базиса налету / forpes.ru

09.06.2022 12:11

SergeySib 0 741 Источник

Этот пост написан автором в далёком 2013 году как дополнение к статье из 2006 года [4]. Тем не менее, предпосылки к отказу от хранения предрассчитанных касательных справедливы и в 2022, а вывод формул изложен понятно и без излишеств.

Мотивация к разработке этого метода заключается в создании касательного базиса налету во фрагментном шейдере, что по иронии диаметрально противоположно мотивации из статьи [2] от 1997 года:

В статье представлен метод бамп-маппинга с минимальным требованием к железу помимо того, которое нужно для освещения по Фонгу. Мы избавимся от затратного попиксельного восстановления касательного пространства...

Так как за окном более не 1997 год, вычисление касательного пространства во фрагментах потенциально даёт некоторый выигрыш:

Понижается сложность инструментария для создания ассетов.
Экономится шина и память для выборки и хранения вершин.
Экономятся интерполяторы атрибутов между геометрическими и фрагментной стадиями.
Отпадает необходимость преобразовывать касательные в вершинном шейдере
Становится проще накладывать карты нормалей при нелинейных деформациях

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

Обычно нормалмаппинг подаётся в статьях не совсем правильно, что я проиллюстрирую простой метафорой на C++. Пусть есть класс Vector3 для векторов, и есть класс Covector3 для ковекторов. Последний будет копией класса обычных векторов за исключением того что он по другому ведёт себя при трансформациях (за подробным введением в ковекторы и сопряжённые пространства смотрите эту статью). Как вы возможно знаете, нормаль - это на самом деле ковектор, а значит так мы её и объявим. Теперь вообразим такую функцию:

Vector3 tangent;
Vector3 bitangent;
Covector3 normal;
Covector3 perturb_normal( float a, float b, float c )
{
    return
        a * tangent +
        b * bitangent +
        c * normal;
        // ^^^^ compile-error: type mismatch for operator +
}

Эта функция замешивает векторы и ковекторы в одно выражение, что ведёт к ошибке компиляции. Если normal объявлена типом Covector3, то tangent и bitangent должны быть того же типа, иначе из них нельзя сделать единый базис. В реальном мире в коде шейдера всё это будет объявлено типом vec3 и будет работать, хоть и не всегда.

Математическая ошибка компиляции

Прискорбно, но "касательное пространство" для нормалмаппинга было введено авторами статьи [2] через вышеупомянутое несоответствие. Всё работает правильно до тех пор пока касательный базис ортогонален. Но когда дело доходит до восстановления касательного базиса во фрагментном шейдере, приходится иметь дело с неортогональной проекцией на экран. Именно по этой причине в исходной статье я ввёл $\vec{T}$ (который должен быть дуальным касательным, "кокасательным") и $\vec B$ (который должен называться ко-би-касательной, как бы упорото это ни звучало) как ковекторы, иначе алгоритм не заработал бы. Таким образом, TBN-базис должен называться не касательным, а кокасательным.

Экскурс в историю: возмущённые нормали Блинна

В этом разделе я покажу как определение $\vec T$ и $\vec B$ ковекторами естественным образом следует из оригинальной статьи Блинна [1] о бамп-маппинге. Блинн рассматривает криволинейную параметрическую поверхность (например, участок поверхности Безье), на которой он определяет векторы $\vec{p}_u$ и $\vec{p}_v$ как производные координаты точки поверхности $\vec p$ по и .

В этом контексте он пишет нижние индексы как краткое обозначение производных, то есть $\vec{p}_u = \frac{\partial\vec p}{\partial u}$ . Нормаль к поверхности он вводит $\vec N = \vec{p}_u\times \vec{p}_v$ и функцию высот . Наконец, он получает первое приближение формулы возмущённой нормали:

$\vec{N}' \simeq \vec{N} + \frac{f_u \vec{N} \times \vec{p}_v + f_v \vec{p}_u \times \vec{N}}{|\vec{N}|}$

Обратите внимание на слагаемые $\vec{N} \times \vec{p}_v$ и $\vec{p}_u \times \vec{N}$ . Они перпендикулярны $\vec{p}_u$ и $\vec{p}_v$ в касательной плоскости, и все вместе образуют базис для смещения и . Кроме того, они ковекторы, что видно из их поведения при трансформациях. А значит их сложение с нормалью не приводит к несоответствию типов. Если мы разделим эти слагаемые на $|\vec N|$ и обратим их знак, то придём к следующему определению $\vec T$ и $\vec B$ :

$\begin{align*} \vec{T} &= -\frac{\vec{N} \times \vec{p}_v}{|\vec{N}|^2} = \nabla u, & \vec{B} &= -\frac{\vec{p}_u \times \vec{N}}{|\vec{N}|^2} = \nabla v, \end{align*}$ $\vec{N}' \simeq \hat{N} - f_u \vec{T} - f_v \vec{B} ,$

где крышечка обозначает единичную нормаль. Это определение совпадает с таковым в [4]. $\vec T$ можно понимать как нормаль к плоскости заданного , а $\vec B$ — как нормаль к плоскости заданного . В итоге у нас есть три взаимно-перпендикулярных вектора (точнее, ковектора) $\vec T$ , $\vec B$ и $\vec N$ , составляющие базис сопряжённого касательного пространства. А ещё $\vec T$ и $\vec B$ это градиенты и . Модуль градиентов определяет величину изменения высоты, о чём поговорим ниже.

Забудьте то чему вас учили раньше

Многие авторы ошибочно берут $\vec T$ и $\vec B$ за $\vec{p}_u$ и $\vec{p}_v$ , что верно лишь пока векторы ортогональны. Давайте забудем "касательные", возьмём "кокасательные" и повторим историческое развитие с этой точки зрения: Пирси с соавторами [2] предрасчитывают значения и (изменение высоты на единицу изменения текстурных координат) и сохраняет их в текстуру. Они называют это "картой нормалей", но фактически это "карта крутизны уклонов", которая недавно была переизобретена под названием карта производных. В такую карту нельзя закодировать горизонтальные нормали, потому что тогда пришлось бы кодировать бесконечное значение крутизны. В дополнение к этой карте нужно ещё хранить некий множитель силы смещения. Килгард [3] вводит карты нормалей в современном понимании как закодированный оператор вращения, полностью избавляясь от аппроксимаций путём явного определения возмущённой нормали:

$\vec{N}' = a \vec{T} + b \vec{B} + c \hat{N} ,$

где , и вычитываются из текстуры. Многие думают, что в карте нормалей хранятся нормали, но это верно лишь отчасти. Идея Килгарда была следующей: так как невозмущённая нормаль - это [0, 0, 1], то достаточно хранить лишь последний столбец матрицы вращения чтобы вычислить возмущённую нормаль. В общем да, карта нормалей хранит векторы, соответствующие возмущённым нормалям, но фактически это закодированный оператор вращения. Трудности начинаются когда нужно сблендить несколько карт нормалей, потому что это интерполяция вращений со всеми вытекающими. Более подробный обзор здесь.

Вывод формул двойственного касательного базиса

Перед нами стоит задача, обратная к решённой в работе Блинна. Из карты нормалей мы знаем возмущённую нормаль, но не знаем двойственный касательный базис. Определим неизвестные (ко)касательные как градиенты текстурных координат и : $\vec{T} = \nabla u$ и $\vec{B} = \nabla v$ , а сами текстурные координаты неявно определим как функции координат вершины :

$\begin{align*} \mathrm{d} u &= \vec{T} \cdot \mathrm{d} \vec{p} , & \mathrm{d} v &= \vec{B} \cdot \mathrm{d} \vec{p} , \end{align*}$

где точка - это скалярное произведение. Градиенты постоянны для интерполированных значений внутри треугольника, потому введём разности $\Delta u_1$ , $\Delta u_2$ , $\Delta v_1$ , $\Delta v_2$ , а так же $\Delta \vec{p}_1$ и $\Delta \vec{p}_2$ . Искомые базисные векторы должны удовлетворять условиям:

$\begin{align*} \Delta u_1 &= \vec{T} \cdot \Delta \vec{p_1} , & \Delta v_1 &= \vec{B} \cdot \Delta \vec{p_1} , \\ \Delta u_2 &= \vec{T} \cdot \Delta \vec{p_2} , & \Delta v_2 &= \vec{B} \cdot \Delta \vec{p_2} , \\ 0 &= \vec{T} \cdot \Delta \vec{p_1} \times \Delta \vec{p_2} , & 0 &= \vec{B} \cdot \Delta \vec{p_1} \times \Delta \vec{p_2}. \end{align*}$

Первые две строки следуют из определений, а последняя из условия что (ко)касательные перпендикулярны нормали. Последняя строка нужна потому что без неё задача недоопределена. Отсюда можно выразить $\vec T$ в матричной форме:

$\vec{T} = \begin{pmatrix} \Delta \vec{p_1} \\ \Delta \vec{p_2} \\ \Delta \vec{p_1} \times \Delta \vec{p_2} \end{pmatrix}^{-1} \begin{pmatrix} \Delta u_1 \\ \Delta u_2 \\ 0 \end{pmatrix}.$

$\vec B$ и $\Delta v$ находятся аналогично.

Переходим к коду шейдера

Полученный результат выглядит устрашающе потому что нужно обращать матрицу в каждом фрагменте! Однако, используя некоторые симметрии, код обращения можно свести к минимуму. Ниже приведён пример функции на GLSL для обращения матрицы 3х3. Похожий код на HLSL появлялся в книге, но с тех пор я его сильно оптимизировал. Смотрите как из веткорных произведений можно получить определитель и присоединённую матрицу:

mat3 inverse3x3( mat3 M )
{
    mat3 M_t = transpose( M ); 
    float det = dot( cross( M_t[0], M_t[1] ), M_t[2] ); 
    mat3 adjugate = mat3( cross( M_t[1], M_t[2] ), cross( M_t[2], M_t[0] ), cross( M_t[0], M_t[1] ) ); 
    return adjugate / det;
}

Строки матрицы из предыдущей части можно подставить в код, раскрыть скобки, привести подобные, получив новое выражение для $\vec T$ . Определитель превращается в $\left| \Delta \vec{p_1} \times \Delta \vec{p_2} \right|^2$ , а присоединённую матрицу можно записать с помощью двух новых переменных $\Delta \vec{p_1}_\perp$ и $\Delta \vec{p_2}_\perp$ :

$\vec{T} = \frac{1}{\left| \Delta \vec{p_1} \times \Delta \vec{p_2} \right|^2} \begin{pmatrix} \Delta \vec{p_2}_\perp \\ \Delta \vec{p_1}_\perp \\ \Delta \vec{p_1} \times \Delta \vec{p_2} \end{pmatrix}^\mathrm{T} \begin{pmatrix} \Delta u_1 \\ \Delta u_2 \\ 0 \end{pmatrix} ,$ $\begin{align*} \Delta \vec{p_2}_\perp &= \Delta \vec{p_2} \times \left( \Delta \vec{p_1} \times \Delta \vec{p_2} \right) , \\ \Delta \vec{p_1}_\perp &= \left( \Delta \vec{p_1} \times \Delta \vec{p_2} \right) \times \Delta \vec{p_1} . \end{align*}$

Внимательный читатель возможно догадался, что $\Delta \vec{p_1}_\perp$ и $\Delta \vec{p_2}_\perp$ это перпендикуляры к сторонам треугольника в плоскости треугольника. И конечно же они ковекторы, образующие хороший базис сопряжённого касательного пространства. Для упрощения дальнейшего изложения заметим следующее:

Последняя строка матрицы неважна, потому что умножается на 0.
Остальные строки матрицы содержат перпендикуляры ( $\Delta \vec{p_1}_\perp$ и $\Delta \vec{p_2}_\perp$ ), которые после транспонирования умножаются на численные производные текстурных координат.
Вместо нормали к грани ( $\Delta \vec{p_1} \times \Delta \vec{p_2}$ ) можно взять интерполированную нормаль к вершине, что проще и даёт лучший визуальный результат.
Определитель (выражение $\left| \Delta \vec{p_1} \times \Delta \vec{p_2} \right|^2$ ) можно отбросить, о чём ниже в разделе про инвариантность при масштабировании.

В итоге код получается простым и надёжным:

mat3 cotangent_frame( vec3 N, vec3 p, vec2 uv )
{
    // Численные производные
    vec3 dp1 = dFdx( p );
    vec3 dp2 = dFdy( p );
    vec2 duv1 = dFdx( uv );
    vec2 duv2 = dFdy( uv );
    
    // Решаем линейную систему
    vec3 dp2perp = cross( dp2, N );
    vec3 dp1perp = cross( N, dp1 );
    vec3 T = dp2perp * duv1.x + dp1perp * duv2.x;
    vec3 B = dp2perp * duv1.y + dp1perp * duv2.y;
    
    // Собираем масштабно-инвариантный базис
    float invmax = inversesqrt( max( dot(T,T), dot(B,B) ) );
    return mat3( T * invmax, B * invmax, N );
}

Инвариантность при масштабировании

В выведенном выше выражении оставался определитель $\left| \Delta \vec{p_1} \times \Delta \vec{p_2} \right|^2$ . Из-за этого дуальные касательные $\vec T$ и $\vec B$ изменяются обратно пропорционально масштабу модели. Так получается потому что они градиенты. Если масштаб треугольника увеличивается, а всё остальное остаётся неизменным, то изменение текстурных координат на единицу изменения координаты вершин становится меньше. Из-за этого уменьшается $\vec{T} = \nabla u = \left( \frac{\partial u}{\partial x}, \frac{\partial u}{\partial y}, \frac{\partial u}{\partial z} \right)$ и (аналогично) $\vec B$ по отношению к $\vec N$ . В итоге, уменьшается возмущение нормали при увеличении масштаба треугольника, как будто карта высот растянулась.

Хотя такое поведение абсолютно логично и корректно, оно мешает накладывать нормалмапы на геометрию различного масштаба. Для решения этой проблемы (как видно из кода), я отбрасываю определитель и нормирую $\vec T$ и $\vec B$ по длине большего из них. Этот хак сохраняет относительные длины $\vec T$ и $\vec B$ , и даже скошенный или растянутый базис сопряжённого касательного пространства обрабатывается корректно.

Бесперспективная оптимизация

В качестве агрессивной оптимизации я попробовал положить $\Delta \vec{p_1} = \Delta \vec{p_2}_\perp$ и $\Delta \vec{p_2} = \Delta \vec{p_1}_\perp$ . Это значит что треугольник прямоугольный, и перпендикуляр к одному катету совпадает со вторым катетом. Во фрагментном шейдере это условие выполняется когда проекция на экран не вносит перспективных искажений. В [4] есть хороший рисунок с демонстрацией этого факта. Такая оптимизация экономит ещё два векторных произведения, но по-моему сильно вредит качеству если фактически перспективное искажение есть.

Собираем всё воедино

Для полноты картины я покажу как с помощью сопряжённого касательного пространства возмутить интерполированную вершинную нормаль. Этим занимается функция perturb_normal, используя вектор из фрагмента в камеру вместо координаты вершины. Это нормально, потому что при вычислении разности положение камеры сокращается.

vec3 perturb_normal( vec3 N, vec3 V, vec2 texcoord ) 
{ 
    // N - интерполированная нормаль,
    // V - интерполированный вектор из вершины в камеру
    vec3 map = texture2D( mapBump, texcoord ).xyz;
#ifdef WITH_NORMALMAP_UNSIGNED
    map = map * 255./127. - 128./127.; 
#endif 
#ifdef WITH_NORMALMAP_2CHANNEL 
    map.z = sqrt( 1. - dot( map.xy, map.xy ) ); 
#endif 
#ifdef WITH_NORMALMAP_GREEN_UP 
    map.y = -map.y; 
#endif
    mat3 TBN = cotangent_frame( N, -V, texcoord ); 
    return normalize( TBN * map ); 
}

varying vec3 g_vertexnormal; 
varying	vec3 g_viewvector;// camera pos - vertex pos 
varying vec2 g_texcoord;   

void main() 
{
    vec3 N = normalize( g_vertexnormal );   
#ifdef WITH_NORMALMAP 
    N = perturb_normal( N, g_viewvector, g_texcoord ); 
#endif   
    // ... 
}

Зелёная ось

OpenGL и DirectX считают началом координат начало массива данных текстуры. Текстурная координата (0, 0) - это угол пикселя, на который указывает указатель на данные изображения. Большинство 3д-пакетов устроены иначе: они считают началом координат нижний левый угол uv-развёртки. Это означает что начало текстурных координат находится в углу первого пикселя последней строки изображения.

Поиск по картинкам в Гугле показывает что нет никакого общепринятого соглашения о зелёном канале в картах нормалей. У некоторых зелёный канал указывает вверх, а у некоторых - вниз. Наши художники предпочитают направление вверх по двум причинам: во-первых, 3ds Max работает именно с таким форматом, а во-вторых, он выглядит естественнее, как при освещении зелёным светом сверху. Это облегчает восприятие карты нормалей невооружённым глазом.

Оглядываясь назад

Исходная статья была написана на уровне демонстрации идеи. Хотя алгоритм был протестирован и работал, он был тяжеловат для того времени. Сегодня ситуация изменилась. Я использую этот алгоритм в реальных проектах и больше не парюсь с касательными в вершинных атрибутах. Мне больше не важно, правильно ли экспортируются касательные из Макса или Майи. Художники тоже не чувствуют потери какой-то части своего пайплайна, потому что всё довольно естественно: есть геометрия, есть текстурные координаты и карта нормалей, и оно просто работает.

Нет никакого касательного базиса в нормалмаппинге. Касательный базис, содержащий нормаль, логически неверен. Под него маскируется двойственный ("кокасательный") базис, когда он ортогонален. А когда он становится неортогональным, он перестаёт работать. Используйте двойственный базис.

Список литературы

[1] James Blinn, “Simulation of wrinkled surfaces”, SIGGRAPH 1978

[2] Mark Peercy, John Airey, Brian Cabral, “Efficient Bump Mapping Hardware”, SIGGRAPH 1997

[3] Mark J Kilgard, “A Practical and Robust Bump-mapping Technique for Today’s GPUs”, GDC 2000

[4] Christian Schüler, “Normal Mapping without Precomputed Tangents”, ShaderX 5, Chapter 2.6, pp. 131 – 140

[5] Colin Barré-Brisebois and Stephen Hill, “Blending in Detail”

Вычисление касательного базиса налету +5