Реализация физически корректных объемных облаков как в игре Horizon Zero Dawn

Раньше облака в играх рисовались обычными 2D спрайтами, которые всегда повернуты в направлении камеры, но последние годы новые модели видеокарт позволяют рисовать физически корректные облака без заметных потерь в производительности. Считается, что объемные облака в игры принесла студия Guerrilla Games вместе с игрой Horizon Zero Dawn. Конечно, такие облака умели рендерить и раньше, но студия сформировала что-то вроде промышленного стандарта на исходные ресурсы и используемые алгоритмы, и в настоящее время любая реализация объемных облаков так или иначе этому стандарту соответствует.

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

Tone mapping, sRGB


Перед началом работы с освещением важно сделать две вещи:

  1. Перед выводом финального изображения на экран применить хотя бы простейший tone mapping:
    tunedColor=color/(1+color)
    Это необходимо потому что вычисляемые значения цветов будут намного больше единицы.
  2. Убедиться что конечный фреймбуффер, в который вы рисуете, и который будет представлен на экране, имеет формат sRGB. Если с активацией sRGB режима проблемы, преобразование можно сделать вручную в шейдере:
    finalColor=pow(color, vec3(1.0/2.2))
    Формула подойдет для большинства случаев, но не в 100% в зависимости от монитора. Важно, что sRGB преобразование всегда делается в последнюю очередь.


Модель освещения


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

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

_11t9wizftnmwilerlsi5fvcrm4.png


Поглощение

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

56734cec6154e1fab8bb2d1538a8dffa.svg


где a3e20411f145906908c5cd2b53a36dbe.svg — оставшийся после поглощения свет в точке 5ccc8d80d010abd90b5c546d05856d3a.svg. 5ccc8d80d010abd90b5c546d05856d3a.svg — точка на отрезке AB на расстоянии 2a31024ea1803c34a47496e24a53a1ef.svg от A.

Рассеивание

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

04bc172f53fa572aebf0399275b5009f.svg


где 5e30b7badd64b9da62be750e48573c7d.svg — доля света, не изменившего свое направление, после рассеивания в точке 5ccc8d80d010abd90b5c546d05856d3a.svg.

Поглощение и рассеивание необходимо объединить:

8d7d5af1a0506ee776c4c2f850e862a9.svg


Функцию 33f2974e9c84331f8b92d8b7ed5aca97.svg называют затуханием (attenuation или extinction). А функцию eead4d0d4f08873f9aa49b69a7c46fd8.svg — передаточной функцией. Она показывает сколько света останется при прохождении из точки A в точку B.

Что касается 610ddbee8d80f847b44fa75c9a136489.svg и cd18c240521dcfb71084a340a504e41f.svg: b76ea9fd577eb1c2cee545caa4153bd3.svg, где C — некая константа, которая может иметь разное значение для каждого канала в RGB, 9ee4949d2897f56ac90e563db1eb96cd.svg — плотность среды в точке 5ccc8d80d010abd90b5c546d05856d3a.svg.

Теперь усложним задачу. Свет движется из в точки A в точку B, он затухает в процессе движения. В точке X часть света рассеивается в разных направлениях, одно из направлений соответствует наблюдателю в точке О. Далее порция рассеянного света движется из точки X в точку О и снова затухает. Путь интересующего нас света AXO.

zob_troftvd8-ggp1ga7wesnrlk.png


Потерю света при движении из A в X мы знаем: 25d8c9f76986f59eb045b801f19d4e0a.svg, также как и знаем потерю света из X в O — это 25b3a98e370f0a95b734794b74bee57a.svg. Но что на счет доли света, которая будет рассеяна в направлении наблюдателя?

Усиливающее рассеивание

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

7e4be1bd9a1d8f05b9d1d47e1bb27ecc.svg


где e638781e456c920d7864e85c8b14a746.svg означает взятие интеграла по сфере, 1ad0bb23eb22871d78d2b5cc67d58292.svg — фазовая функция, da97bc3c51ddad728182deecd39ba0d0.svg — свет пришедший из направления 9c6a44e8ed4cb10046e00346c6b58b15.svg.

Посчитать свет со всех направлений достаточно сложно, однако, мы знаем, что основную порцию света несет наш оригинальный луч AB. Формулу можно существенно упростить:

ef84aa22259c648ea9925502712250f6.svg


где 2cbbcb347a44c276c1095ac5bb3f8242.svg — угол между лучом света и лучом наблюдателя (т.е. угол AXO), 3091ee38b2de66a8aa89af86fe87eb65.svg — исходное значение интенсивности света. Суммируя все вышесказанное, получим формулу:

642fb727356c245fc6f053dbc9affde8.svg


где 3091ee38b2de66a8aa89af86fe87eb65.svg — входящий свет, f7afbd138efabfc5bf5dad3e26e6c0fb.svg — свет, дошедший до наблюдателя.

Усложним задачу еще чуть-чуть. Допустим, свет испускается источником света типа directional light, т.е. солнцем:

wzr2ggl_1oshhmhn3kkqalhtkoo.png


Происходит все тоже самое как в предыдущем случае, но много раз. Свет из точки А1 рассеивается в точке Х1 в сторону наблюдателя в точке О, свет из точки А2 рассеивается в точке Х2 в сторону наблюдателя в точке О и т.д. Мы видим, что свет достигший наблюдателя равен сумме:

83819124f0beb6bdd5adc3644a79653b.svg


Или более точное интегральное выражение:

e017b05d2026f311f43bc53da48c130e.svg


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

Небо


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

k0mmtzjcir3wnqxbcjzsouln4yo.png


И даже не одному типу рассеивания, а целым двум: релеевскому рассеиванию и рассеиванию Ми. Первое вызвано молекулами воздуха, а второе аэрозолем воды.

Суммарная плотность воздуха (либо аэрозоля), через которую пройдет луч света, двигаясь из точки A в точку B:

0bcd80e5d5b35e43fef8336e6c66cde0.svg

где 14b01abf043545d56ac10d61e69792e8.svg — масштабирующая высота, h — текущая высота.

Простое решение интеграла будет таким:

d3fb07e750a91c2a253fdcf1de74de72.svg


где dh — размер шага, с которым берется выборка высоты.

Теперь посмотрим на рисунок и воспользуемся формулой, выведенной в предыдущей части «модель освещения»:

grn0fr0f-dly-jjpcjpg_pxof9a.png


Наблюдатель смотрит из O в O». Мы хотим собрать весь свет, который достигнет точек X1, X2, …, Xn, будет в них рассеян, а затем достигнет наблюдателя:

540c25470466173df2e5cf7880390ff9.svg


где c87cb3ca4e0a0f0fc9d1a9d5c485c116.svg интенсивность излучаемого солнцем света, 6c3b25b5f22b9a03c093ef2f25a06339.svg — высота в точке ee62f8bed8d0d0e9b3d3cd2c81841ded.svg; в случае с небом константа С, находящаяся в функции cd18c240521dcfb71084a340a504e41f.svg обозначается как ceaf36742def9a68dbd173a0ac9c9aeb.svg.

Решение интеграла может быть таким:

5be44a9fcdb92469fee29b924a957556.svg


Эта формула справедлива как для релеевского рассеивания так и для рассеивания Ми. В итоге значения света для каждого из рассеиваний просто складываются:

34fcb4a9100e4b9c60a7c7eca2c489ec.svg


Релеевское рассеивание

71ea55e8043a70776eb70f29eb84d877.svg

67903f8e987346a525fb22e02e718aad.svg (содержит значения для каждого канала RGB)

4c902f3848ae3cb24347ee0b367132d6.svg

Результат:

nkxs8zt8hyh909tjilktrkrrvri.png


Рассеивание Ми

72b94bcc19f99b0e164c6a48e4408b44.svg

c9a52e00c5e6169a8ad778f902e49595.svg (значения для всех каналов RGB одинаковы)

208fcf3bf216540af2bfa93a186a3a19.svg

Результат:

fujvsbtduq-il-1reylizx1kosa.png


Количество выборок на отрезке 2a1abf7c35ec6c9a9d35514977d15981.svg и на отрезке 392559f75e1c271b8bbfbc6255e70d4d.svg можно брать 32 и выше. Значения радиуса Земли — 6371000 м, толщина атмосферы — 100000 м.

Что со всем этим делать:

  1. В каждом пикселе экрана вычисляем направление наблюдателя V
  2. Берем положение наблюдателя O равным {0, 6371000, 0}
  3. Находим 28dcc9fe786e9216afc4ab3ae2196ebe.svg как результат пересечения луча, берущим начало в точке O, и направлением V и сферы с центром в точке {0,0,0} и радиусом 6471000
  4. Отрезок 392559f75e1c271b8bbfbc6255e70d4d.svg делим на 32 участка равной длины
  5. Для каждого участка вычисляем релеевское рассеивание и рассеивание Ми, и все складываем. При этом для вычисления 2a1abf7c35ec6c9a9d35514977d15981.svg нам также также потребуется разделить отрезок 0e6fb04e8525f5a3448df984993ab3c4.svg на 32 равных участка в каждом случае. 60e970a20b30f36e7c21fe71c3d4dadb.svg можно считать через переменную, значение которой увеличивается на каждом шаге в цикле.


Итоговый результат:

nogbtlqbuzffixmq7xmt1w4j_ja.png


Модель облаков


Нам потребуется несколько видов шума в 3D. Первый — это тайлящийся fractal Brownian motion (fBm) шум Перлина:

Результат для 2D среза:

ldh6v6zglre4jhon_n8qu193o8k.png


Второй — тайлящийся fBm шум Вороного.

Результат для 2D среза:

aarcs2hvles3h89jqfc1lfrufjk.png


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

float fbmTiledWorley3(...)
{
    return clamp((1.0-fbmTiledVoronoi3(...))*1.5-0.25, 0.0, 1.0);
}


Результат сразу напоминает облачные структуры:

hapfemxrqrxuqfpckajp53j4zxc.png


Для облаков нужно получить две специальных текстуры. Первая имеет размер 128×128х128 и отвечает за низкочастотный шум, вторая имеет размер 32×32х32 и отвечает за высокочастотный шум. Каждая текстура использует только один канал в формате R8. В некоторых примерах используют 4 канала R8G8B8A8 для первой текстуры и три канала R8G8B8 для второй, а затем каналы смешивают в шейдере. Я не вижу в этом смысла, потому что смешивание можно произвести заранее тем самым получив большее попадание в cache coherence.

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

float remap(float value, float minValue, float maxValue, float newMinValue, float newMaxValue)
{
    return newMinValue+(value-minValue)/(maxValue-minValue)*(newMaxValue-newMinValue);
}


Начнем подготовку текстуры с низкочастотным шумом:
R-канал — тайлящийся fBm шум Перлина
G-канал — тайлящийся fBm шум Ворлея
B-канал — тайлящийся fBm шум Ворлея с меньшим масштабом
A-канал — тайлящийся fBm шум Ворлея с еще меньшим масштабом

ms_9z-hofvyvya5-zuh6vpi6da0.png


Смешивание производим таким образом:

finalValue=remap(noise.x, (noise.y * 0.625 + noise.z*0.25 + noise.w * 0.125)-1, 1, 0, 1)


Результат для 2D среза:

64qhrl3wl2vrsbedao4fk9zr3xg.png


Теперь готовим текстуру с высокочастотным шумом:
R-канал — тайлящийся fBm шум Ворлея
G-канал — тайлящийся fBm шум Ворлея с меньшим масштабом
B-канал — тайлящийся fBm шум Ворлея с еще меньшим масштабом

oxm6inosejxwrkx2kdz7zv17eea.png
finalValue=noise.x * 0.625 + noise.y*0.25 + noise.z * 0.125;


Результат для 2D среза:

lgmrmc2ohi3ayicwbji49fdhqb4.png


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

fuoim8wbtdzgqoc2lma8g6gb0ow.png


R-канал — покрытие облаков на малой высоте
G-канал — покрытие облаков на большой высоте
B-канал — максимальная высота облаков
A-канал — плотность облаков

Теперь мы готовы к тому, чтобы составить функцию, которая будет возвращать плотность облаков в зависимости от координат 3D пространства.

На входе точка пространства с координатами в км

vec3 position

Сразу добавляем смещение ветром

position.xz+=vec2(0.2f)*ufmParams.time;

Получаем значения погодной карты

vec4 weather=textureLod(ufmWeatherMap, position.xz/4096.0f, 0);

Получаем процент высоты (от 0 до 1)

float height=cloudGetHeight(position);

Добавляем небольшое скругление облаков снизу:

float SRb=clamp(remap(height, 0, 0.07, 0, 1), 0, 1);

Делаем линейное уменьшение плотности до 0 с увеличением высоты согласно B-каналу погодной карты:

float SRt=clamp(remap(height, weather.b*0.2, weather.b, 1, 0), 0, 1);

Объединяем результат:

float SA=SRb*SRt;

Снова добавляем скругление облаков снизу:

float DRb=height*clamp(remap(height, 0, 0.15, 0, 1), 0, 1);

Также добавляем скругление облаков сверху:

float DRt=height*clamp(remap(height, 0.9, 1, 1, 0), 0, 1);

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

float DA=DRb*DRt*weather.a*2*ufmProperties.density;

Объединяем низкочастотный и высокочастотный шум из наших текстур:

float SNsample=textureLod(ufmLowFreqNoiseTexture, position/48.0f, 0).x*0.85f+textureLod(ufmHighFreqNoiseTexture, position/4.8f, 0).x*0.15f; 

Во всех документах, которые я читал, объединение происходит другим способом, но мне понравился такой вариант.

Определяем величину покрытия (% неба, занятый облаками), которое задается через gui, также используются R-, G- каналы погодной карты:

float WMc=max(weather.r, clamp(ufmProperties.coverage-0.5, 0, 1)*weather.g*2);

Вычисляем финальную плотность:

float d=clamp(remap(SNsample*SA, 1-ufmProperties.coverage*WMc, 1, 0, 1), 0, 1)*DA;


Функция целиком:

float cloudSampleDensity(vec3 position)
{
	position.xz+=vec2(0.2f)*ufmParams.time;

	vec4 weather=textureLod(ufmWeatherMap, position.xz/4096.0f+vec2(0.2, 0.1), 0);
	float height=cloudGetHeight(position);
	
	float SRb=clamp(remap(height, 0, 0.07, 0, 1), 0, 1);
	float SRt=clamp(remap(height, weather.b*0.2, weather.b, 1, 0), 0, 1);
	float SA=SRb*SRt;
	
	float DRb=height*clamp(remap(height, 0, 0.15, 0, 1), 0, 1);
	float DRt=height*clamp(remap(height, 0.9, 1, 1, 0), 0, 1);
	float DA=DRb*DRt*weather.a*2*ufmProperties.density;
	
	float SNsample=textureLod(ufmLowFreqNoiseTexture, position/48.0f, 0).x*0.85f+textureLod(ufmHighFreqNoiseTexture, position/4.8f, 0).x*0.15f; 
	
	float WMc=max(weather.r, clamp(ufmProperties.coverage-0.5, 0, 1)*weather.g*2);
	float d=clamp(remap(SNsample*SA, 1-ufmProperties.coverage*WMc, 1, 0, 1), 0, 1)*DA;
	
	return d;
}


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

tfnzbzvd2oucw2h_oed9pbm_8ku.png

Интегрирование


Атмосфера Земли делится на два слоя: внутренний и внешний, между которыми могут находится облака. Эти слои можно представить сферами, а можно и плоскостями. Я остановился на сферах. Для первого слоя я взял радиус сферы 6415 км, для второго слоя радиус 6435 км. Радиус Земли округлил до 6400 км. От условной толщины «облачной» части атмосферы (20 км) в дальнейшем будут зависеть некоторые параметры.

0zir8r9qtvwlh9xk85i-hznhqa0.png


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

5kneljw4unaofsxxgbq12ealebw.png


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

edsrlyis3-2or1nztj6q-osmats.png


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

5prrdfocofvzsjobvawnairijp0.png


Я в итоге остановился на 4 выборках, которые лежат на одной линии, однако последняя берется с увеличенным в 6 раз шагом. Размер шага 20 км * 0.01, что равно 200 м.

Функция довольно проста:

float cloudSampleDirectDensity(vec3 position, vec3 sunDir)
{
	//определяем размер шага
	float avrStep=(6435.0-6415.0)*0.01;
	float sumDensity=0.0;
	for(int i=0;i<4;i++)
	{
		float step=avrStep;
		//для последней выборки умножаем шаг на 6
		if(i==3)
			step=step*6.0;
		//обновляем позицию
		position+=sunDir*step;
		//получаем значение плотности, вызывая функцию, которая уже 
                //рассматривалась ранее
		float density=cloudSampleDensity(position)*step;
		sumDensity+=density;
	}
	return sumDensity;
}


Теперь можно перейти к более сложной части. Определяем наблюдателя на поверхности Земли в точке {0, 6400,0} и находим пересечение луча наблюдения со сферой радиусом 6415 км и центром {0,0,0} — получаем стартовую точку S.

25jezjeea4j3bl51jneypfjetsi.png


Ниже базовый вариант функции:

vec4 mainMarching(vec3 viewDir, vec3 sunDir)
{
	vec3 position;
	crossRaySphereOutFar(vec3(0.0, 6400.0, 0.0), viewDir, vec3(0.0), 6415.0, position);
	
	float avrStep=(6435.0-6415.0)/64.0;
	
	for(int i=0;i<128;i++)
	{
		position+=viewDir*step;
		if(length(position)>6435.0)
			break;
	}
	
	return vec4(0.0);
}


Размер шага определяем как 20 км / 64. Т.е. в случае со строго вертикальным направлением луча наблюдателя мы сделаем 64 выборки. Однако, когда это направление более горизонтально, то выборок будет несколько больше, поэтому в цикле не 64 шага, а 128 — с запасом.

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

vec3 color=vec3(0.0);
float transmittance=1.0;

…
//это все внутри цикла

//получаем плотность облака для текущей выборки
float density=cloudSampleDensity(position)*avrStep;
//получаем суммарную плотность облака, через которую
//прошел луч солнца
float sunDensity=cloudSampleDirectDensity(position, sunDir);

//по знакомой формуле находим значение света				
float m2=exp(-ufmProperties.attenuation*sunDensity);
float m3=ufmProperties.attenuation2*density;
float light=ufmProperties.sunIntensity*m2*m3;

//добавляем найденный цвет и обновляем значение прозрачности
color+=sunColor*light*transmittance;
transmittance*=exp(-ufmProperties.attenuation*density);

…

return vec4(color, 1.0-transmittance);


ufmProperties.attenuation — есть не что иное как C в 610ddbee8d80f847b44fa75c9a136489.svg, а ufmProperties.attenuation2 — C в cd18c240521dcfb71084a340a504e41f.svg. ufmProperties.sunIntensity — интенсивность излучения солнца. sunColor — цвет излучения солнца.

Результат:

vcftz9vsswkxpqoggirezpfpygm.png


В глаза сразу бросается недостаток — сильное затенение. Но сейчас мы будем исправлять отсутствие усиленного освещения вблизи солнца. Так получилось, потому что мы не добавили фазовую функцию. Для расчета рассеивания света, проходящего через облака, используется фазовая функция Хеньи-Гринштейна, которые открыли ее в 1941 году для аналогичных вычислений в газовых скоплениях в космосе:

4b0b2e33fbd07036b14317d40cb1eac0.svg


Здесь следует сделать отступление. По каноничной модели освещения фазовая функция должна быть одна. Однако на деле получаемый результат никого не устраивает и фазовых функций все используют две, да еще и объединяют их значения особым образом. Я также остановился на двух фазовых функциях, но их значения просто складываю. Первая фазовая функция имеет g близкое к 1 и позволяет сделать яркое освещение вблизи солнца. Вторая фазовая функция имеет g близкое к 0.5 и позволяет сделать постепенное уменьшение освещенности по всей небесной сфере.

Обновленный код:

//находим cos(theta)
float mu=max(0, dot(viewDir, sunDir));
				
float m11=ufmProperties.phaseInfluence*cloudPhaseFunction(mu, ufmProperties.eccentrisy);
float m12=ufmProperties.phaseInfluence2*cloudPhaseFunction(mu, ufmProperties.eccentrisy2);
float m2=exp(-ufmProperties.attenuation*sunDensity);
float m3=ufmProperties.attenuation2*density;
float light=ufmProperties.sunIntensity*(m11+m12)*m2*m3;


ufmProperties.eccentrisy, ufmProperties.eccentrisy2 — это значения g

Результат:

ipqrfivyafyxnmdswnn_nbpvfdk.png


Теперь можно начать борьбу со слишком сильным затенением. Оно присутствует, потому что мы не учли свет от окружающих облаков и неба, который есть в реальной жизни.

Я решил эту проблему так:

return vec4(color+ambientColor*ufmProperties.ambient, 1.0-transmittance);


Где ambientColor — цвет неба в направлении луча наблюдения, ufmProperties.ambient — настроечный параметр.

Результат:

ct8hjiyk7brz5pv8vuwjeji33bc.png


Теперь нужно решить последнюю проблему. В реальной жизни чем более горизонтального направления придерживается взгляд, тем больше мы видим некий туман или дымку, которые не позволяют нам видеть совсем далекие объекты. Это также нужно отразить в коде. Я взял обычный косинус угла направления взгляда и экспоненциальную функцию. На основании этого вычисляется некий коэффициент blending, который позволяет сделать линейную интерполяцию между результирующим цветом и цветом фона.

float blending=1.0-exp(-max(0.0, dot(viewDir, vec3(0.0,1.0,0.0)))*ufmProperties.fog);
blending=blending*blending*blending;
return vec4(mix(ambientColor, color+ambientColor*ufmProperties.ambient, blending), 1.0-transmittance);


ufmProperties.fog — для настройки вручную.

02td9ho6kelgum68wldvcf-zymu.png

Итоговая функция:

vec4 mainMarching(vec3 viewDir, vec3 sunDir, vec3 sunColor, vec3 ambientColor)
{
	vec3 position;
	crossRaySphereOutFar(vec3(0.0, 6400.0, 0.0), viewDir, vec3(0.0), 6415.0, position);
	
	float avrStep=(6435.0-6415.0)/64.0;
	
	vec3 color=vec3(0.0);
	float transmittance=1.0;
	
	for(int i=0;i<128;i++)
	{
		float density=cloudSampleDensity(position)*avrStep;
		if(density>0.0)
		{
			float sunDensity=cloudSampleDirectDensity(position, sunDir);
			float mu=max(0.0, dot(viewDir, sunDir));
				
			float m11=ufmProperties.phaseInfluence*cloudPhaseFunction(mu, ufmProperties.eccentrisy);
			float m12=ufmProperties.phaseInfluence2*cloudPhaseFunction(mu, ufmProperties.eccentrisy2);
			float m2=exp(-ufmProperties.attenuation*sunDensity);
			float m3=ufmProperties.attenuation2*density;
			float light=ufmProperties.sunIntensity*(m11+m12)*m2*m3;
		
			color+=sunColor*light*transmittance;
			transmittance*=exp(-ufmProperties.attenuation*density);
		}
		position+=viewDir*avrStep;

		if(transmittance<0.05 || length(position)>6435.0)
			break;
	}
	
	float blending=1.0-exp(-max(0.0, dot(viewDir, vec3(0.0,1.0,0.0)))*ufmProperties.fog);
	blending=blending*blending*blending;
	return vec4(mix(ambientColor, color+ambientColor*ufmProperties.ambient, blending), 1.0-transmittance);
}


Видео с демонстрацией:


Оптимизация и возможные улучшения


После реализации базового алгоритма рендеринга следующая проблема заключается в том, что он слишком медленно работает. Мой вариант выдавал 25 fps в full hd на radeon rx 480. Два следующих подхода для решения проблемы предложили сами Guerrilla Games.

Рисуем то что действительно видно

Экран делится на тайлы размером 16×16 пикселей. Сначала отрисовывается обычное 3D-шное окружение. Оказывается, что большая часть неба закрыта горами или крупными объектами. Соответственно требуется выполнить расчет только в тех тайлах, облака в которых ничем не загорожены.

Репроекция

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

Частичное обновление

Идея с репроекцией мне не нравится тем, что при резком повороте камеры может получиться, что рендерить облака придется для трети экрана, что может вызвать лаг. Не знаю, как с этим разобрались Guerrilla Games, но по крайней мере, в Horizon Zero Dawn при управлении джойстиком камера движется плавно и проблем с резкими скачками там нет. Поэтому в качестве эксперимента я придумал свой подход. Облака рисуются в кубическую карту, в 5 граней, т.к. низ нас не интересует. Сторона кубической карты имеет пониженное разрешение, равное ⅔ от высоты экрана. Каждая грань кубической карты делится на тайлы размером 8×8. Каждый кадр на каждой грани обновляется только один из 64 пикселей в каждом тайле. Это дает заметные артефакты при резких изменениях, но т.к. Облака довольно статичны, то такой трюк незаметен. В итоге radeon rx 480 выдает 500 fps в full hd под вулкан и 330 fps под opengl. Radeon hd 5700 series выдает 109 fps в full hd под opengl (vulkan не поддерживает).

Использование mip-уровней

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

Высотные облака

Чтобы симулировать наличие перисто-высотных и перисто-кучевых облаков в Guerrilla Games при интеграции последние выборки делаются не из 3D текстур, о которых я рассказывал, а из специальной 2D текстуры.

udf2vjpfht0hhpbf9zptwfmadhg.jpeg


Curl-шум

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

-nqqpgmixyogyjxqrdyhm4oxlpc.png


Божественные лучи

qdwqxvzveiwuxgu2v84aw5fhjwo.png


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

apvj-assr02-vaov3o1k5frfawk.png


Теперь нужно применить радиальное сглаживание.

jdlg1hp7jajzhlbom4ocmn3vl70.png


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

Полезные ссылки


Презентации Guerrilla Games
d1z4o56rleaq4j.cloudfront.net/downloads/assets/Nubis-Authoring-Realtime-Volumetric-Cloudscapes-with-the-Decima-Engine-Final.pdf? mtime=20170807141817
killzone.dl.playstation.net/killzone/horizonzerodawn/presentations/Siggraph15_Schneider_Real-Time_Volumetric_Cloudscapes_of_Horizon_Zero_Dawn.pdf
www.youtube.com/watch? v=-d8qT5–1LOI

Реализация в GPU Pro 7
vk.com/doc179245989_437393482? hash=a9af5f665eda4edf58&dl=806d4dbdac0f7a761c

Небо
www.scratchapixel.com/lessons/procedural-generation-virtual-worlds/simulating-sky/simulating-colors-of-the-sky

Подход Frostbite
media.contentapi.ea.com/content/dam/eacom/frostbite/files/s2016-pbs-frostbite-sky-clouds-new.pdf
www.shadertoy.com/view/XlBSRz

© Habrahabr.ru