Расшифровка чёрных ящиков из habr-статьи «Процедурная растительность на OpenGL и GLSL».

На хабре есть прекрасная статья от 2016 года, рассказывающая как отрендерить на сцене динамическую растительность, красиво развивающуюся под действием «ветра». Там же, в самом начале статьи, приложена ссылка на github автора, и он крайне советует читать код по мере изучения материала статьи. Я же в данной статье постараюсь дать подробные комментарии к его шейдеру управления тесселяции (далее просто шейдер тесселяции) и приведу немного ссылок по тесселяции в OpenGL, которые могут оказаться полезными.

[spoiler title=’Бесполезные рассуждения’ style=’green’ collapse_link=’true’]Забегая в конец статьи — материал прекрасный. Много полезной (реально полезной) и новой информации, а особенно хорошо то, что есть пример кода и описание его работы от человека, который в этом, судя по всему, неплохо шарит (в отличии от меня). На шейдеры тесселяции в рунете статьей — по пальцам пересчитать, а тут такой подарок. Но есть и один минус — автор по ходу повествования приводит лишь примитивные части кода из своих шейдеров, или же те части, которые были потом им доработаны перед релизом. Тем, кто параллельно смотрит в код  — крайне сложно порой понимать, что же в них происходит (например я очень много путался с константах и не понимал, нафига же они нужны и как тут оказались). Как человек, никогда до этого не работавший с шейдерами тесселяции, но немного поработавший со стандартным для OpenGL версии 3.* вершинными и фрагментыми шейдерами, я не испытал сложностей с шейдером контроля тесселяции (TCS) — в ~35 строках кода сложно запутаться, да и плюсом этот момент автор отлично прокомментировал. Поэтому под разбор попал только шейдер управления тесселяцией (TES)[/spoiler]

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

• Перевод (как я понял) урока по тесселяции с сайта OGLDev с примитивными примерами и некоторым количеством кода. Урок внёс некоторую ясность в процесс побега вершин от вершинного шейдера по катакомбам тесселяции до фрагментного, но ясности не хватало.
• Далее был найден туториал от steps3d по шейдерам тесселяции. Общая картина стала налаживаться.
• Само собой некоторые аспекты объясняются в статье с habr, которую, собственно, и разбираем.
• Но самый подробный разбор со схемой был найден в книге Computer Graphics Through OpeNG, 2-ого издания. В ней приводится 4-ый участник тесселяционной оргии, о котором до этого не было сказано. Да и в целом визуализация написанного — всегда приятна.

Теперь же перейдём к исходному коду шейдера с полным разбором того, что и откуда берётся. Подразумевается, что вы прочитали статью и ознакомились в общих чертах с тем, как работает тесселяция с помощью OpenGL.

Как было сказано выше (в спойлере) шейдер контроля тесселяции проблем не вызывает. Поэтому начинаем с

шейдерного шейдера. Или шейдера оценки тесселяции, или же просто шейдера тесселяции. Но версия гугла мне нравится больше.

//выходной тип патча - изолинии с равным расстоянием
layout(isolines, equal_spacing) in;
//к слову, шейдер вызывается для каждой точки из группы Control Points
out vec3 normal; //будем выводить нормаль точки
out vec3 binormal; //бинормаль точки
out float stemThickness; //толщину стебля в данной координате 
out float along; //текущую координату вдоль стебля
flat out float stemIdx; //номер стебля в кусте
//flat тут означает, что значение не будет интерполировано

Далее уже веселее:

//включён ли ветер. 
//Если не передаём переменную через uniform с CPU, тогда по-умолчанию - включён
uniform float windEnabled = 1.0f; 
//включена ли турбулентность. К слову, турбулентность и ветер тут не имеют 
// никакого отношения к физике, просто автор подогнал параметры модели таким 
// образом, чтобы поведение модели было схоже с действием ветра и турбулентности 
// в реальной жизни
uniform float shiveringEnabled = 1.0f; 

//значение того, насколько сильно будут колыхаться стебли (или псевдо-турбулентность)
uniform float windTurbulence = 0.2; 

//значение скорости движения стеблей
uniform float windSpeed = 3000.0f; 

//направление ветра (в нашем случае дует исключительно по Оси X)
uniform vec3 windDirection = vec3(1.0f, 0.0f, 0.0f);

Дальше по-порядку методы:

// Классический 2D шум Перлина
float cnoise2dab(const in vec2 x, const in float a, const in float b) {
    return fma((cnoise2d(x) + 1.0f) * 0.5f, b-a, a);
}

Как писал автор в статье на habr, шум он не реализовывал и взял готовое решение у Стефана Густавсона. Как можно видеть, в функции cnoise2dab вызывается метод cnoize2d от аргумента x — он описан в файле classicnoize2D.glsl и, собственно, вычисляет шум.

Дальше:

// Вычисление толщины стебля в зависимости от координаты абстрактного патча,
// на которой мы в текущий момент находимся
float thickness(const in float x) {
    return (1.0f - x) / 0.9f;
}

Координаты абстрактного патчта нормализованы и располагаются в отрезке [0; 1]. Как я понимаю, координата из которой начинается стебель — находится в позиции 0. Соответственно по данной формуле у основания стебель будет наиболее широкий, а при движении вдоль стебля будет плавно уменьшать свою толщину.

Дальше:

float flexibility(const in float x) {
    return x * x;
}

Функция гибкости. Реализована для того, чтобы величину смещения текущей точки патча (а смещается она от ветра, турбулентности и шума) домножить на значение функции от аргумента — позиции текущей точки по оси x — т.е. вдоль стебля. Т.к. у основания позиция этой точки равняется 0 — основание не будет смещаться от ветра и прочих источников. Согласитесь, если стебель не рос из земли а «плавал» над ней — было бы немного странно.

Переходим к функции main().

// Coordinate along the stem for point position computtations.
stemThickness = thickness(gl_TessCoord.x);

//координата точки патча, для которой вызван шейдер
along = gl_TessCoord.x;

//номер стебля в кусте
stemIdx = gl_TessCoord.y;

//ID примитива. Вычисляется как отношение между ID патча (gl_PrimitiveID) к 
//количеству кустов на сцене (numPrimitives - передаётся в shader.glsl)
float primitiveIdx = gl_PrimitiveID / float(numPrimitives);

// Текущая координата вершины патча, для которой вызван шейдер
float t = gl_TessCoord.x;

Далее параметры, которые не описаны в статье и приходится понимать, откуда они и зачем:

//частота вращения точек патча относительно исходных 
//от эффекта турбулентности
//почему именно такая? я так думаю, что экспериментально
//поиграйтесь с ней в коде и поймёте, что это оптимальное значение
const float shiverRotFreq = 1526.0f;

// периодичность сдвига стеблей от эффекта турбулентности
// в зависимости с частотой, что мы ввели строчкой выше
const float shiverDriftFreq = shiverRotFreq / (2.0f * M_PI);

Далее интересно — берём случайное число, которое было загружено на GPU посредством текстуры. Как это делается — описано в статье с habr. Как это делается с пояснениями — описано мной тут.

//берём рандомный шум из текстуры с равномерно распределёнными случайными числами
float shiverPhase = texture(urandom01, gl_TessCoord.y).r * 2.0f * M_PI;

//радиус отклонения точки патча
const float shiverRadius = 0.1f;

// текущий номер кадра (если турбулентность вообще включена) + ID примитива
// и умноженное на период колебаний для текущей точки патча
// используется затем как вторая компонента вектора, переданного 
// в аргументы функции, генерирующей шум Перлина
float pIdxSfn = fma(sfn, shiveringEnabled, primitiveIdx) * shiverDriftFreq;
// так же для шума перлина, по факту принимает либо номер текущего кадра, 
// либо 0, если турбулентность отключена
float seSfn = sfn * shiveringEnabled;
//рандомный номер стебля
float randY = rand(gl_TessCoord.y);

Далее идёт тёмный лес, который автор позаимствовал у крутого парня, написавшего классический шум Перлина. Собственно, генерируется 3 вектора с шумами:

    vec2 v1 = vec2 (
        cnoise2dab( vec2((randY + 27.0f + seSfn) * shiverDriftFreq, pIdxSfn) , 0.0f, 0.5f),
        cnoise2dab( vec2((randY + 31.0f + seSfn) * shiverDriftFreq, pIdxSfn) , 0.5f, 1.0f)
    );

    vec2 v2 = vec2(
        cnoise2dab(vec2((randY + 99.0f + seSfn) * shiverDriftFreq, pIdxSfn), 0.5f, 1.0f),
        cnoise2dab(vec2((randY + 51.0f + seSfn) * shiverDriftFreq, pIdxSfn), 0.0f, 0.6f)
    );

    vec2 v3 = vec2(
        cnoise2dab(vec2((randY + 46.0f + seSfn) * shiverDriftFreq, pIdxSfn), 0.75f, 1.3f),
        cnoise2dab(vec2((randY + 76.0f + seSfn) * shiverDriftFreq, pIdxSfn), 0.75f, 1.3f)
    );

Теперь приступаем к тому, чтобы изогнуть нашу линию как кубическую кривую Безье. Автор приводит формулу, приведу и я её здесь, чтобы не надо было бегать по вкладкам:

И поехали в код. Так как мы работаем со скелетом, мы можем предположить (и тем самым очень облегчить себе жизнь), что скелет располагается в 2D-плоскости (нормали потом повернём отдельно). Таким образом все вектора для каждой из опорных точек P0 — P3  кривой Безье состоят из двух координат — x и y. t, в нашем случае — это координата текущей точки на абстрактном патче.

Т.к. наш куст растёт из земли, то координата точки P0 = (0;0). Остальные рассчитываются как интерполяция(mix) между вектором шума(v*), сдвигом битов у переменных частоты, периода и радиуса вращения стебля влево на значение шума (rotate) и между этими двумя параметрами при интерполировании встаёт переменная, определяющая включена ли турбулентность.

vec2 p0 = vec2(0.0f, 0.0f);
vec2 p1 = mix(v1, rotate(shiverRotFreq, shiverPhase, shiverRadius, v1), shiveringEnabled);
vec2 p2 = mix(v2, rotate(shiverRotFreq, shiverPhase, shiverRadius, v2), shiveringEnabled);
vec2 p3 = mix(v3, rotate(shiverRotFreq, shiverPhase, shiverRadius, v3), shiveringEnabled);

Далее деформируем изолинию с помощью формулы кривой Безье и смещаем радиус на рандомное значение, взятое из зашумленной текстуры:

const float stemOffset = 0.3f;
const float stemRandAmplitude = 1.0f;
float t1 = t - 1.0f;
vec3 position;

// Deform isoline as cubic bezier curve.
position.xy = -p0 * (t1 * t1 * t1) + p3 * (t * t * t) + p1 * t * (t1 * t1) * 3.0f -
    p2 * (t * t) * t1 * 3.0f;
// Offset by noisy circle radius.
position.x +=
    fma(texture(urandom01, gl_TessCoord.y + primitiveIdx).r, stemRandAmplitude, stemOffset);
position.z = 0.0f;

По комментариям автора, в принципе, всё понятно.

Затем идёт участок кода, хорошо описанный в самой статье — поворот нормали для текущей точки (зачем её поворачивать — описано там же 🙂 ):

float alpha = gl_TessCoord.y * 2.0f * M_PI;
float cosAlpha = cos(alpha);
float sinAlpha = sin(alpha);
mat3 circDistribution = mat3(
    cosAlpha, 0.0f, -sinAlpha,
    0.0f,     1.0f,      0.0f,
    sinAlpha, 0.0f, cosAlpha);
position = circDistribution * position;

normal = normalize(circDistribution * vec3(
    p0.y * (t1 * t1) * -3.0f + p1.y * (t1 * t1) * 3.0f - p2.y * (t * t) * 3.0f +
        p3.y * (t * t) * 3.0f - p2.y * t * t1 * 6.0f + p1.y * t * t1 * 6.0f,
    p0.x * (t1 * t1) *  3.0f - p1.x * (t1 * t1) * 3.0f + p2.x * (t * t) * 3.0f -
        p3.x * (t * t) * 3.0f + p2.x * t * t1 * 6.0f - p1.x * t * t1 * 6.0f,
    0.0f
));

 

Затем получаем бинормаль и переводим координаты в мировые. В самом конце инициализируем вектор позиции, который является выходным в TES:

binormal = (circDistribution * vec3(0.0f, 0.0f, 1.0f));

// переводим в мировые координаты
position += gl_in[0].gl_Position.xyz;

vec3 windState = position;
windState.x += sfn * windSpeed;
vec3 windDisplacement = windDirection *
    flexibility(t) * (cnoise3d(windState * windTurbulence) + 1.0f) * 0.5f;

gl_Position = vec4(position + windDisplacement * windEnabled, 1.0f);

Таким образом, более-менее код стал понятен.

Большое спасибо автору за этот пример и в целом за статью — она прекрасна.

Ещё раз продублирую на неё ссылку.