[Из песочницы] Пишем настоящий шум Перлина
По поисковому запросу шум перлина сразу попадается этот перевод на Хабре. Как справедливо заметили в комментариях к публикации, речь идёт вовсе не о шуме Перлина. Возможно, автор перевода и сам был не в курсе.
Чем выгодно отличается шум Перлина, легко заметить, если сравнить картинки.
Обычный шум (из той самой статьи): 
Шум Перлина: 
И увеличением количества октав первую картинку ко второй никак не приблизишь. Я не буду описывать достоинства шума Перлина и область его применения, а постараюсь объяснить как он реализован. Думаю, это будет полезно многим программистам, ведь хакерские исходники Кена Перлина не много объясняют даже при наличии комментариев.
Рассмотрим двухмерный вариант. Пока напишем только класс-заготовку. Входящие данные — двухмерный вектор или два числа с плавающей точкой: x, y.
Возвращаемое значение — число от -1.0 до 1.0:
public class Perlin2d
{
public float Noise(float x, float y)
{
throw new NotImplementedException();
}
}
Пару слов об интерполяции
Идея обычного сглаженного шума в том, что есть дискретная сетка псевдослучайных значений, и для запрашиваемой точки происходит интерполяция между узлами сетки (чем ближе точка к какому-нибудь узлу сетки, тем больше его значение соответствует значению узла).
Здесь в третьем условном квадрате точка в самом центре после интерполяции будет иметь значение 3:
Рассмотрим детальнее как там получается тройка. Координаты точки: x:2.5, y:0.5
Целочисленные координаты точки (верхний левый угол квадрата): x:2, y:0
получаются округлением в сторону 0 (функция floor).
Локальные координаты точки внутри квадрата получаются вычитанием: x = 2.5 – 2 = 0.5, y = 0.5 – 0 = 0.5
Берём значение левого верхнего угла квадрата (1) и верхнего правого (2). Интерполируем верхнюю грань используя локальную координату x (0.5). Линейная интерполяция выглядит так:
static float Lerp(float a, float b, float t)
{
// return a * t + b * (1 - t); можно переписать с одним умножением (раскрыть скобки, взять в другие скобки):
return a + (b - a) * t;
}
Берём значение левого нижнего угла квадрата (2) и нижнего правого (7). Интерполируем нижнюю грань используя всё ту же локальную координату x (0.5).
Результаты: верхняя: 1.5 нижняя: 4.5
Теперь осталась интерполяция верхней и нижней с использованием локальной координаты y (тоже 0.5): 1.5 * 0.5 + 4.5 * (1 – 0.5) = 3
Билинейная интерполяция самая простая, но и результат не самый привлекательный.
Другие варианты интерполяции подразумевают модифицирование локальной координаты (параметра t) перед интерполяцией. Получаются более плавные переходы возле граничных значений (0 и 1).
В шуме Перлина задействован первый вариант, он даёт достаточно сильное искривление.
static float QunticCurve(float t)
{
return t * t * t * (t * (t * 6 - 15) + 10);
}
...
// комбинирование с функцией линейной интерполяции:
Lerp(a, b, QuinticCurve(t))
Главная идея и отличие шума Перлина
Всё очень просто:
1. В узлах сетки — псевдослучайные вектора (двухмерные для двухмерного шума, трехмерные для трехмерного и так далее), а не псевдослучайные числа.
2. Интерполируем между скалярными произведениями a) векторов от вершин квадрата до точки внутри квадрата (куба в трехмерном варианте) и b) псевдослучайных векторов (при описании шума Перлина их называют градиентными векторами).
В своём улучшенном варианте шума Кен Перлин использует всего 12 градиентных векторов. Для двухмерного варианта требуется всего 4 — по количеству граней (у квадрата их 4). Вектора направлены (условно из центра куба/квадрата) в сторону каждой из граней и не нормализованы.
Вот они:
{ 1, 0 }
{ -1, 0 }
{ 0, 1 }
{ 0,-1 }
Итак, каждому узлу сетки соответствует один из четырёх векторов. Пусть вектор у нас будет массивом float-ов.
float[] GetPseudoRandomGradientVector(int x, int y)
{
int v = // псевдо-случайное число от 0 до 3 которое всегда неизменно при данных x и y
switch (v)
{
case 0: return new float[]{ 1, 0 };
case 1: return new float[]{ -1, 0 };
case 2: return new float[]{ 0, 1 };
default: return new float[]{ 0,-1 };
}
}
Реализация
Нам понадобится скалярное произведение векторов:
static float Dot(float[] a, float[] b)
{
return a[0] * b[0] + a[1] * b[1];
}
Главный метод:
public float Noise(float fx, float fy)
{
// сразу находим координаты левой верхней вершины квадрата
int left = (int)System.Math.Floor(fx);
int top = (int)System.Math.Floor(fy);
// а теперь локальные координаты точки внутри квадрата
float pointInQuadX = fx - left;
float pointInQuadY = fy - top;
// извлекаем градиентные векторы для всех вершин квадрата:
float[] topLeftGradient = GetPseudoRandomGradientVector(left, top );
float[] topRightGradient = GetPseudoRandomGradientVector(left+1, top );
float[] bottomLeftGradient = GetPseudoRandomGradientVector(left, top+1);
float[] bottomRightGradient = GetPseudoRandomGradientVector(left+1, top+1);
// вектора от вершин квадрата до точки внутри квадрата:
float[] distanceToTopLeft = new float[]{ pointInQuadX, pointInQuadY };
float[] distanceToTopRight = new float[]{ pointInQuadX-1, pointInQuadY };
float[] distanceToBottomLeft = new float[]{ pointInQuadX, pointInQuadY-1 };
float[] distanceToBottomRight = new float[]{ pointInQuadX-1, pointInQuadY-1 };
// считаем скалярные произведения между которыми будем интерполировать
/*
tx1--tx2
| |
bx1--bx2
*/
float tx1 = Dot(distanceToTopLeft, topLeftGradient);
float tx2 = Dot(distanceToTopRight, topRightGradient);
float bx1 = Dot(distanceToBottomLeft, bottomLeftGradient);
float bx2 = Dot(distanceToBottomRight, bottomRightGradient);
// готовим параметры интерполяции, чтобы она не была линейной:
pointInQuadX = QunticCurve(pointInQuadX);
pointInQuadY = QunticCurve(pointInQuadY);
// собственно, интерполяция:
float tx = Lerp(tx1, tx2, pointInQuadX);
float bx = Lerp(bx1, bx2, pointInQuadX);
float tb = Lerp(tx, bx, pointInQuadY);
// возвращаем результат:
return tb;
}
В качестве бонуса:
public float Noise(float fx, float fy, int octaves, float persistence = 0.5f)
{
float amplitude = 1; // сила применения шума к общей картине, будет уменьшаться с "мельчанием" шума
// как сильно уменьшаться - регулирует persistence
float max = 0; // необходимо для нормализации результата
float result = 0; // накопитель результата
while (octaves-- > 0)
{
max += amplitude;
result += Noise(fx, fy) * amplitude;
amplitude *= persistence;
fx *= 2; // удваиваем частоту шума (делаем его более мелким) с каждой октавой
fy *= 2;
}
return result/max;
}
И последнее — использование таблицы со случайными числами. В коде Кена Перлина такая таблица прописана вручную и достаются оттуда значения совсем по-другому. Здесь можно экспериментировать и от этого немало зависит равномерность шума и отсутствие в нём явных паттернов.
Я сделал
class Perlin2D
{
byte[] permutationTable;
public Perlin2D(int seed = 0)
{
var rand = new System.Random(seed);
permutationTable = new byte[1023];
rand.NextBytes(permutationTable); // заполняем случайными байтами
}
private float[] GetPseudoRandomGradientVector(int x, int y)
{
// хэш-функция с Простыми числами, обрезкой результата до размера массива со случайными байтами
int v = (int)(((x * 1836311903) ^ (y * 2971215073) + 4807526976) & 1023);
v = permutationTable[v]&3;
switch (v)
{
...
& 3 здесь обрезает любое int32 число до 3, читайте об операции AND на википедии
Операция типа % 3 тоже сработала бы, но намного медленней.
class Perlin2D
{
byte[] permutationTable;
public Perlin2D(int seed = 0)
{
var rand = new System.Random(seed);
permutationTable = new byte[1023];
rand.NextBytes(permutationTable);
}
private float[] GetPseudoRandomGradientVector(int x, int y)
{
int v = (int)(((x * 1836311903) ^ (y * 2971215073) + 4807526976) & 1023);
v = permutationTable[v]&3;
switch (v)
{
case 0: return new float[]{ 1, 0 };
case 1: return new float[]{ -1, 0 };
case 2: return new float[]{ 0, 1 };
default: return new float[]{ 0,-1 };
}
}
static float QunticCurve(float t)
{
return t * t * t * (t * (t * 6 - 15) + 10);
}
static float Lerp(float a, float b, float t)
{
return a + (b - a) * t;
}
static float Dot(float[] a, float[] b)
{
return a[0] * b[0] + a[1] * b[1];
}
public float Noise(float fx, float fy)
{
int left = (int)System.Math.Floor(fx);
int top = (int)System.Math.Floor(fy);
float pointInQuadX = fx - left;
float pointInQuadY = fy - top;
float[] topLeftGradient = GetPseudoRandomGradientVector(left, top );
float[] topRightGradient = GetPseudoRandomGradientVector(left+1, top );
float[] bottomLeftGradient = GetPseudoRandomGradientVector(left, top+1);
float[] bottomRightGradient = GetPseudoRandomGradientVector(left+1, top+1);
float[] distanceToTopLeft = new float[]{ pointInQuadX, pointInQuadY };
float[] distanceToTopRight = new float[]{ pointInQuadX-1, pointInQuadY };
float[] distanceToBottomLeft = new float[]{ pointInQuadX, pointInQuadY-1 };
float[] distanceToBottomRight = new float[]{ pointInQuadX-1, pointInQuadY-1 };
float tx1 = Dot(distanceToTopLeft, topLeftGradient);
float tx2 = Dot(distanceToTopRight, topRightGradient);
float bx1 = Dot(distanceToBottomLeft, bottomLeftGradient);
float bx2 = Dot(distanceToBottomRight, bottomRightGradient);
pointInQuadX = QunticCurve(pointInQuadX);
pointInQuadY = QunticCurve(pointInQuadY);
float tx = Lerp(tx1, tx2, pointInQuadX);
float bx = Lerp(bx1, bx2, pointInQuadX);
float tb = Lerp(tx, bx, pointInQuadY);
return tb;
}
public float Noise(float fx, float fy, int octaves, float persistence = 0.5f)
{
float amplitude = 1;
float max = 0;
float result = 0;
while (octaves-- > 0)
{
max += amplitude;
result += Noise(fx, fy) * amplitude;
amplitude *= persistence;
fx *= 2;
fy *= 2;
}
return result/max;
}
}
Результат: 
