Доброго времени суток. Не так давно мне удалось реализовать один довольно интересный алгоритм затенения, о чем я и хочу рассказать.

И так, есть система частиц, находящаяся на освещаемой сцене.

Необходимо добавить затенение от частиц и по возможности самозатенение самих частиц.

Есть две основных проблемы, не позволяющие использовать классические карты теней. Первая из них — это тот факт, что тень отбрасывается полупрозрачным объектом, а значит и сама тень не будет полной во многих затеняемых точках. Вторая проблема заключается в том, что частицы поглощают свет на протяжении всей глубины рисуемого объекта (дым, облако, пыль), а не поверхностью. А это значит, что при расчете самозатенения необходимо учитывать изменение светового потока внутри самой системы частиц.

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

Описание алгоритма

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

А сейчас разберемся, что представляет из себя искомая функция.

Для удобства будем считать, что мы имитируем дым. Мы можем представить его как огромное множество полностью непрозрачных мелких частиц, распределенных по объему. Здесь под частицами я понимаю именно физические микрочастицы, а не метод отрисовки. Луч света, проходящий через этот объем, может либо наткнуться на частицу и полностью поглотиться, либо беспрепятственно пройти через этот объем. Причем этот процесс имеет случайный характер. То есть в качестве функции освещенности мы можем принять вероятность беспрепятственного прохождения света через систему частиц до искомой точки. Или, другими словами, в качестве функции затенения мы можем принять вероятность того, что луч от источника света пройдет меньшее расстояние, чем расстояние до освещаемой точки. Отлично! Нам нужна функция распределения случайной величины, и мы знаем как её получить.

Итак, метод базируется на Variance Shadow Maps (можно почитать тут и тут). Конкретнее нас интересует оценка случайной величины через её два первых момента.

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

Пока что выглядит сложно и запутанно. Будем разбираться поэтапно.

Построение карты теней

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

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

Итак, рассмотрим рендер в один единственный тексель теневой карты. Пусть на отрисовку попала одна единственная частица, и пусть значение непрозрачности (opacity) будет a.

Обозначим через d расстояние от источника света до спрайта частицы, а через L – расстояние, на котором происходит поглощение света. Тут надо понимать, что в данном конкретном случае свет будет всегда поглощаться на расстоянии d(или не поглощаться вовсе), но вообще-то L – это случайная величина, распределение которой нам и нужно будет оценить.

Доля лучей, поглощенных частицей, составляет a. Доля лучей, прошедших беспрепятственно, составляет b = 1 – a.

Спрайт расположен перпендикулярно направлению камеры, поэтому среднее расстояние, на котором произошло поглощение M[L] = d, а среднее значение квадрата этого расстояния равно M[L²] = d².

Немного усложним модель, вспомнив, что частицы имитируют не плоскую фигуру, а некоторый объем. Пусть свет равномерно поглощается на участке, равном размеру частицы. Размер частицы обозначим через s.

В этом случае изменится только M[L²]

Теперь рассмотрим случай, когда на рендер попали две частицы.

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

Обозначим непрозрачность спрайтов через a₁ и a₂, а расстояния от источника света до спрайтов через d₁ и d₂.

Доля света, прошедшего через первый спрайт и добравшегося до второго спрайта, равна b₁ = 1 - a_1. А доля света, прошедшего через оба спрайта, равна

Средние значения расстояний(и их квадратов), на которых происходит поглощение света для каждого из спрайтов по отдельности, вычисляются так же как и в предыдущем случае. Обозначим их M₁[L], M₁[L²] и M₂[L], M₂[L²] для первого и второго спрайтов соответственно.

Так как доля света, поглощенного первым спрайтом, равна a₁; а доля света, поглощенного вторым спрайтом, равна (1 – a₁) * a₂, то средние значения для обеих частиц будут равны:

Обозначим a₁ + (1 – a₁) * a₂ как w и умножим на него левые и правые части полученных равенств:

Обратите внимание, если мы рисуем сначала частицу №2, а потом поверх неё рисуем частицу №1, то величины M[L]*w, M[L²] * w и a₁ + (1 – a₁) * a₂ могут быть получены с помощью альфа блендинга. Таким образом мы можем накапливать эти значения в карте теней, а после по ним восстановить M[L] и M[L²].

Также стоит обратить внимание, что w = 1 – b₀. То есть w – это доля света, которая была поглощена обеими частицами.

Случай, когда на отрисовку попадают 3 и более частиц, не отличается от описанного выше. Достаточно отсортировать частицы от дальней к ближней.

Таким образом, сортируя частицы и используя альфа блендинг мы можем накапливать в карте теней значения M[L]*w, M[L²]*w и w.

Затенение

Итак, мы находимся на этапе расчета освещения. Я не буду здесь касаться материалов и тонкостей работы с источниками освещения. Описываемый алгоритм занимается только тем, что определяет, какая доля светового потока достигла некоторой точки в пространстве.

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

Обозначим через h расстояние от затеняемой точки до источника света.

Первое, что нам необходимо, это прочитать значения из карты теней. Напомню, что там мы накопили значения M[L]*w, M[L²]*w и w, где w – это доля света, поглощенная частицами на всем пути прохождения света; а L – это длина светового луча до его поглощения частицами.

Вычислим дисперсию:

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

Не сильно мудрствуя, я использовал аппроксимацию функцией smoothstep на отрезке в 6 сигм:

От искомого коэффициента затенения нас отделяет только тот факт, что часть света может беспрепятственно проходить через всю систему частиц. Его доля равна 1 – w. В этом случае коэффициент затенения будет равен:

Особенности реализации

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

Прежде всего, я использовал API Vulkan, поэтому вся терминология будет взята из него. Шейдерный язык, естественно, GLSL.

Метод используется в дополнение к классическим картам теней. То есть тени от непрозрачной геометрии рассчитываются отдельно от теней частиц. Тем не менее используется общая карта теней в формате VK_FORMAT_R16G16B16A16_UNORM. Канал R используется для теней от непрозрачных объектов, GBA каналы используются для теней от системы частиц. Так как используется формат с фиксированной точкой, то все длины приводятся в диапазон [0;1].

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

Построение карты теней происходит в два прохода. На первом проходе в канал R отрисовывается вся непрозрачная геометрия. На втором проходе в каналы BGA отрисовываются частицы.

Параметры блендинга для второго прохода:

  VkPipelineColorBlendAttachmentState blendingState
  {
    VK_TRUE,                              // blendEnable
    VK_BLEND_FACTOR_SRC_ALPHA,            // srcColorBlendFactor
    VK_BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,  // dstColorBlendFactor
    VK_BLEND_OP_ADD,                      // colorBlendOp
    VK_BLEND_FACTOR_ONE,                  // srcAlphaBlendFactor
    VK_BLEND_FACTOR_ONE_MINUS_SRC_ALPHA,  // dstAlphaBlendFactor
    VK_BLEND_OP_ADD,                      // alphaBlendOp
    VK_COLOR_COMPONENT_G_BIT |            // colorWriteMask;
      VK_COLOR_COMPONENT_B_BIT |
      VK_COLOR_COMPONENT_A_BIT
  };

Средние значения дистанций до источника света вычисляются в геометрическом шейдере:

float distance = -gl_in[0].gl_Position.z;
distance -= nearFar.near;				// Приводим дистанцию и размер
distance /= nearFar.nearFar;			// К диапазону [0;1]
float normalizedSize = size / nearFar.nearFar;

// Вычисляем среднее значение квадрата дистанции поглащения
float nearDistance3 = distance - normalizedSize / 2.f;
nearDistance3 = nearDistance3 * nearDistance3 * nearDistance3;
float farDistance3 = distance + normalizedSize / 2.f;
farDistance3 = farDistance3 * farDistance3 * farDistance3;
float sqDistance2 = (farDistance3 - nearDistance3) / 3.f / normalizedSize;

vec2 avgDistances = vec2(distance, sqDistance2);

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

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

float getShadowFactor(int layer,
                      vec2 shadowCoords,
                      float normalizedDistanceToLight,
                      float slope)
{
  …
  float opaqueFactor = getOpaqueShadowFactor(layer,
                                             centerCoords,
                                             normalizedDistanceToLight,
                                             slope);

  float transparentFactor = getTransparentShadowFactor( layer,
                                                        centerCoords,
                                                        normalizedDistanceToLight);

  return min(opaqueFactor, transparentFactor);
}

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

Важным здесь является то, что мы сначала вычисляем тень от непрозрачных объектов(getOpaqueShadowFactor), а потом от частиц (getTransparentShadowFactor) и выбираем из них наиболее темную.

Код getTransparentShadowFactor приведу полностью:

float getTransparentShadowFactor( int layer,
                                    vec2 centerCoords,
                                    float normalizedDistanceToLight)
{
  vec3 variadicValues = textureLod(shadowMap[layer], centerCoords, 0).gba;
  if(variadicValues.z == 0.f) return 1.f;

  float avgDistance = variadicValues.x / variadicValues.z;
  float avgSqDistance = variadicValues.y / variadicValues.z;
  float variance2 = max(avgSqDistance - avgDistance * avgDistance, 0.f);

  float deviation = normalizedDistanceToLight - avgDistance;
  float limit = 3.f * sqrt(variance2);
  float blackout = smoothstep(-limit, limit, deviation);

  blackout *= variadicValues.z;

  return 1.f - blackout;
}

Здесь variance2 – это дисперсия. Обратите внимание, что из-за ошибок дискретизации иногда дисперсия получается отрицательной величиной, поэтому мы вынуждены использовать функцию max. deviation – это отклонение дистанции до источника от среднего расстояния поглощения.

Результаты

Прежде всего соберем тестовую сцену. Нам понадобится источник света и объекты, на которые будет отбрасываться тень.

Модель баобаба я взял из opengameart.org, большое спасибо YouriNikolai. Модель распространяется по лицензии CC-BY-SA 4.0. Использована оригинальная геометрия, текстуры заменены на монотонные материалы.

На сцену добавлен источник направленного освещения, затенение с использованием shadow maps. Размер карты теней 2048*2048. Также на сцене присутствует рассеянный монотонный свет от фона.

Все тесты проводились на NVIDIA GeForce RTX 3060 Ti. Размер поверхности рендера 1615*962 пикселя. Мультисэмплинг не используется.

Заключение

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