Представленное в этой статье исследование было проведено исследователем графики и лектором Мэттью Делером в DAE Research (HOWEST University of Applied Sciences). Результаты были использованы для получения степени магистра в Breda University of Applied Sciences под научным руководством кандидата технических наук Жако Биккера. Эта статья является сокращённой версией исследования и может служить введением в тему. Она содержит примеры кода, вкратце объясняющие алгоритм, что, как мы надеемся, поможет разработчикам и исследователям в воссоздании прототипа в собственном фреймворке/движке. Публикацию/дипломную работу можно найти на ResearchGate, а полный пример кода — по этой ссылке.

Контекст

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


Симуляция глобального освещения часто основывается на уравнении рендеринга, предложенном Каджийей [1]. Хотя это интегральное уравнение описывает в общем виде распространение света, само оно нерешаемо, поскольку бесконечно рекурсивно. Несмотря на невозможность аналитического решения уравнения, симуляцию распространения света можно аппроксимировать при помощи методик интегрирования Монте-Карло. Интегрирование Монте-Карло может создавать реалистичные изображения, но часто оно вычислительно затратно. Также качество изображения зависит от количества сделанных сэмплов и распределения этих сэмплов. Хотя множество студий улучшало степень сходимости методов Монте-Карло, например, при помощи таких методик, как importance sampling и next event estimation, для обеспечения хорошей аппроксимации по-прежнему необходимо много сэмплов. На оборудовании текущего поколения невозможно обеспечивать интерактивные частоты кадров в приложениях реального времени при получении множества сэмплов на пиксель (samples per pixel, SPP) для вычисления отражённого освещения.

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

При офлайн-рендеринге эту проблему часто решают иным способом. Вместо сохранения информации о затенении в экранные буферы она сохранятся в 3D-точки, находящиеся в виртуальной сцене [2]. Другими словами, информация сохраняется в ближайшей 3D-точке, являющейся частью облака 3D-точек в мировом пространстве. Это разрешает некоторые из ранее описанных проблем. При этом нет необходимости использовать данные движения и данные можно сохранять в постоянной структуре данных. Хотя это очень сильная методика, она часто требует этапа предварительной обработки для создания набора 3D-точек, а этот набор точек занимает много памяти.

Последние несколько лет исследователи изучали способы переноса этой методики в приложения реального времени, чтобы объём занимаемой ею памяти был при этом минимальным. EA SEED [3] реализовала методику, использующую гибридный конвейер рендеринга для создания во время выполнения приложения облака точек. При помощи алгоритма заполнения пробелов в экранном пространстве создаются surfels (surface elements, элементы поверхностей). Несмотря на чрезвычайную перспективность этой методики, она всё равно требует дополнительную информацию, относящуюся к самому сохраняемому surfel (например, позицию, радиус, нормаль к поверхности и т.п.).

В нашем исследовании мы стремимся отказаться от сохранения этой дополнительной информации при помощи схемы хэширования, реализованной Nvidia в методике path space filtering [4], [5]. Методика Nvidia, по сути, выполняет косвенную вокселизацию с фиксированной дискретизацией, а мы воспользуемся результатами вокселизации для создания функции случайного шума, генерирующей 3D-точки. Таким образом мы попытаемся совместить мощь случайных наборов точек с представленной Nvidia методикой хэширования для снижения объёма занимаемой памяти и отказа от затратных операций поиска в памяти, положившись на имплицитность. По сравнению с методикой Nvidia, наша методика не использует фиксированную дискретизацию, снижая таким образом количество потенциальных артефактов дискретизации. Затем мы воссоздаём готовое изображение фильтрацией кэшированных данных из структуры данных, которая также является хэш-таблицей, и делаем это в мировом пространстве, а не используем фильтры воссоздания экранного пространства.

В следующем разделе представлено краткое описание алгоритма вместе с кодом шейдера. Полный прототип, разработанный на основе фреймворка Microsoft MiniEngine, выложен на Github.

Алгоритм

Алгоритм состоит из четырёх этапов. Во время Implicit Point Generation Pass (прохода генерации косвенных точек) для каждого пикселя определяется ближайшая случайная косвенная 3D-точка. Описывающие точку данные (т.е. уникальный ID) сохраняются во временный 2D-буфер. В процессе определения ближайшей точки параллельно можно вычислять информацию о затенении. В прототипе мы используем трассировку лучей с аппаратным ускорением для вычисления ambient occlusion (AO), которое обладает схожими с отражённым освещением свойствами. После нахождения ближайшей точки для каждого пикселя мы накапливаем всю информацию затенения во временной хэш-таблице на этапе Accumulation Pass (прохода накопления). Это необходимо, потому что на одну косвенную 3D-точку своей информацией затенения может влиять множество пикселей. При накоплении мы используем атомарные операции (атомарное сложение), гарантирующие правильность сложения всей информации. После этого прохода у нас получается небольшая хэш-таблица, в которой каждый элемент уникален. На этапе World Caching Pass (прохода кэширования мира) эта хэш-таблица объединяется с постоянной хэш-таблицей. Так как каждый элемент уникален, атомарные операции не требуются. После обновления постоянной структуры данных мы можем воссоздать готовое изображение на этапе Visualization Pass (прохода визуализации), отфильтровав наши данные в мире при помощи модифицированной интерполяции Шепарда [6].


Различные проходы в основном необходимы для обеспечения правильной синхронизации меду данными в GPU, а само ядро алгоритма довольно просто. 3D-локация из сцены, воссоздаваемая из буфера глубин, дискретизируется и хэшируется, по сути, превращаясь из 3D-представления в уникальное 1D-представление. Этот уникальный ID затем используется в качестве порождающего значения (seed value) для генератора случайного 3D-шума, снова создавая 3D-локацию или точку. Затем эта точка хэшируется для создания уникального 1D-индекса, который используется для поиска элемента в постоянной хэш-таблице. Мы пропускаем этот последний этап, потому что вместо этого используем seed value для поиска уникального элемента в хэш-таблице. Разработчик может легко добавить этот этап хэширования, чтобы проверить, соответствует ли он реализации хэш-таблицы, использованной в его фреймворке. Также будет полезно добавить дополнительное хэширование для повышения занятости хэш-таблицы, что сильно влияет на общую производительность.

Генерация косвенных точек

Первая задача этого прохода — дискретизация 3D-локации, которая нам в данный момент интересна. Как сказано, ранее, 3D-локацию пикселя можно получить, воссоздав её из буфера глубин. У Мэтта Петтинео есть очень хорошая статья об этом. После дискретизации локации мы используем её для создания уникального (1D) идентификатора. Мы реализуем это, представляя каждую ось в виде значения с фиксированной запятой, сохраняя бит со знаком в самом старшем бите. Затем мы передаём эти три беззнаковые целочисленные значения в функцию генерации хэша. На основании проведённого Jarzynski et al. [7] исследования мы решили использовать для этого этапа PCG.

uint GetSeedWithSignedBits(in float3 discretePos, in uint fixedPointFractionalBits)
{
    //Convert float value to fixed point representation for hash functions
    //Putting the signs of all 3 axis in the 3 most significant bits
    const uint fracScale = 1 << fixedPointFractionalBits;
    
    const float3 rawFractionalPart = frac(discretePos);
    const float3 rawIntegerPart = trunc(discretePos);
    
    const uint3 fp = uint3(rawFractionalPart * fracScale); //fractional part
    const uint3 ip = (uint3(abs(rawIntegerPart)) << fixedPointFractionalBits) & 0x1FFFFFFF;
    const uint xp = clamp(sign(rawIntegerPart.x) * -1, 0, 1) << 31;
    const uint yp = clamp(sign(rawIntegerPart.y) * -1, 0, 1) << 30;
    const uint zp = clamp(sign(rawIntegerPart.z) * -1, 0, 1) << 29;
    
    const uint3 fic = fp | ip;
    const uint3 inp = uint3(xp | fic.x, yp | fic.y, zp | fic.z);
    return pcg_nested(inp);
}

Получив уникальный идентификатор, мы используем его как порождающее значение (seed value) для генерации 3D-точек. Так как мы работаем с 3D-объёмом, одно seed value используется для генерации восьми косвенных точек, по одной для каждого субвоксела. Для генерации каждой точки мы просто выполняем инкремент базового seed value, а затем вызываем функцию. Это даёт нам паттерн, похожий на заполняющую пространство кривую. Важно, чтобы используемая функция генерации была прогрессивной, так как существуют различные типы генераторов. В исследовании, проведённом Christensen et al. [8], есть хороший обзор разных типов.

//Generate a sample in the voxel within a certain quadrant.
//The sample is a normalized value within the voxel itself.
float3 GenerateSampleInVoxelQuadrant(in uint seed, in uint level, in uint quadrant)
{
    //Generate normalized sample, based on seed and offset based on quadrant (progressive required)
    const float3 sample = hash3(seed + quadrant);
    const float3 offsets = (float3)OffsetVectorFromQuadrant(quadrant);
    const float cellDelta = 0.5f;
    //Return remapped to level
    return (offsets * cellDelta) + (sample / (1.f / cellDelta));
}

Ниже показано тело кода шейдера для этого прохода. В него добавлен один не упомянутый выше аспект: мы не выбираем ближайшую точку, а используем вероятность на основе барицентрических координат для случайного выбора одной из трёх ближайших точек. Это позволяет обеспечивать более плавный переход между дискретными точками, как показано на изображении ниже. Из-за использования методики со случайными точками также важно использовать соседние seed value, представляющие соседние косвенные вокселы, потому что мы не знаем, где могут находиться случайные косвенные точки. Может оказаться, что случайная точка в соседнем косвенном вокселе ближе к позиции в мире, чем случайная точка, сгенерированная для текущего косвенного воксела, в котором мы находимся. Для обеспечения максимальной точности мы итеративно обходим все соседние вокселы, благодаря чему проверяем 27 косвенных вокселов.

uint GetIndexByProbability(in float3 position, in ClosestPointSample closestSamples[3])
{  
    //Create probability values based on "area" (actually just distances)
    const float3 distances = float3(
        dot(position - closestSamples[0].sample, position - closestSamples[0].sample),
        dot(position - closestSamples[1].sample, position - closestSamples[1].sample),
        dot(position - closestSamples[2].sample, position - closestSamples[2].sample));
    float3 areas = float3(
        distances[1] * distances[2],
        distances[0] * distances[2],
        distances[0] * distances[1]);
    const float totalArea = areas.x + areas.y + areas.z;
    areas /= totalArea; //normalize

    //Get random normalized value, seed based on position
    const float randomValue = randomNormalizedFloat(GetSeed(position, 15));
    
    //Get index based on probability
    return saturate(uint(randomValue / areas.x)) + saturate(uint(randomValue / (areas.x + areas.y)));
}

//Main function
ClosestPointSample GetGeneratedSample(in float3 position, in uint discreteCellSize, in uint level, 
    in uint voxelConnectivity, in uint maxLevels)
{  
    const float FLT_MAX = 3.402823466e+38f;  
    ClosestPointSample closestSamples[3];
    [unroll]
    for (int i = 0; i < 3; ++i)
    {
        closestSamples[i].sample = float3(FLT_MAX, FLT_MAX, FLT_MAX);
        closestSamples[i].seed = 0;
    }
    
    //Discretize position based on discrete cell size and current level
    const float deltaOnLevel = 1.f / (1 << level);
    const float cellSizeBasedOnLevel = discreteCellSize * deltaOnLevel;
    const float3 discretePosition = Discretize(position, cellSizeBasedOnLevel);
    
    //For that discrete position, interate over all the neighbouring voxels(based on the size of cell size)
    for (uint v = 0; v < (voxelConnectivity + 1); ++v)
    {
        //Get the offset and calculate the neighbouring absolute position
        const int3 neighbourOffsets = NeighbourOffsets[v];
        const float3 neighbourPosition = discretePosition + (neighbourOffsets * cellSizeBasedOnLevel);
        
        //Get the top level discrete position of the neighbour position. This is due to the fact 
        //that we want to find the quadrant in a normalized fashion and the fact that the neighbour 
        //position might be in the same or another voxel as the position!
        const float3 neighbourDiscreteTopPosition = Discretize(neighbourPosition, discreteCellSize);
        const float3 normalizedRelativePosition = 
            abs((neighbourDiscreteTopPosition - neighbourPosition) / (int) discreteCellSize);
        const float deltaNextLevel = deltaOnLevel * 0.5f;
        const float3 centerNormalizedRelativePosition = 
            normalizedRelativePosition + float3(deltaNextLevel, deltaNextLevel, deltaNextLevel);
        
        //Create samples in the current sub voxel & find closest!
        const float3 acp = neighbourDiscreteTopPosition 
            + (centerNormalizedRelativePosition * discreteCellSize);
        const uint absoluteSeed = GetSeedWithSignedBits(acp, maxLevels);
        
        [unroll]
        for (uint i = 0; i < 8; ++i)
        {
            uint currentSeed = absoluteSeed + i;
            uint currentLevel = level;
        
            const float3 normalizedSample = GenerateSampleInVoxelQuadrant(absoluteSeed, level, i);
            float3 absoluteSample = neighbourPosition + (normalizedSample * cellSizeBasedOnLevel);
          
            GetClosestSampleTriplet(position, absoluteSample, currentSeed, closestSamples);
        }
    }

    const uint index = GetIndexByProbability(position, closestSamples);
    return closestSamples[index];
}

Обновление структур данных

Следующий этап — это «проход накопления». Во время этого прохода мы просто накапливаем всю информацию затенения при помощи атомарной операции сложения. Этот проход можно изменять в зависимости от конкретной реализации хэш-таблиц. Пример, основанный на реализации Дэвида Фаррелла, можно найти в полном исходном коде прототипа.

Более важным является этап объединения накопленных данных и постоянной хэш-таблицы. Так как информация AO хранится в интервале [0, 1], мы можем сохранить информацию затенения в виде с фиксированной запятой, где основная часть битов используется для дробной части. На этапе накопления это невозможно сделать, потому что потенциально мы можем накапливать множество нормализованных значений. Другими словами, запись с фиксированной запятой в обеих структурах данных слегка различается. На этапе слияния мы учитываем это и преобразуем значения в соответствующий интервал. Затем мы также изменяем на константу масштаб готового значения, которое сохраняем в постоянную хэш-таблицу. Это нужно для предотвращения переполнения данных на случай, если в один элемент записывается множество данных AO.

[numthreads(8, 8, 1)]
void main(uint3 DTid : SV_DispatchThreadID)
{
    const uint2 screenDimensions = uint2(1920, 1080);
    
    //Input Data
    const uint index = uint(DTid.x + (DTid.y * screenDimensions.x));
    if (index >= AccumulationHashTableElementCount)
        return;
    
    //Store in World Hash Table
    const KeyData accumulatedData = HashTableLookupBySlotID(AccumulationBuffer, 
        AccumulationHashTableElementCount, index);
    if (accumulatedData.key == 0) //No data accumulated, exit
        return;
    
    //Because no additional hashing in hashtable, the key is the seed value used
    const uint se = accumulatedData.key;
    const KeyData cachedData = HashTableLookup(WorldHashTable, WorldHashTableElementCount, se);
    
    //Perform constant rescale to prevent overflow
    const float constantOverflowMultiplier = 0.9f;
    const float scaledAccumulatedCount = accumulatedData.count * constantOverflowMultiplier;
    const float cachedCount = FromFixedPoint(cachedData.count, WorldHashTableCountFractionalBits);
    
    const float totalCount = scaledAccumulatedCount + cachedCount;
    const float cachedValueTerm = (cachedCount / totalCount) 
        * FromFixedPoint(cachedData.value, WorldHashTableValueFractionalBits);
    const float accumulationValueTerm = (scaledAccumulatedCount / totalCount) 
        * (FromFixedPoint(accumulatedData.value, AccumulationHashTableValueFractionalBits) 
            / accumulatedData.count);
    
    const uint valueToStore = ToFixedPoint(cachedValueTerm + accumulationValueTerm, 
        WorldHashTableValueFractionalBits);
    const uint countToStore = ToFixedPoint(totalCount, WorldHashTableCountFractionalBits);

    //Store in World Hash Table
    HashTableInsert(WorldHashTable, WorldHashTableElementCount, se, valueToStore, countToStore);
}

Визуализация

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

//Visualize the pixel by interpolating in world space using a modified Shepard interpolation.
//Single level as it doesn't guarantee data for previous levels.
float WorldSpaceShepardInterpolationSingleLevel(in uint discreteCellSize, in uint level,
    in float3 pixelWorldPosition, in uint maxLevels)
{
    //Shepard Interpolation Variables
    const float searchRadius = float(discreteCellSize) / float(1 << level);
    float valueSum = 0.f;
    float weightSum = 0.f;
    
    //For that discrete position, interate over all the neighbouring voxels
    //(based on the size of cell size)
    for (uint v = 0; v < (VoxelConnectivity + 1); ++v)
    {
        //Seed Generation similar to Implicit Point Generation
		//...
        for (uint i = 0; i < 8; ++i)
        {
			//...
            //Interpolation - Due to going over all points per subvoxel, duplicates are not possible
            const KeyData cachedData = 
                HashTableLookup(WorldHashTable, WorldHashTableElementCount, currentSample.seed);
            if (cachedData.key != 0)
            {
                const float sampleCount = FromFixedPoint(cachedData.count, 
                    WorldHashTableCountFractionalBits);
                const float dist = distance(pixelWorldPosition, currentSample.sample);
                const float weight = clamp(searchRadius - dist, 0.f, searchRadius) * sampleCount;
                
                valueSum += FromFixedPoint(cachedData.value, WorldHashTableValueFractionalBits) * weight;
                weightSum += weight;
            }
        }
    }
    
    if (weightSum > 0)
        valueSum /= weightSum;
    
    return valueSum;
}

Результаты и дальнейшее развитие

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


Хотя ещё многое можно усовершенствовать (например, использовать адаптивную плотность точек на основе уровня детализации или сложности локального затенения, переменную частоту обновления, поддержку динамических изменений при помощи временного интегрирования на основе экспоненциальных средних значений перемещений, и т.п.), первые результаты выглядят интересно. В нашем исследовании мы сравнили результаты с частичной реализацией техники path space filtering компании Nvidia. Метрики, результаты, наши выводы и подробности можно найти в дипломной работе.