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

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

UPD. Картинка обновлена. Изменения описаны в конце статьи.



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

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:


Поглощение

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


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

Рассеивание

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


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

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


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

Что касается и : , где C — некая константа, которая может иметь разное значение для каждого канала в RGB, — плотность среды в точке .

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


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

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

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


где означает взятие интеграла по сфере, — фазовая функция, — свет пришедший из направления .

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


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


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

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


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


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


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

Небо

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


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

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

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

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

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


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


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

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

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


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



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



Результат:


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



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



Результат:


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

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

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

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

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

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

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


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

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


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

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

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


Для облаков нужно получить две специальных текстуры. Первая имеет размер 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 шум Ворлея с еще меньшим масштабом


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

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

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


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

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

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


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


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;
}

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

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

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



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


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


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


Я в итоге остановился на 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.


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

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 в , а ufmProperties.attenuation2 — C в . ufmProperties.sunIntensity — интенсивность излучения солнца. sunColor — цвет излучения солнца.

Результат:


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


Здесь следует сделать отступление. По каноничной модели освещения фазовая функция должна быть одна. Однако на деле получаемый результат никого не устраивает и фазовых функций все используют две, да еще и объединяют их значения особым образом. Я также остановился на двух фазовых функциях, но их значения просто складываю. Первая фазовая функция имеет 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

Результат:


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

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

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

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

Результат:


Осталось решить последнюю проблему. В реальной жизни чем более горизонтального направления придерживается взгляд, тем больше мы видим некий туман или дымку, которые не позволяют нам видеть совсем далекие объекты. Это также нужно отразить в коде. Я взял обычный косинус угла направления взгляда и экспоненциальную функцию. На основании этого вычисляется некий коэффициент 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 — для настройки вручную.


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

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 текстуры.


Curl-шум

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


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


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


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


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

Обновление. V2.

Второй итерацией я внес небольшие изменения и улучшил внешний вид облаков.

Мультирассеивание

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

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

Напомню формулу интегрирования при единственном рассеивании:


где

Формула для множественного рассеивания будет такой:


где

где a, b, c подбираются вручную, для начала можно взять значения 0.5 для всех, а N взять равным 7.

Мультирассеивание требуется только для облаков, для неба оно не внесет существенных изменений.

Исправление цвета неба

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

#define RAY_SCATTERING vec3(5.8 / 1000000.0, 13.5 / 1000000.0, 33.1 / 1000000.0)
#define MIE_SCATTERING vec3(2.0 / 1000000.0)
#define OZONE_SCATTERING vec3(0.0)
#define RAY_EXTINCTION RAY_SCATTERING
#define MIE_EXTINCTION (MIE_SCATTERING * 1.1)
#define OZONE_EXTINCTION (vec3(3.426, 8.298, 0.356) * 0.06 / 100000.0)

Результат

В принципе это все. Далее будут результаты того что получилось после внесения изменений. С мультирассеиванием без исправления неба:





С мультирасеиванием и с исправленным небом:

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

Презентации 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