Скорее всего, вам доводилось слышать о градиентном шуме, вернее, о той его версии, которая называется шум Перлина и описывает одну конкретную реализацию, сопряжённую с различными оптимизациями на уровне ЦП. Поскольку это невероятно мощный инструмент для творческой работы, он используется практически везде: при создании визуальных эффектов, видеоигр, процедурно-математического искусства и т.д. Да, как следует настроить его — порой тонкая работа, и неисправная реализация на первый взгляд всё равно может выглядеть хорошо или интересно. В конце концов, «смотрится неплохо, а я художник, я так вижу».  

Чтобы глубже и результативнее понять градиентный шум, мы сначала изучим его одномерную версию (в литературе этот случай обычно не рассматривается), а затем медленно пойдём вверх по лестнице измерений в сторону усложнения задачи. Эту тему мы будем рассматривать с точки зрения графического процессора (GPU), а не с точки зрения обычного ЦП. Все примеры кода и анимации, приведённые в этой статье, реализованы на WebGL2/GLSL (надеюсь, это будет не слишком сильно сказываться на производительности). Примеры должны работать на большинстве современных устройств.

Перед тем, как начать — необходимая оговорка: большая часть этого материала совершенно не нова. В этой статье изложен многонедельный опыт изучения соответствующего материала и сопутствующих экспериментов. Математическая основа взята с отличного сайта, который ведёт Иньиго Килес, а также из других разрозненных источников, имеющихся в Интернете. Но при всей ценности и насыщенности этих ресурсов, материал на них зачастую излагается слишком быстро — авторы пропускают многие детали, так как считают их очевидными. В данном посте я постараюсь заполнить эти пробелы.

Хеширующая функция и псевдослучайные значения

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

Реализация Перлина основана на таблице перестановок. Это удобно, если работать с обычным процессором (ЦП), но не так просто, если используется шейдер. При работе с графическим процессором (GPU) большинство практиков полагаются на различные ухищрения, при которых задействуются операции над числами с плавающей точкой или не самые оптимальные битовые операции. Часто такие приёмы не срабатывают при исследовании полного 32-разрядного диапазона входных значений.

Примечание

Мы не можем воспользоваться генератором псевдослучайных чисел, поскольку работа такого генератора опирается на состояние, в котором нам требуется детерминизм (для одной из координат нам всё время требуется случайное значение). Существует подкласс генераторов псевдослучайных чисел, так называемые LCG (линейные конгруэнтные генераторы), где состояние фактически содержится в возвращаемом значении. Но даже тогда мы не можем повторно использовать возвращённое значение, поскольку требуется поиск — искать будем случайное значение, ассоциированное с нашей целочисленной координатой. Если попытаемся скормить PRNG/LCG нашу координату вместо ожидаемого состояния, то, фактически, при каждом вызове мы будем засевать новое случайное значение, а это может привести к существенному сдвигу нормального распределения. Именно для того, чтобы этого избежать, нам требуется хешировать целые числа.

Итак, нам нужна хеширующая функция, а сами мы должны удовлетвориться 32 разрядами, так как работаем на графическом процессоре. К счастью, долго искать не приходится, поскольку ещё в 2018 году Крис Уэллонс  нашёл очень хорошую функцию для этой цели, назвав её lowbias32. Впоследствии эту функцию доработал TheIronBorn.

Таким образом, первый элемент конструкции, которым мы воспользуемся — это хеширующая функция:

uint hash(uint x) {
    x = (x ^ (x >> 16)) * 0x21f0aaadU;
    x = (x ^ (x >> 15)) * 0x735a2d97U;
    return x ^ (x >> 15);
}

Она сама по себе всем хороша, но от обычного 32-разрядного беззнакового целого как такового пользы немного. Нам понадобится нормированное число с плавающей точкой. Можно было бы, не мудрствуя, разделить на 0xffffffff, но тогда у нас получится неоднородное распределение (честно говоря, это не так важно, как кажется). Вместо этого приспособим применительно к числам с плавающей точкой  приём, продемонстрированный Виньей Себастьяно при работе с числами двойной точности и напишем код на GLSL:

float u2f(uint x) { return float(x >> 8U) * uintBitsToFloat(0x33800000U); }

Скомбинировав две эти функции в u2f(hash(x)), позволяющую сопоставить любую 32-разрядную координату x с числом с плавающей точкой в интервале [0,1).

Чтобы приспособить к этому примеру функцию h(x), рассмотренную выше, а также приблизить её по форме к «сигналу», можно центрировать её около  0, перестроив из интервала [0,1) в [−1,1): h(x)=u2f(hash(x))*2.0-1.0.

Расширяем хеширующую функцию на дополнительные измерения

При работе в 2 или более измерениях наша координатная сетка будет получать в качестве входных координат более 1 значения, но нам нужно придумать, как сообщать их нашей хеширующей функции, у которой всего один параметр. Один из способов этого добиться — использовать многоуровневый хеш, построенный по принципу «исключающего ИЛИ».

uint hash(uvec2 x) { return hash(x.x ^ hash(x.y)); }  // для ввода в 2D 
uint hash(uvec3 x) { return hash(x.x ^ hash(x.yz)); } // для ввода в 3D 

Воспользовавшись координатами пикселей в качестве ввода, можно проверить хеширующую функцию в 2D:

Предупреждение

На одном из этапов работы выполним промежуточное преобразование целого со знаком, а затем перейдём к операциям с беззнаковыми числами — это требуется, чтобы избежать проблем, связанных с отрицательными координатами. Как сказано в спецификации GLSL 4.60: «В случае преобразования отрицательного значения с плавающей точкой в uint возникает неопределённое поведение». При таком преобразовании проблемы у нас возникнут лишь в том  случае, когда координаты выйдут за пределы диапазона 32-разрядных знаковых чисел.

Базовый сигнал и белый шум

Кроме пошагового подхода также можно присваивать значения лишь целочисленным координатами и проводить между ними линейную интерполяцию. Попробуем сделать это в одномерном случае:

float value(int x) { return u2f(hash(uint(x)))*2.0 - 1.0; } // h(x)

float vnoise1_linear(float p) {
    int i = int(floor(p));                // целочисленные координаты или решётка
    float f = fract(p);                   // x-позиция между 2 окружающими значениями
    return mix(value(i), value(i+1), f);  // линейная интерполяция между 2 значениями
}

Задав этой функции в качестве значения координату x, получим амплитуду (высоту) сигнала:

Примечание

Возможно, это пока не очевидно, но красота рассматриваемого примера заключается в том, что этот «бесконечный» сигнал является детерминированным и, соответственно, поддаётся поиску. Действительно, можно продвинуться как вперёд, так и назад к любой реальной позиции p — и сразу же узнать, как в ней выглядит сигнал. Это принципиально важно для процедурного программирования.

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

• Кубический эрмитов сплайн f(t)=3t2−2t3 (также применяемый в функции GLSL smoothstep()). Именно им исходно воспользовался Кен Перлин, создавая первую реализацию шума Перлина.

• Более современную (и сложную) кривую пятого порядка f(t)=6t5−15t4+10t3 Кен Перлин предложил в 2002 году, чтобы скорректировать разрывы, возникающие в производных второго порядка f"(t).

Зафиксируем, что в GLSL:

float fade_quintic(float t) { return ((6.0*t-15.0)*t+10.0)*t*t*t; }
float fade_hermite(float t) { return (3.0-2.0*t)*t*t; }

В оставшейся части этой статьи продолжим работать с кривой пятого порядка:

#define fade fade_quintic

float vnoise1(float p) {
    int i = int(floor(p));
    float f = fract(p);
    float a = fade(f);
    return mix(value(i), value(i+1), a);
}

Одномерный градиентный шум

Всё-таки, сигнал, сгенерированный таким образом, может нам не подходить из-за резких изменений в кривизне/частоте; он слишком «нестабилен». Так что не будем использовать случайные значения непосредственно в качестве шума, а интерпретируем их как градиент: у нас получится градиентный шум.

Примечание

Педантично подходя к вопросу, признаем, что в 1D нельзя говорить о градиентах; следует рассматривать их как наклоны или угловые коэффициенты. Но всё равно будем считать их градиентами, чтобы обеспечить единообразие с ситуациями, возникающими в 2 и более измерениях.

В каждой решётке будем присваивать градиентам y=0  с регулярным интервалом, чтобы их случайные градиенты влияли на окружающую часть кривой. Может казаться, что это сложно реализовать на практике, но на самом деле мы обойдёмся 2 операциями умножения и 1 операцией вычитания дополнительно к получению численного шума (как минимум, в одномерном случае):

float grad(int x) { // целочисленная решётка для интервала [-1,1)
    return u2f(hash(uint(x))) * 2.0 - 1.0;
}

float noise1(float p) {
    int i = int(floor(p));
    float g0 = grad(i);
    float g1 = grad(i + 1);

    float f = fract(p);
    float v0 = g0 * f;
    float v1 = g1 * (f - 1.0);

    float a = fade(f);
    return mix(v0, v1, a);
}

С геометрической точки зрения, v0​ и v1​ — это координаты по оси y, которые мы получаем 2 наклона (маленькие палочки на каждой решётке) вокруг актуальной целевой точки p и находим, где их пересекает вертикальная линия x=p. Затем мы выполняем между ними гладкую интерполяцию, воспользовавшись нашей функцией затухания.

Примечание

У градиентного шума есть важное свойство: в каждой решётке его функция проходит через 0, что указывает на постоянство его частоты. Иными словами,  noise1​(p)=0 во всех случаях, когда p — целое число.

Расширяем пример до двух измерений

Если перенести этот пример в 2D, то в каждой точке решётки будет храниться 2-компонентный вектор градиентов. Чтобы получить значение шума в точке p=(x,y), не ограничиваясь простым умножением, вычисляем скалярное произведение каждого из градиентных векторов с вектором из угла решётки в (x,y), а затем билинейно интерполируем четыре этих скалярных произведения, пользуясь 2D-функцией затухания.

#define bmix(a,b,c,d,x,y) mix(mix(a,b,x),mix(c,d,x),y) // билинейная интерполяция

float noise2(vec2 p) {
    ivec2 i = ivec2(floor(p));
    vec2 g0 = grad(i);
    vec2 g1 = grad(i + ivec2(1, 0));
    vec2 g2 = grad(i + ivec2(0, 1));
    vec2 g3 = grad(i + ivec2(1, 1));

    vec2 f = fract(p);
    float v0 = dot(g0, f);
    float v1 = dot(g1, f - vec2(1.0, 0.0));
    float v2 = dot(g2, f - vec2(0.0, 1.0));
    float v3 = dot(g3, f - vec2(1.0, 1.0));

    vec2 a = fade(f);
    return bmix(v0, v1, v2, v3, a.x, a.y);
}

Поскольку теперь градиенты представлены в 2D, эту случайную градиентную функцию необходимо расширить. Нам нужен генератор нормированных векторов, то есть, механизм для создания векторов единичной длины, которые могли бы указывать в любом направлении.

Первый вариант решения — вызывать хеширующую функцию дважды (обычно применительно к самой себе), чтобы получать случайные координаты (x,y) в квадрате. Нормировав их, мы переносим их на круг:

vec2 grad(ivec2 x) { // решётка ivec2 для случайного нормированного 2D-вектора (нормированная точка квадрата)
    uint h1 = hash(uvec2(x));
    uint h2 = hash(h1);
    return normalize(vec2(u2f(h1), u2f(h2)) * 2.0 - 1.0);
}

Но может получиться даже лучше — при этом нам понадобится всего один хеш и немного тригонометрии:

const float TAU = 6.283185307179586;

vec2 grad(ivec2 x) { // решётка ivec2 для случайного нормированного 2D-вектора (точка круга)
    float angle = u2f(hash(uvec2(x))) * TAU;
    return vec2(cos(angle), sin(angle));
}

Этого достаточно, чтобы получить градиентный шум в 2D:

Расширяем пример до трёх измерений

Аналогично, чтобы построить градиентный шум в 3D, потребуется внести 2 следующих изменения:

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

• При использовании случайных единичных векторов в качестве градиентов в 3D их нужно равномерно распределить на сфере, а не на круге

vec3 grad(ivec3 x) { // решётка ivec3 для случайного нормированного единичного 3D-вектора (точка сферы)
    uint h0 = hash(uvec3(x));
    uint h1 = hash(h0);
    // используем первое случайное значение для полярного угла (широта)
    float c = 2.0*u2f(h0) - 1.0, // c = cos(theta) = cos(acos(2x-1)) = 2x-1
          s = sqrt(1.0 - c*c);   // s = sin(theta) = sin(acos(c)) = sqrt(1-c*c)
    float phi = TAU * u2f(h1);   // используем второе случайное значение для азимута (широта)
    return vec3(cos(phi) * s, sin(phi) * s, c);
}

#define tmix(a,b,c,d,e,f,g,h,x,y,z) mix(bmix(a,b,c,d,x,y),bmix(e,f,g,h,x,y),z) // трилинейная интерполяция

float noise3(vec3 p) {
    ivec3 i = ivec3(floor(p));
    vec3 g0 = grad(i);
    vec3 g1 = grad(i + ivec3(1, 0, 0));
    vec3 g2 = grad(i + ivec3(0, 1, 0));
    vec3 g3 = grad(i + ivec3(1, 1, 0));
    vec3 g4 = grad(i + ivec3(0, 0, 1));
    vec3 g5 = grad(i + ivec3(1, 0, 1));
    vec3 g6 = grad(i + ivec3(0, 1, 1));
    vec3 g7 = grad(i + ivec3(1, 1, 1));

    vec3 f = fract(p);
    float v0 = dot(g0, f);
    float v1 = dot(g1, f - vec3(1, 0, 0));
    float v2 = dot(g2, f - vec3(0, 1, 0));
    float v3 = dot(g3, f - vec3(1, 1, 0));
    float v4 = dot(g4, f - vec3(0, 0, 1));
    float v5 = dot(g5, f - vec3(1, 0, 1));
    float v6 = dot(g6, f - vec3(0, 1, 1));
    float v7 = dot(g7, f - vec3(1, 1, 1));

    vec3 a = fade(f);
    return tmix(v0, v1, v2, v3, v4, v5, v6, v7, a.x, a.y, a.z);
}

Примечание

Сферическое представление не имеет никакого отношения к выбору случайного единичного 3D-вектора на сфере. В данной реализации шум применяется к любой 3D-геометрии. В качестве примера мы выбрали сферу просто потому, что это самая простая трёхмерная фигура, на которой можно построить функцию. Если бы мы работали в 2D, то также могли бы сообщить функции трёхмерного шума аргумент p=(x,y,t), где t — текущее время. Так у нас отобразилась бы «дымчатая» атмосферная картинка.

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

Фрактальное броуновское движение (fBm)

Приём fBm основан на следующей идее: суммировать множество «октав» шума, чтобы на их основе собрать более чёткий узор. Как правило, для этого нужно поднять частоту, удвоив её (речь идёт о коэффициенте «лакунарности») и вдвое уменьшив амплитуду (этот коэффициент характеризует «усиление» или «устойчивость»):

Алгоритм практически одинаков во всех измерениях, это просто сумма сигналов. Например, в 2D можно записать:

const float LACUNARITY = 1.98;
const float GAIN = 0.51;

float fbm(vec2 p, int octaves) {
    float sum = 0.0;
    float amp = 1.0, freq = 1.0;
    for (int i = 0; i < octaves; i++) {
        sum += amp * noise2(p * freq);
        freq *= LACUNARITY;
        amp *= GAIN;
    }
    return sum;
}

Примечание

Не будем брать для лакунарности и усиления ровно 2,0 и 0,5 — это нужно, чтобы нарушить корреляции.

А в 3D для полноты картины (если не говорить о трассировщике лучей, при помощи которого мы описываем сферу, весь код шума остаётся прежним, он просто вызывает noise3(p*freq) и использует vec3 p в качестве ввода):

Производные

Производные (характеризующие темп изменения сигнала) бывают полезны во многих ситуациях. Например, вы уже могли заметить, что 1D-сигналы, рассматриваемые в этой статье, окрашены тем светлее, чем круче их наклон, и тем темнее, чем более он пологий. Для интерполяции между двумя цветами применяется производная.

Возвращаясь к случаю 1D: чтобы отобразить кривую с правильной толщиной линии, я пользуюсь производными, описанными у Иньиго, где он описывает фокус с расстоянием до кривой):

float dist = abs(v - p.y) / sqrt(1.0 + d*d); // v: значение кривой, p: позиция, d: производная кривой

Воспользовавшись такой величиной расстояния со знаком как любой другой, можем гладко отобразить кривую:

В более высоких измерениях можно воспользоваться производными, чтобы получить освещение: действительно, вооружившись производными, можно вычислить нормаль, которая затем будет использоваться для отражений:

// Ламбертов приём для работы с освещением
vec3 normal = normalize(vec3(-d, 1.0)); // нормаль от частичных производных d в 2D 
float a = 5.0*TAU/8.0; // свет падает с юго-запада
vec3 light_direction = normalize(vec3(cos(a), sin(a), 2.0));
float lighting = max(dot(normal, light_direction), 0.0);
col *= lighting;

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

float fbm2e(vec2 p) {
    float sum = 0.0;
    float amp = 1.0, freq = 1.0;
    vec2 d = vec2(0.0);
    for (int i = 0; i < octaves; i++) {
        vec3 n = noise2d(p * freq);    // скорректированная noise2(), возвращающая частичные производные в .xy
        d += n.xy;                     // накопленные производные без without frequency scaling
        float w = 1.0/(1.0+dot(d,d));  // вес на основе градиента/уклона 
        sum += amp * n.z * w;          // значение, убавленное с учётом веса
        freq *= LACUNARITY;
        amp *= GAIN;
    }
    return sum;
}

Интуитивно понятно, что по этой формуле каждый новый слой «приглушается», когда градиенты показывают крутой подъём. Так рисунок не будет излишне «изрезан» более высокими шумами, и весь рисунок будет смотреться резче.

В любом случае, это просто подборка некоторых вещей, которые можно сделать при помощи производных, но, вероятно, есть и другие, о которых я забыл. Суть в том, что они особенно полезны, поэтому дальше мы изучим производные подробнее.

Сравнение числовых и аналитических производных

Производные можно получать двумя способами:

• Числовой метод, пользуясь которым мы вычисляем скорость изменения между 2 близкими точками. Он «работает» с любыми кривыми (если мы располагаем достаточной прецизионностью для такой работы) но порой хромает в отношении точности (а иногда и в отношении скорости).

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

Числовой метод реализуем на GPU благодаря функциям thanks dFdx и dFdy. Эти функции используют данные локального фрагмента и соседних с ним фрагментов для вычисления производной указанного значения. Это одна из тех редких функций, которые сообщают информацию кроссфрагментно. Как вы уже догадываетесь, для работы с ней требуется синхронизация, что неизбежно влияет на производительность.

Но предположим, к примеру, что мы берём нашу последнюю сцену с 2D-шумом и сравниваем две производные: числовую и аналитическую:

Слева показан стандартный 2D-градиентный шум с фрактальным броуновским движением, а справа — производные шума: сначала аналитическая версия, а затем числовая. Вторая версия должна выглядеть зернистой.

Числовые производные

Числовые производные были получены следующим образом:

float w = 2.0/resolution.y;           // размер пикселя
float v = noise2(p);                  // значение шума в позиции p
vec2 d = vec2(dFdx(v), dFdy(v)) / w;  // частичные числовые производные

Примечание

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

Аналитические производные

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

Но нам также понадобятся производные интерполяционных функций для каждого измерения:

Вот же они:

Далее идёт доказательство частичных производных билинейных и трилинейных интерполяционных функций, и вы можете их пропустить, если вас не привлекает перспектива оказаться по колено в математических символах.

Примечание

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

или ещё более устрашающей версии с вертикальной чертой я воспользуюсь ∂xf(x=17,y=u+v,z=v2).

Наконец, имея эти частичные производные, можно выяснить, как вычислить производные самого шума. Полезно знать, что производная fract(p) (в большинстве случаев) равна 1, а производная floor(p) (обычно) равна 0, поэтому здесь всё достаточно аккуратно упрощается. На этот раз избавлю вас от подробностей:

Скорректируем наши функции градиентного шума, чтобы они возвращали производные вместе со значением шума как таковым:

vec2 noise1d(float p) {
    int i = int(floor(p));
    float g0 = grad(i);
    float g1 = grad(i + 1);

    float f = fract(p);
    float v0 = g0 * f;
    float v1 = g1 * (f - 1.0);

    float a = fade(f);
    float v = mix(v0, v1, a);

    float g = mix(g0, g1, a);
    float d = v1 - v0;                    // производная mix относительно интерполянта 
    float da = ((f-2.0)*f+1.0)*30.0*f*f;  // fade'(t): производная интерполянта пятой степени

    return vec2(g + d*da, v);
}
vec3 noise2d(vec2 p) {
    // [...]
    float v = bmix(v0, v1, v2, v3, a.x, a.y);

    vec2 g = bmix(g0, g1, g2, g3, a.x, a.y);
    vec2 d = mix(vec2(v1,v2)-v0, v3-vec2(v2,v1), a.yx); // производные bmix относительно интерполянта 
    vec2 da = ((f-2.0)*f+1.0)*30.0*f*f;

    return vec3(g + d*da, v);
}
vec4 noise3d(vec3 p) {
    // [...]
    float v = tmix(v0, v1, v2, v3, v4, v5, v6, v7, a.x, a.y, a.z);

    vec3 g = tmix(g0, g1, g2, g3, g4, g5, g6, g7, a.x, a.y, a.z);
    vec3 d = bmix(vec3(v1,v2,v4)-v0, vec3(v3-v2,v3-v1,v5-v1),  // производные tmix 
                  vec3(v5-v4,v6-v4,v6-v2), v7-vec3(v6,v5,v3),  // относительно интерполянта
                  a.yxx, a.zzy);
    vec3 da = ((f-2.0)*f+1.0)*30.0*f*f;

    return vec4(g + d*da, v);
}

Примечание

Производные находятся в первой компоненте, так что, имея n=noise3d(p), можем записать vec3 d=n.xyz для частичных производных x/y/z, а не прибегать к более неуклюжим vec3 d=n.yzw.

Чтобы получить эквивалентные производные для численного шума, просто устанавливаем  g=0 в окончательном выражении производных.

Можно размотать вложенное выражение mix, чтобы (попробовать) ускорить его работу, но я считаю mix-выражения простыми, элегантными и, пожалуй, более стабильными в числовом отношении.

Использование производных

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

float freq = 0.1;
vec4 n = noise3d(x*freq);
float v = n.a;            // значение
vec2 d = n.xyz * freq;    // частичные производные

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

Аналогично, предположим, что хотим переместить v из интервала  [−1,1] в [0,1]:

v = (v + 1.0) / 2.0; // [-1,1] -> [0,1]
d = d / 2.0;

Интуитивно понятно, что, если ужать значение вполовину от его исходного размера, то его производные (наклоны) также станут более пологими.

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

vec4 fbm3d(vec3 p) {
    vec4 sum = vec4(0.0);
    float amp = 1.0, freq = 1.0;
    for (int i = 0; i < octaves; i++) {
        vec4 n = noise3d(p * freq);
        sum.xyz += amp * n.xyz * freq; // производные
        sum.a += amp * n.a; // значение
        freq *= LACUNARITY;
        amp *= GAIN;
    }
    return sum;
}

Всё усложняется, когда шум корректируется с учётом самих производных, как мы видели выше в примере с эрозией, особенно потому, что при этом будут применяться производные второго порядка.

Аналогично, существует техника с поворотом p на каждой итерации fBm, что призвано снизить корреляцию между слоями шума. Такой поворот представляет собой линейное преобразование, но требует и тщательных корректировок в сами производные. В качестве домашнего задания предложу вам самим реализовать эти преобразования правильно. Напомню: очень важен строгий учёт производных.

Что дальше

В этой вводной статье мы успели исследовать лишь верхушку айсберга. Так, вам будет интересно узнать, что существуют и альтернативные шумы, например, OpenSimplex, где происходит интерполяция с позициями решётки на симплексной (в 2D получаются треугольники, в 3D — тетраэдры), а не на прямоугольной решётке. У этого шума есть свои достоинства — например, при fBm он даёт меньше артефактов при выборе направления, тем самым избавляя нас от необходимости делать ситуативные повороты.

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

Возможно, вам также будет интересно исследовать, как согласованно масштабировать шум, чтобы он помещался в контролируемом интервале значений. Спойлер: это на редкость сложная тема.

Кроме того, сейчас мы остановились на 3D, но и 4D изучить полезно. Например, мы хотим изменить форму 3D-шума с учётом времени:  p=(x,y,z,t).