Участник:Vlamakarenko/Трассировка лучей

авторы: В.А.Макаренко, Р.А.Габдуллин

1 Свойства и структура алгоритма

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

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

1.2 Математическое описание алгоритма

1.2.1 Некоторые понятия из радиометрии

Энергетическая яркость

BRDF

Энергия излучения - энергия, переносимая оптическим излучением. Обозначается через  [math]Q[/math].

Поток излучения - мощность, переносимая оптическим излучением через какую-либо поверхность.  [math]\Phi = \frac{dQ}{dt}[/math].

Плотность потока излучения - отношение потока излучения, проходящего через малый участок поверхности, к площади этой поверхности.   [math]\frac{d\Phi}{dA}[/math],   где   [math]dA[/math] - дифференциал поверхности.

Облученность - отношение потока излучения, падающего на малый участок поверхности, к площади этой поверхности.   [math]E=\frac{d\Phi}{dA}[/math],   где   [math]dA[/math] - дифференциал поверхности.

Энергетическая светимость - отношение потока излучения, испускаемого малым участком поверхности, к площади этой поверхности.  [math]M=\frac{d\Phi}{dA}[/math],   где   [math]dA[/math] - дифференциал поверхности.

Сила излучения - мощность излучения источника света в некотором направлении. Равна отношению потока излучения, исходящего от источника, к малому телесному углу.   [math]I=\frac{d\Phi}{d\omega}[/math].

Энергетическая яркость - отношение потока излучения через малую площадку поверхности и распространяющегося в малом телесном угле, к площади проекции этой площадки на плоскость, перпендикулярную направлению распространения, и величине телесного угла.  [math]L=\frac{d^2\Phi}{dA\cdot d\omega \cdot \cos \theta}[/math]

Энергетическую яркость, падающую на элемент поверхности, будем назвать падающей яркостью и обозначать [math]L_i[/math]. Энергетическую яркость, отраженную элементом поверхности, будем называть отраженной яркостью и обозначать [math]L_o[/math].

Из определений следует, что:

1) [math]dE_i(\overline{p}, \overline{\omega}_i) = L_i(\overline{p}, \overline{\omega}_i)\cos \theta_i \, d\omega_i[/math]
2) [math]E(\overline{p}) = \int_{\Omega_i} L_i(\overline{p}, \overline{\omega}_i)\cos \theta_i \, d\omega_i[/math],  где интегрирование проводится по направлениям [math]\overline{\omega}_i[/math] таким, что [math](\overline{n}, \overline{\omega}_i)[/math] > 0,  то есть по полусфере.
([math]\overline{n}[/math] - нормаль к поверхности в точке [math]\overline{p}[/math],   [math]\Omega_i = 2\pi^{+}[/math]).

Пусть [math]\overline{p}[/math] - точка поверхности, [math]\overline{n}[/math] - нормаль в этой точке, [math]\overline{\omega}_i, \overline{\omega}_o[/math] - единичные вектора такие, что [math](\overline{n}, \overline{\omega}_i) \gt 0, \,\, (\overline{n}, \overline{\omega}_o) \gt 0 [/math]
Будем считать, что   [math]dL_o(\overline{p}, \overline{\omega}_o)[/math] пропорционален   [math]dE_i(\overline{p},\overline{\omega}_i)[/math]:
[math]dL_o(\overline{p}, \overline{\omega}_o) \sim dE_i(\overline{p},\overline{\omega}_i)[/math].
Функция   [math]f_r(\overline{p}, \overline{\omega}_i, \overline{\omega}_o) = \frac{dL_o(\overline{p}, \overline{\omega}_o)}{dE_i(\overline{p},\overline{\omega}_i)} = \frac{dL_o(\overline{p}, \overline{\omega}_o)}{L_i(\overline{p}, \overline{\omega}_i)\cos \theta_i d\omega_i}[/math]   называется двулучевой функцией отражательной способности (BRDF).
В точке [math]\overline{p}[/math] может быть задано несколько BRDF - множество [math]\{f_{r_k}(\overline{p}, \overline{\omega}_i, \overline{\omega}_o)\}[/math].
Тогда итоговая BRDF просто равна сумме BRDF из этого множества: [math]f_r(\overline{p}, \overline{\omega}_i, \overline{\omega}_o) = \sum_k f_{r_k}(\overline{p}, \overline{\omega}_i, \overline{\omega}_o)[/math].
[math]L_o(\overline{p}, \overline{\omega}_o) = \int_{\Omega_i} f_r(\overline{p}, \overline{\omega}_i, \overline{\omega}_o) L_i(\overline{p}, \overline{\omega}_i)\cos \theta_i d\omega_i[/math].   - уравнение отражений.
Важная часть алгоритма трассировки лучей - определить отраженную яркость по падающим яркостям, то есть посчитать этот интеграл. Обычно для вычисления интеграла применяют метод Монте-Карло.

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

Уравнение рендеринга - это интегральное уравнение, которое определяет количество излученного в определенном направлении света как сумму отраженного и собственного излучения. Физической основой для уравнения рендеринга является закон сохранения энергии.

Уравнение отражений является частью уравнения рендеринга, поэтому остается только добавить собственное излучение объектов [math]L_e(\overline{p}, \overline{\omega}_o)[/math].

Уравнение рендеринга имеет вид [math]L_o(\overline{p}, \overline{\omega}_o) = L_e(\overline{p}, \overline{\omega}_o) + \int_{\Omega_i} f_r(\overline{p}, \overline{\omega}_i, \overline{\omega}_o) L_i(\overline{p}, \overline{\omega}_i)\cos \theta_i d\omega_i[/math].

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

1.2.1.2 Поиск приближенного решения уравнения рендеринга

Введем оператор поиска первого пересечения. [math]r_c(\overline{p}, \overline{\omega})[/math] - оператор, который находит точку первого пересечения с объектами луча, выпускаемого из точки [math]\overline{p}[/math] по направлению [math]\overline{\omega}[/math].

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

[math] L_i(\overline{p},\overline{\omega}_i) = L_o(r_c(\overline{p}, \overline{\omega}_i), -\overline{\omega}_i) [/math].

и уравнение рендеринга превратится в уравнение Фредгольма второго рода

[math] L_o(\overline{p},\overline{\omega}_o) = L_e(\overline{p},\overline{\omega}_o) + \int_{\Omega_i} f_{r}(\overline{p},\overline{\omega}_i,\overline{\omega}_o)L_o(r_c(\overline{p}, \overline{\omega}_i), -\overline{\omega}_i) \cos\theta_i d\omega_i [/math].

Обозначим интеграл в правой части как [math] K \circ L_o [/math]. Тогда уравнение рендеринга перепишется, как

[math] L_o = L_e + K \circ L_o [/math].

Искать решение будем итерационным методом: выберем какое-либо начальное приближение [math] L_{o}^{(0)} [/math]. Найдем значение правой части [math]L_o^{(1)} = L_e + K \circ L_o^{(0)}[/math], воспользуемся им для поиска следующего приближенного решения [math]L_o^{(2)} = L_e + K \circ L_o^{(1)} = L_e + K \circ ( L_e + L_o^{(0)} ) = L_e + K \circ L_e + K^2 \circ L_o^{(0)} [/math] и т.д.

В конечном итоге получим, что [math]L_o = \sum \limits_{n=0}^\infty K^n \circ L_o^{(0)}[/math].

Этот метод называется методом последовательных приближений.

1.2.2 Определения

Камера и видовая плоскость

Объект - поверхность в трехмерном пространстве.
Материал объекта - множество BRDF.
Сцена - четверка {камера, объекты, источники света, видовая плоскость}.

1.2.2.1 Камера

Камера - "виртуальный глаз", через который мы видим объекты сцены.

Камера задается положением в пространстве и ортонормированным базисом [math](\overline{u}_c, \overline{v}_c, \overline{w}_c)[/math] (локальная система координат). [math]\overline{u}_c[/math] направлен вправо относительно наблюдателя, [math]\overline{v}_c[/math] - вверх, [math]\overline{w}_c[/math] - на наблюдателя. Таким образом, [math](\overline{u}_c, \overline{v}_c, \overline{w}_c)[/math] - правая тройка векторов. Вектор направления камеры в своей системе координат: [math](0, 0, -1)[/math].

1.2.2.2 Видовая плоскость

Перед камерой на расстоянии [math]d[/math] поставим пиксельную сетку [math]w \times h[/math] с размером пикселя [math]s \times s[/math] так, чтобы камера была направлена в ее центр. Итоговое изображение будет получено из этой сетки покраской пикселей в нужный цвет.

[math]p(i, j)[/math], где [math]i=\overline{0,h-1}, \quad j=\overline{0,w-1}[/math] - пиксел [math]i[/math]-ой строки снизу [math]j[/math]-го столбца слева.

Заметим, что с увеличением [math]d[/math] уменьшается угол обзора, а [math]s[/math] влияет на масштаб.

1.2.2.3 Источники света

Типы источников света:
1) Point light - точка, излучающая свет одиноково во всех направлениях.
2) Directional light - бесконечно удаленный источник света, задается единичным направлением на себя.
3) Area light - объект, излучающий свет одинаково во всех направлениях.

1.2.3 Алгоритм трассировки лучей

1.2.3.1 Генерация луча

Перед камерой на расстоянии [math]d[/math] поставим пиксельную сетку [math]w \times h[/math] с размером пикселя [math]s \times s[/math] так, чтобы камера была направлена в центр сетки. Итоговое изображение будет получено из этой сетки покраской пикселей в нужный цвет.

Пусть [math]p(i, j)[/math], где [math]i=\overline{0,h-1}, \quad j=\overline{0,w-1}[/math] - пиксел [math]i[/math]-ой строки снизу [math]j[/math]-го столбца слева. Проведем пучок лучей из камеры через пиксел в локальной системе координат камеры. Начало каждого луча в точке [math]\overline{pos}=(0,0,0)[/math]. Найдем вектора направления [math]\overline{dir}_n[/math]:

[math]\overline{dir}_n.x = s \cdot (j + \delta_n.x -\frac{w}{2} + \frac{1}{2}) [/math]

[math]\overline{dir}_n.y = s \cdot (i + \delta_n.y -\frac{h}{2} + \frac{1}{2}) [/math]

[math]\overline{dir}_n.z = d[/math],

где [math] \delta_n = (\delta_n.x, \delta_n.y) \in [0;1]^2[/math] - смещения начала n-ого луча.

Нормируем вектора [math]\overline{dir}_n[/math] и переведем лучи в глобальную систему координат.

1.2.3.2 Пересечение луча с объектами сцены

Все точки луча [math]ray(\overline{p}, \overline{d})[/math] можно представить в виде:

[math]h(x,y,z)=\overline{p}+t \cdot \overline{d}, \quad t \in [0;\infty) [/math].

Луч пересекает объект [math]O_i[/math] тогда и только тогда, когда [math]t_i = \inf\{t \gt 0 \, | \quad \overline{p}+t \cdot \overline{d} \in O_i\} \, \lt \infty[/math]

Если луч пересекает объект [math]O_i[/math], то [math]h_i(x,y,z) = \overline{p}+t_i \cdot \overline{d}[/math] - их первое пересечение, [math]t_i[/math] - расстояние от начала луча до точки пересечения.

1.2.3.3 Вычисление энергитической яркости [math]L_o[/math], выходящей из точки пересечения

Найдем [math]T = \{t_i\}[/math] - параметры пересечения с объектами (см. предыдущий пункт).

[math]t_{min} = \min \{t_i\}, \quad O_{min} = \{O_i | t_i = t_{min}\}[/math]

Если [math]t_{min} = \infty[/math], то луч не пересекает ни один объект и нужно покрасить соответствующий пиксел в черный цвет (цвет фона).

Если [math]t_{min} \lt \infty[/math], то нужно найти выходящее излучение (цвет) из точки пересечения в камеру. Выходящее излучение зависит от свойств материала объекта и входящего излучения в точку пересечения.

Чтобы вычислить входящее излучение, находится приближенное значение интеграла [math] \int_{\Omega_i} f_{r}(\overline{p},\overline{\omega}_i,\overline{\omega}_o)L_o(r_c(\overline{p}, \overline{\omega}_i), -\overline{\omega}_i) \cos\theta_i d\omega_i [/math] методом Монте-Карло.

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

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

1.2.4 Улучшение визуальных качеств изображения

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

1.3 Вычислительное ядро алгоритма

Основные вычисления связаны с поиском пересечения луча с объектами сцены.

1.4 Макроструктура алгоритма

Для каждого пикселя изображения генерируется пучок лучей, проходящих через этот пиксел.

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

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

1.5 Схема реализации последовательного алгоритма

Здесь приведен псевдокод ключевых моментов алгоритма.

//Покраска пикселей изображения

    for (uint r = 0; r < h; r++) //строки
    {
        for (uint c = 0; c < w; c++) //столбцы
        {
            image(r,c) = renderPixel(r, c); //вычисляем цвет в пикселе
        }
    }

    //renderPixel(r, c)

    pixel_color = {0,0,0};

    for (uint j = 0; j < num_rays_per_pixel; j++)
    {
        pp = generatePoint(); //генерируем произвольное смещение внутри пикеля

        //вычисляем координаты точки на полотне
        x = s * (c - 0.5 * w + pp.x);
        y = s * (r - 0.5 * h + pp.y);

        //cоздаем луч
        Ray ray({0,0,0}, normalize({x,y,-d}));
        Ray worldRay = convertToWorldCoords(ray);
        
        pixel_color += traceRay(worldRay, 0) * 1.0 / num_rays_per_pixel; //вычисляем цвет для каждого луча и усредняем
    }

    //traceRay(ray, depth)

    if (depth > max_depth) //превысили максимальную глубину луча, дальше не идем
    {
        return black; //возвращаем черный цвет
    }

    Intersection is = hitObjects(ray); //нашли первое пересечение с объектами сцены и записали информацию о нем в is:
                                       //материал объекта, расстояние, нормаль и т. п.

    if (is) //если есть пересечение
    {
        is.ray = ray; //запомнили луч
        is.depth = depth; // и глубину

        return shade(is); //расчитали цвет в точке пересечения в зависимости
                          //от материала объекта, используя информацию о пересечении
                                        
    }

    return getBackgroundColor(ray); //вернули цвет фона

//shade(is)

//цвета
direct_illumination = {0,0,0}; //прямое излучение
indirect_illumintaion = {0,0,0}; //непрямое излучение

for (uint i = 0; i < light_count; i++)
{
    if (light[i].isVisible(is)) //если источник света виден из точки пересечения
    {
         direct_illumination = getIllumination(light[i], is); //рассчитываем прямое излучение от i-го источника света
    }
}


if (isReflective(is.object))
{
      indirect_illumination = traceRay(reflecedRay(is.ray), is.depth + 1);
}

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

1.6 Последовательная сложность алгоритма

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

[math]w[/math] - ширина полотна

[math]h[/math] - высота полотна

[math]depth_{max}[/math] - максимальная глубина луча в рекурсии

[math]N_{obj}[/math] - количество объектов сцены

[math]N_{lights}[/math] - количество источников света сцены

[math]N_{rays}[/math] - число лучей на один пиксел

На каждой итерации ищется пересечение луча со всеми объектами: [math]O(N_{obj})[/math]

После нахождения точки пересечения, выясняется, виден ли каждый источник света из этой точки. Для этого нужно найти первое пересечение теневого луча (начало в точке пересечения, направлен к источнику света) и сравнить расстояние до точки пересечения с расстоянием до источника света. Для всего этого требуется [math]O(N_{lights} \cdot N_{obj})[/math] операций.

Для каждого луча понадобится [math]O(depth_{max}(N_{obj} + N_{obj} \cdot N_{lights})) = O(depth_{max}\cdot N_{obj} \cdot N_{lights})[/math] операций.

Всего лучей, выходящих из камеры [math]w \cdot h \cdot N_{rays}[/math].

В итоге получаем [math]O(depth_{max}\cdot N_{obj} \cdot N_{lights} \cdot w \cdot h \cdot N_{rays})[/math] операций.

1.7 Информационный граф

Зеленые вершины - генерация лучей для данного пикселя.
Желтые вершины - трассировка лучей данного пикселя.

1.8 Ресурс параллелизма алгоритма

Из информационного графа видно, что обработка каждого пикселя происходит независимо от других, а затем происходит соединение всех пикселей в итоговое изображение. Пусть [math]N_p \geqslant 2[/math] - количество процессов. Всем процессам, кроме первого поручим обработать свою часть изображения (каждый процесс обработает [math]\Big[\frac{w \cdot h}{N_p - 1}\Big] [/math] пикселей, а оставшуюся часть работы поручим произвольному процессу). Первому процессу поручим принять работу и собрать итоговое изображение (используем MPI_Gatherv).

Время работы: [math]O\Big(\frac{depth_{max}\cdot N_{obj} \cdot N_{lights} \cdot w \cdot h \cdot N_{rays}}{N_p}\Big)[/math].

1.9 Входные и выходные данные алгоритма

На вход алгоритму подается описание сцены (описание камеры, видовой плоскости, объектов и их материалов, источников света).
Выходные данные алгоритма - готовое изображение.

2 Программная реализация алгоритма

2.4 Масштабируемость алгоритма и его реализации

2.4.1 Реализация с помощью MPI

Исследование проводилось на суперкомпьютере "Ломоносов".

Компилятор C++
Intel 13.0.1 + impi 4.1.0

1)
Размер изображения 700 на 700 пикселей, 256 лучей на пиксел, глубина рекурсии 3.
Сцена:

Число процессов	Время (с)
1	479.915
2	492.614
4	313.172
8	170.866
16	82.375
32	43.812
64	22.318
128	12.195
256	8.793
512	6.954
1024	6.741

2)
Размер изображения 1340 на 1340 пикселей, 1024 лучей на пиксел, глубина рекурсии 5.
Сцена:

Число процессов	Время (с)
1	600.270
2	610.583
4	243.261
8	136.025
16	67.517
32	33.467
64	17.234
128	9.361

Алгоритм перестает масштабироваться, когда время работы параллельной части алгоритма (визуализация) мало по сравнению с накладными расходами MPI и временем, которое расходуется на инициализацию сцены (загрузка сцены из файла, построение оптимизирующих структур данных для ускорения визуализации 3d-моделей).

2.4.2 Реализация с помощью NVIDIA CUDA

Использовалось следующее программное обеспечение:

OS
Microsoft Windows 10 Education Edition
Ubuntu 16.10 Yakkety Yak

Компиляторы C++
nvcc из пакета NVIDIA CUDA SDK 8.0.44 + Microsoft Visual C++ 2015(Windows)
nvcc из пакета NVIDIA CUDA SDK 8.0.44 + gcc-4.9(Linux)

Результаты тестирования получены на следующих машинах:

CPU	RAM	GPU	OS
Core i5-3210m 3.1GHz	6Gb DDR3 1333MHz	NVIDIA GT 640m 2.0Gb DDR3	Ubuntu 16.10
Core i5-6600 3.9GHz	16Gb DDR4 2133MHz	NVIDIA GTX 1070 8.0Gb GDDR5	Ubuntu 16.10

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

Лучей в окне	GTX 1070	GT 640m
512	29.8 c	58.0 c
1024	15.2 c	35.5 c
2048	8.2 c	25.2 c
4096	4.6 c	24.5 c
8192	3.0 c	23.5 c
16384	3.1 c	23.4 c
32768	2.8 c	23.2 c
65536	2.5 c	23.1 c
131072	2.4 c	23.0 c

Видно, что с ростом размера окна "исчезает" рост производительности. Это связано с накладными расходами на хранение данных в памяти GPU.

2.7 Существующие реализации алгоритма

RenderMan — очень известный пакет программ для рендеринга 3D-анимации, который был создан компанией Pixar в 1986 году. Данный пакет использовался во множестве известных фильмов и мультфильмов, например: Аватар(2009), Гарри Поттер и Орден Феникса(2007), В поисках Немо(2003).

V-Ray — система рендеринга, разработанная компанией Chaos Group в 2000 году. VRay — это рейтрейсный рендерер, в котором присутствует несколько алгоритмов просчета глобального освещения.

NVIDIA IRay — технология создания фотореалистичных изображений с использованием графических ускорителей, созданная компанией NVIDIA.

mental ray — профессиональная система рендеринга и визуализации изображений, разработанная компанией mental images(является дочерней компанией NVIDIA) в 1986 году. Использован в фильмах Терминатор 3: Восстание машин(2003), Трон: Наследие(2010).

Hydra Renderer — это плагин для рендера реалистичных изображений в 3D Studio Max. Разработан при поддержке Лаборатории компьютерной графики и мультимедиа факультета ВМК МГУ.

3 Литература

1. Suffern, Kevin G. "Ray Tracing from the Ground Up / Kevin G. Suffern", 2007

2. Scratchapixel 2.0. Learn Computer Graphics From Scratch!

3. Параллельные вычисления на GPU. Архитектура и программная модель CUDA: Учеб. пособие / А.В. Боресков и др. - М.: Издательство Московского университета, 2012. - 336 с.

4 Полученные изображения

• 
• 
• 
• 
• 
• 
• 
•