Трассировка пути на GPU в Unity — часть 2

«Нет ничего хуже чёткого образа размытой концепции». – фотограф Энсел Адамс

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

Основываясь на уже имеющемся у нас коде, мы итеративно решим уравнение рендеринга, сформулированное Джеймсом Каджия в 1986 году и преобразуем наш рендерер в трассировщик пути, способный передавать вышеупомянутые эффекты. Мы снова будем использовать C# для скриптов и HLSL для шейдеров. Код выложен на Bitbucket.

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

Показанное ниже изображение отрендерено с помощью карты Graffiti Shelter с сайта HDRI Haven. Другие изображения в этой статье отрендерены с использованием карты Kiara 9 Dusk.

Уравнение рендеринга

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

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

Первая часть — называется двумерной функцией распределения отражательной способности объектов (bidirectional reflectance distribution function, BRDF). Эта функция визуально описывает тип материала, с которым мы имеем дело: металл или диэлектрик, тёмный или яркий, блестящий или матовый. BRDF определяет долю освещения, приходящего из , которая отражается в направлении . На практике это реализуется с помощью трёхкомпонентного вектора со значениями красного, зелёного и синего цвета в интервале .

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

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

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

---

¹Не забывайте об этом! Мы часто будем говорить о косинусе, и при этом всегда будем иметь в виду скалярное произведение. Так как , а мы имеем дело с направлениями (единичными векторами), то скалярное произведение и является косинусом в большинстве задач компьютерной графики.

На помощь приходит Монте-Карло

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

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

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

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

Предварительные условия

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

Накопление

По каким-то причинам Unity не передаёт мне HDR-текстуру как destination в OnRenderImage. У меня сработал формат R8G8B8A8_Typeless, поэтому точность быстро становится слишком низкой для накопления большого количества сэмплов. Чтобы справиться с этим, давайте добавим в скрипт на C# private RenderTexture _converged. Это будет наш буфер, накапливающий с высокой точностью результаты перед отображением их на экране. Текстуру мы инициализируем/освобождаем так же, как _target в функции InitRenderTexture. В функции Render удваиваем блиттинг:

Graphics.Blit(_target, _converged, _addMaterial);
Graphics.Blit(_converged, destination);

Детерминированные сцены

При внесении изменений в рендеринг для оценки эффекта полезным бывает сравнение с предыдущими результатами. Пока у нас при каждом перезапуске режима Play или рекомпиляции скрипта будет получаться новая случайная сцена. Чтобы избежать этого, добавим в скрипт на C# public int SphereSeed и следующую строку в начало SetUpScene:

Random.InitState(SphereSeed);

Теперь мы можем вручную задавать seed сцены. Введите любое число и отключайте/снова включайте компонент RayTracingMaster, пока не получите подходящую сцену.

Для изображений-примеров использовались следующие параметры: Sphere Seed 1223832719, Sphere Radius [5, 30], Spheres Max 10000, Sphere Placement Radius 100.

Случайность шейдера

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

float2 _Pixel;
float _Seed;

float rand()
{
    float result = frac(sin(_Seed / 100.0f * dot(_Pixel, float2(12.9898f, 78.233f))) * 43758.5453f);
    _Seed += 1.0f;
    return result;
}

Инициализируем _Pixel непосредственно в CSMain как _Pixel = id.xy, чтобы каждый пиксель мог использовать разные случайные значения. _Seed инициализируется из C# в функции SetShaderParameters.

RayTracingShader.SetFloat("_Seed", Random.value);

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

Сэмплирование полусферы

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

float3 SampleHemisphere(float3 normal)
{
    // Однородная выборка направлений полусферы
    float cosTheta = rand();
    float sinTheta = sqrt(max(0.0f, 1.0f - cosTheta * cosTheta));
    float phi = 2 * PI * rand();
    float3 tangentSpaceDir = float3(cos(phi) * sinTheta, sin(phi) * sinTheta, cosTheta);

    // Преобразование направления в пространство мира
    return mul(tangentSpaceDir, GetTangentSpace(normal));
}

Направления генерируются для полусферы, центрированной относительно положительной оси Z, поэтому нам необходимо преобразовать их, чтобы они центрировались относительно нужной нормали. Мы генерируем касательную и бинормаль (два вектора, ортогональных к нормали и ортогональных друг другу). Сначала мы выбираем вспомогательный вектор для генерирования касательной. Для этого мы берём положительную ось X, и возвращаемся к положительной Z только если нормально (приблизительно) сонаправлена с осью X. Затем мы можем использовать векторное произведение для генерирования касательной, а потом и бинормали.

float3x3 GetTangentSpace(float3 normal)
{
    // Выбираем вспомогательный вектор для векторного произведения
    float3 helper = float3(1, 0, 0);
    if (abs(normal.x) > 0.99f)
        helper = float3(0, 0, 1);

    // Генерируем векторы
    float3 tangent = normalize(cross(normal, helper));
    float3 binormal = normalize(cross(normal, tangent));
    return float3x3(tangent, binormal, normal);
}

Рассеяние по Ламберту

Теперь, когда у нас есть однородные случайные направления, мы можем приступить к реализации первой BRDF. Для рассеянного отражения наиболее часто используется BRDF по Ламберту, которая на удивление проста: , где — это albedo поверхности. Давайте вставим её в наше уравнение рендеринга Монте-Карло (я пока не буду учитывать излучаемость) и посмотрим, что произойдёт:

Давайте сразу вставим это уравнение в шейдер. В функции Shade заменим код внутри конструкции if (hit.distance < 1.#INF) следующими строками:

// Рассеянное затенение
ray.origin = hit.position + hit.normal * 0.001f;
ray.direction = SampleHemisphere(hit.normal);
ray.energy *= 2 * hit.albedo * sdot(hit.normal, ray.direction);
return 0.0f;

Новое направление отражённого луча определяется с помощью нашей функции однородных выборок по полусфере. Энергия луча умножается на соответствующую часть показанного выше уравнения. Так как поверхность не испускает никакого освещения (а только отражает свет, получаемый напрямую или косвенно от неба), мы возвращаем здесь 0. Не забывайте, что AddShader усредняет выборки, поэтому нам не нужно беспокоиться о . Функция CSMain уже содержит умножение на (следующий отражённый луч), поэтому нам остаётся не так много работы.

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

float sdot(float3 x, float3 y, float f = 1.0f)
{
    return saturate(dot(x, y) * f);
}

Давайте подведём итог тому, что пока делает наш код. CSMain генерирует первичные лучи камеры и вызывает Shade. При пересечении с поверхностью эта функция в свою очередь генерирует новый луч (однородно случайный в полусфере вокруг нормали) и учитывает BRDF материала и косинус в энергии луча. При пересечении луча с небом мы сэмплируем HDRI (наш единственный источник освещения) и возвращаем освещение, которое умножается на энергию луча (т.е. результат всех предыдущих пересечений, начиная с камеры). Это простой сэмпл, который смешивается со сходящимся результатом. В результате в каждом сэмпле учитывается воздействие .

Настало время проверить всё в работе. Так как у металлов нет рассеянного отражения, давайте пока их отключим в функции SetUpScene скрипта на C# (но всё равно вызывая здесь Random.value для сохранения детерминированности сцены):

bool metal = Random.value < 0.0f;

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

Зеркальное отражение по Фонгу

Вполне неплохо для всего нескольких строк кода (и небольшой доли математики). Давайте улучшим картинку, добавив зеркальные отражения с помощью BRDF по Фонгу. Исходная формулировка Фонга имела свои проблемы (отсутствие взаимоотражений и сохранения энергии), но к счастью другие люди устранили их. Усовершенствованная BRDF показана ниже. — это направление идеально отражённого света, а — это показатель Фонга, управляющий шероховатостью (roughness):

Интерактивный двухмерный график показывает, как выглядит BRDF по Фонгу при для луча, падающего под углом 45°. Попробуйте изменять значение .

Вставим это в наше уравнение рендеринга Монте-Карло:

И наконец сложим это с уже имеющейся BRDF по Ламберту:

И вот как они выглядят в коде вместе с рассеянием по Ламберту:

// Затенение по Фонгу
ray.origin = hit.position + hit.normal * 0.001f;
float3 reflected = reflect(ray.direction, hit.normal);
ray.direction = SampleHemisphere(hit.normal);
float3 diffuse = 2 * min(1.0f - hit.specular, hit.albedo);
float alpha = 15.0f;
float3 specular = hit.specular * (alpha + 2) * pow(sdot(ray.direction, reflected), alpha);
ray.energy *= (diffuse + specular) * sdot(hit.normal, ray.direction);
return 0.0f;

Заметьте, что мы заменили скалярное произведение на немного отличающееся, но эквивалентное (отражённое вместо ). Теперь снова включим металлические материалы в функции SetUpScene и проверим, как это работает.

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

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

Улучшенное сэмплирование

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

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

// Вычисляем вероятности рассеянного и зеркального отражений
hit.albedo = min(1.0f - hit.specular, hit.albedo);
float specChance = energy(hit.specular);
float diffChance = energy(hit.albedo);
float sum = specChance + diffChance;
specChance /= sum;
diffChance /= sum;

// Рулеткой выбираем путь луча
float roulette = rand();
if (roulette < specChance)
{
    // Зеркальное отражение
    ray.origin = hit.position + hit.normal * 0.001f;
    ray.direction = reflect(ray.direction, hit.normal);
    ray.energy *= (1.0f / specChance) * hit.specular * sdot(hit.normal, ray.direction);
}
else
{
    // Рассеянное отражение
    ray.origin = hit.position + hit.normal * 0.001f;
    ray.direction = SampleHemisphere(hit.normal);
    ray.energy *= (1.0f / diffChance) * 2 * hit.albedo * sdot(hit.normal, ray.direction);
}

return 0.0f;

energy является небольшой вспомогательной функцией, усредняющей цветовые каналы:

float energy(float3 color)
{
    return dot(color, 1.0f / 3.0f);
}

Так мы создали более красивый трассировщик лучей Уиттеда из предыдущей части, но теперь с настоящим рассеянным затенением (а это подразумевает мягкие тени, ambient occlusion, рассеянное глобальное освещение):

Выборка по важности

Давайте ещё раз взглянем на базовую формулу Монте-Карло:

Как видите, мы делим влияние каждой выборки (сэмпла) на вероятность выбора этого конкретного сэмпла. Пока мы использовали однородные выборки по полусфере, поэтому имели постоянную . Как мы видели выше, это далеко от оптимального, например, в случае BRDF по Фонгу, которая велика в очень малом количестве направлений.

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

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

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

Выборка косинусов

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

В этой статье Томаса Пуле описывается, как это делается. Мы добавим в нашу функцию SampleHemisphere параметр alpha. Функция определяет показатель выборки косинусов: 0 для однородной выборки, 1 для выборки косинусов, или выше для более высоких показателей Фонга. В коде это выглядит так:

float3 SampleHemisphere(float3 normal, float alpha)
{
    // Сэмплируем полусферу, где альфа определяет тип сэмплирования
    float cosTheta = pow(rand(), 1.0f / (alpha + 1.0f));
    float sinTheta = sqrt(1.0f - cosTheta * cosTheta);
    float phi = 2 * PI * rand();
    float3 tangentSpaceDir = float3(cos(phi) * sinTheta, sin(phi) * sinTheta, cosTheta);

    // Преобразуем направление в пространство мира
    return mul(tangentSpaceDir, GetTangentSpace(normal));
}

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

Ламбертова выборка по важности

Для начала мы усовершенствуем рендеринг рассеянного отражения. В нашем однородном распределении уже используется постоянная BRDF по Ламберту, но мы можем улучить её, добавив косинус. Распределение вероятностей выборки по косинусу (где ) равно , что упрощает нашу формулу Монте-Карло для рассеянного отражения:

// Рассеянное отражение
ray.origin = hit.position + hit.normal * 0.001f;
ray.direction = SampleHemisphere(hit.normal, 1.0f);
ray.energy *= (1.0f / diffChance) * hit.albedo;

Это немного ускорит наше рассеянное затенение. Теперь давайте займёмся настоящей проблемой.

Фонгова выборка по важности

Для BRDF по Фонгу процедура похожа. На этот раз у нас есть произведение двух косинусов: стандартного косинуса из уравнения рендеринга (как и в случае с рассеянным отражением), умноженного на собственный косинус BRDF. Мы займёмся только последним.

Давайте вставим распределение вероятностей из примеров выше в уравнение Фонга. Подробный вывод можно изучить в статье Lafortune and Willems: Using the Modified Phong Reflectance Model for Physically Based Rendering (1994):

// Зеркальное отражение
float alpha = 15.0f;
ray.origin = hit.position + hit.normal * 0.001f;
ray.direction = SampleHemisphere(reflect(ray.direction, hit.normal), alpha);
float f = (alpha + 2) / (alpha + 1);
ray.energy *= (1.0f / specChance) * hit.specular * sdot(hit.normal, ray.direction, f);

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

Материалы

Наконец, давайте расширим нашу генерацию сцены, чтобы создавать изменяющиеся значения гладкости и излучаемости сфер! В struct Sphere из скрипта на C# добавьте public float smoothness и public Vector3 emission. Так как мы изменили размер struct, нам нужно изменить шаг при создании Compute Buffer (4 × количество чисел float, помните?). Сделаем так, чтобы функция SetUpScene вставляла значения для гладкости и излучаемости.

В шейдере прибавьте обе переменные к struct Sphere и struct RayHit, а затем инициализируйте их в CreateRayHit. И, наконец, задайте оба значения в IntersectGroundPlane (жёстко заданной в коде, вставьте любые значения) и IntersectSphere (получающей значения из Sphere).

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

float SmoothnessToPhongAlpha(float s)
{
    return pow(1000.0f, s * s);
}

float alpha = SmoothnessToPhongAlpha(hit.smoothness);

Использование излучаемости выполняется возвратом значения в Shade:

return hit.emission;

Результаты

Глубоко вдохните. расслабьтесь и подождите, пока изображение превратится в такую красивую картину:

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

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

Благодарю за чтение! В третьей части мы (на какое-то время) выйдем из джунглей сэмплирования и затенения, и вернёмся в цивилизацию, чтобы встретиться с джентльменами Моллером и Трумбором. Нам нужно будет поговорить с ними о треугольниках.

Об авторе: Дэвид Кури — разработчик Three Eyed Games, программист Virtual Engineering Lab компании Volkswagen, исследователь компьютерной графики и музыкант в жанре Heavy Metal.

Автор: PatientZero

Источник