Участник:Vlamakarenko/Трассировка лучей
авторы: В.А.Макаренко, Р.А.Габдуллин
Свойства и структура алгоритма
Общее описание алгоритма
Трассировка лучей - это технология визуализации трехмерных сцен путем отслеживания обратной траектории лучей света (от наблюдателя к источнику). Достоинствами данного метода являются реалистичность итоговых изображений, возможность визуализации гладких объектов без аппроксимации полигональными поверхностями, простота реализации отражений, преломлений, взятия в фокус, реалистичного освещения; возможность параллельной обработки лучей.
Математическое описание алгоритма
Некоторые понятия из радиометрии

Энергия излучения - энергия, переносимая оптическим излучением. Обозначается через [math]\displaystyle{ Q }[/math].
Поток излучения - мощность, переносимая оптическим излучением через какую-либо поверхность. [math]\displaystyle{ \Phi = \frac{dQ}{dt} }[/math].
Плотность потока излучения - отношение потока излучения, проходящего через малый участок поверхности, к площади этой поверхности. [math]\displaystyle{ \frac{d\Phi}{dA} }[/math], где [math]\displaystyle{ dA }[/math] - дифференциал поверхности.
Облученность - отношение потока излучения, падающего на малый участок поверхности, к площади этой поверхности. [math]\displaystyle{ E=\frac{d\Phi}{dA} }[/math], где [math]\displaystyle{ dA }[/math] - дифференциал поверхности.
Энергетическая светимость - отношение потока излучения, испускаемого малым участком поверхности, к площади этой поверхности. [math]\displaystyle{ M=\frac{d\Phi}{dA} }[/math], где [math]\displaystyle{ dA }[/math] - дифференциал поверхности.
Сила излучения - мощность излучения источника света в некотором направлении. Равна отношению потока излучения, исходящего от источника, к малому телесному углу. [math]\displaystyle{ I=\frac{d\Phi}{d\omega} }[/math].
Энергетическая яркость - отношение потока излучения через малую площадку поверхности и распространяющегося в малом телесном угле, к площади проекции этой площадки на плоскость, перпендикулярную направлению распространения, и величине телесного угла. [math]\displaystyle{ L=\frac{d^2\Phi}{dA\cdot d\omega \cdot \cos \theta} }[/math]
Энергетическую яркость, падающую на элемент поверхности, будем назвать падающей яркостью и обозначать [math]\displaystyle{ L_i }[/math]. Энергетическую яркость, отраженную элементом поверхности, будем называть отраженной яркостью и обозначать [math]\displaystyle{ L_o }[/math].
Из определений следует, что:
1) [math]\displaystyle{ dE_i(\overline{p}, \overline{\omega}_i) = L_i(\overline{p}, \overline{\omega}_i)\cos \theta_i \, d\omega_i }[/math]
2) [math]\displaystyle{ E(\overline{p}) = \int_{\Omega_i} L_i(\overline{p}, \overline{\omega}_i)\cos \theta_i \, d\omega_i }[/math], где интегрирование проводится по направлениям [math]\displaystyle{ \overline{\omega}_i }[/math] таким, что [math]\displaystyle{ (\overline{n}, \overline{\omega}_i) }[/math] > 0, то есть по полусфере.
([math]\displaystyle{ \overline{n} }[/math] - нормаль к поверхности в точке [math]\displaystyle{ \overline{p} }[/math], [math]\displaystyle{ \Omega_i = 2\pi^{+} }[/math]).
Пусть [math]\displaystyle{ \overline{p} }[/math] - точка поверхности, [math]\displaystyle{ \overline{n} }[/math] - нормаль в этой точке, [math]\displaystyle{ \overline{\omega}_i, \overline{\omega}_o }[/math] - единичные вектора такие, что [math]\displaystyle{ (\overline{n}, \overline{\omega}_i) \gt 0, \,\, (\overline{n}, \overline{\omega}_o) \gt 0 }[/math]
Будем считать, что [math]\displaystyle{ dL_o(\overline{p}, \overline{\omega}_o) }[/math] пропорционален [math]\displaystyle{ dE_i(\overline{p},\overline{\omega}_i) }[/math]:
[math]\displaystyle{ dL_o(\overline{p}, \overline{\omega}_o) \sim dE_i(\overline{p},\overline{\omega}_i) }[/math].
Функция [math]\displaystyle{ 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]\displaystyle{ \overline{p} }[/math] может быть задано несколько BRDF - множество [math]\displaystyle{ \{f_{r_k}(\overline{p}, \overline{\omega}_i, \overline{\omega}_o)\} }[/math].
Тогда итоговая BRDF просто равна сумме BRDF из этого множества: [math]\displaystyle{ 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]\displaystyle{ 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]. - уравнение отражений.
Важная часть алгоритма трассировки лучей - определить отраженную яркость по падающим яркостям, то есть посчитать этот интеграл. Обычно для вычисления интеграла применяют метод Монте-Карло.
Уравнение рендеринга
Уравнение рендеринга - это интегральное уравнение, которое определяет количество излученного в определенном направлении света как сумму отраженного и собственного излучения. Физической основой для уравнения рендеринга является закон сохранения энергии.
Уравнение отражений является частью уравнения рендеринга, поэтому остается только добавить собственное излучение объектов [math]\displaystyle{ L_e(\overline{p}, \overline{\omega}_o) }[/math].
Уравнение рендеринга имеет вид [math]\displaystyle{ 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].
По сути, алгоритмы рендеринга, в том числе и трассировка лучей, занимаются поиском приближенного решения уравнения рендеринга.
Поиск приближенного решения уравнения рендеринга
Введем оператор поиска первого пересечения. [math]\displaystyle{ r_c(\overline{p}, \overline{\omega}) }[/math] - оператор, который находит точку первого пересечения с объектами луча, выпускаемого из точки [math]\displaystyle{ \overline{p} }[/math] по направлению [math]\displaystyle{ \overline{\omega} }[/math].
Воспользуемся предположением, что количество излучения, распространяющегося вдоль одной прямой одинаково в любой точке прямой. Тогда:
[math]\displaystyle{ L_i(\overline{p},\overline{\omega}_i) = L_o(r_c(\overline{p}, \overline{\omega}_i), -\overline{\omega}_i) }[/math].
и уравнение рендеринга превратится в уравнение Фредгольма второго рода
[math]\displaystyle{ 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]\displaystyle{ K \circ L_o }[/math]. Тогда уравнение рендеринга перепишется, как
[math]\displaystyle{ L_o = L_e + K \circ L_o }[/math].
Искать решение будем итерационным методом: выберем какое-либо начальное приближение [math]\displaystyle{ L_{o}^{(0)} }[/math]. Найдем значение правой части [math]\displaystyle{ L_o^{(1)} = L_e + K \circ L_o^{(0)} }[/math], воспользуемся им для поиска следующего приближенного решения [math]\displaystyle{ 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]\displaystyle{ L_o = \sum \limits_{n=0}^\infty K^n \circ L_o^{(0)} }[/math].
Этот метод называется методом последовательных приближений.
Определения

Объект - поверхность в трехмерном пространстве.
Материал объекта - множество BRDF.
Сцена - четверка {камера, объекты, источники света, видовая плоскость}.
Камера
Камера - "виртуальный глаз", через который мы видим объекты сцены.
Камера задается положением в пространстве и ортонормированным базисом [math]\displaystyle{ (\overline{u}_c, \overline{v}_c, \overline{w}_c) }[/math] (локальная система координат). [math]\displaystyle{ \overline{u}_c }[/math] направлен вправо относительно наблюдателя, [math]\displaystyle{ \overline{v}_c }[/math] - вверх, [math]\displaystyle{ \overline{w}_c }[/math] - на наблюдателя. Таким образом, [math]\displaystyle{ (\overline{u}_c, \overline{v}_c, \overline{w}_c) }[/math] - правая тройка векторов. Вектор направления камеры в своей системе координат: [math]\displaystyle{ (0, 0, -1) }[/math].
Видовая плоскость
Перед камерой на расстоянии [math]\displaystyle{ d }[/math] поставим пиксельную сетку [math]\displaystyle{ w \times h }[/math] с размером пикселя [math]\displaystyle{ s \times s }[/math] так, чтобы камера была направлена в ее центр. Итоговое изображение будет получено из этой сетки покраской пикселей в нужный цвет.
[math]\displaystyle{ p(i, j) }[/math], где [math]\displaystyle{ i=\overline{0,h-1}, \quad j=\overline{0,w-1} }[/math] - пиксел [math]\displaystyle{ i }[/math]-ой строки снизу [math]\displaystyle{ j }[/math]-го столбца слева.
Заметим, что с увеличением [math]\displaystyle{ d }[/math] уменьшается угол обзора, а [math]\displaystyle{ s }[/math] влияет на масштаб.
Источники света
Типы источников света:
1) Point light - точка, излучающая свет одиноково во всех направлениях.
2) Directional light - бесконечно удаленный источник света, задается единичным направлением на себя.
3) Area light - объект, излучающий свет одинаково во всех направлениях.
Алгоритм трассировки лучей
Генерация луча
Перед камерой на расстоянии [math]\displaystyle{ d }[/math] поставим пиксельную сетку [math]\displaystyle{ w \times h }[/math] с размером пикселя [math]\displaystyle{ s \times s }[/math] так, чтобы камера была направлена в центр сетки. Итоговое изображение будет получено из этой сетки покраской пикселей в нужный цвет.
Пусть [math]\displaystyle{ p(i, j) }[/math], где [math]\displaystyle{ i=\overline{0,h-1}, \quad j=\overline{0,w-1} }[/math] - пиксел [math]\displaystyle{ i }[/math]-ой строки снизу [math]\displaystyle{ j }[/math]-го столбца слева. Проведем пучок лучей из камеры через пиксел в локальной системе координат камеры. Начало каждого луча в точке [math]\displaystyle{ \overline{pos}=(0,0,0) }[/math]. Найдем вектора направления [math]\displaystyle{ \overline{dir}_n }[/math]:
[math]\displaystyle{ \overline{dir}_n.x = s \cdot (j + \delta_n.x -\frac{w}{2} + \frac{1}{2}) }[/math]
[math]\displaystyle{ \overline{dir}_n.y = s \cdot (i + \delta_n.y -\frac{h}{2} + \frac{1}{2}) }[/math]
[math]\displaystyle{ \overline{dir}_n.z = d }[/math],
где [math]\displaystyle{ \delta_n = (\delta_n.x, \delta_n.y) \in [0;1]^2 }[/math] - смещения начала n-ого луча.
Нормируем вектора [math]\displaystyle{ \overline{dir}_n }[/math] и переведем лучи в глобальную систему координат.
Пересечение луча с объектами сцены
Все точки луча [math]\displaystyle{ ray(\overline{p}, \overline{d}) }[/math] можно представить в виде:
[math]\displaystyle{ h(x,y,z)=\overline{p}+t \cdot \overline{d}, \quad t \in [0;\infty) }[/math].
Луч пересекает объект [math]\displaystyle{ O_i }[/math] тогда и только тогда, когда [math]\displaystyle{ t_i = \inf\{t \gt 0 \, | \quad \overline{p}+t \cdot \overline{d} \in O_i\} \, \lt \infty }[/math]
Если луч пересекает объект [math]\displaystyle{ O_i }[/math], то [math]\displaystyle{ h_i(x,y,z) = \overline{p}+t_i \cdot \overline{d} }[/math] - их первое пересечение, [math]\displaystyle{ t_i }[/math] - расстояние от начала луча до точки пересечения.
Вычисление энергитической яркости [math]\displaystyle{ L_o }[/math], выходящей из точки пересечения
Найдем [math]\displaystyle{ T = \{t_i\} }[/math] - параметры пересечения с объектами (см. предыдущий пункт).
[math]\displaystyle{ t_{min} = \min \{t_i\}, \quad O_{min} = \{O_i | t_i = t_{min}\} }[/math]
Если [math]\displaystyle{ t_{min} = \infty }[/math], то луч не пересекает ни один объект и нужно покрасить соответствующий пиксел в черный цвет (цвет фона).
Если [math]\displaystyle{ t_{min} \lt \infty }[/math], то нужно найти выходящее излучение (цвет) из точки пересечения в камеру. Выходящее излучение зависит от свойств материала объекта и входящего излучения в точку пересечения.
Чтобы вычислить входящее излучение, находится приближенное значение интеграла [math]\displaystyle{ \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] методом Монте-Карло.
Заметим, что каждый луч необходимо продолжить в некотором направлении(а не создавать новый пучок), так как необходимое количество лучей для численного интегрирования можно заранее сгенерировать в камере.
Таким образом, алгоритм моделирует обратное движение пучка лучей света, который пришел в камеру. Пучок лучей движется из камеры и, сталкиваясь с объектами, распространяется в пространство в случайных направлениях, согласно свойствам материалов объектов.
Улучшение визуальных качеств изображения
Заметим, что чем больше лучей выпускается через каждый пиксел изображения, тем лучше получается картинка.
Вычислительное ядро алгоритма
Основные вычисления связаны с поиском пересечения луча с объектами сцены.
Макроструктура алгоритма
Для каждого пикселя изображения генерируется пучок лучей, проходящих через этот пиксел.
Для каждого луча находится точка первого пересечения с объектами сцены, расчитывается входящие излучения и на основе этого - итоговое выходящее.
После этого цвет усредняется по всем лучам, проходящим через данный пиксел.
Схема реализации последовательного алгоритма
Здесь приведен псевдокод ключевых моментов алгоритма.
//Покраска пикселей изображения
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);
}
//зная входящее излучение, расчитываем выходящее в зависимости от свойств материала объекта
//...
Последовательная сложность алгоритма
[math]\displaystyle{ w }[/math] - ширина полотна
[math]\displaystyle{ h }[/math] - высота полотна
[math]\displaystyle{ depth_{max} }[/math] - максимальная глубина луча в рекурсии (помимо первичного луча)
[math]\displaystyle{ N_{obj} }[/math] - количество объектов сцены
[math]\displaystyle{ N_{lights} }[/math] - количество источников света сцены
[math]\displaystyle{ N_{rays} }[/math] - число лучей на один пиксел
На каждой итерации ищется пересечение луча со всеми объектами: [math]\displaystyle{ O(N_{obj}) }[/math]
После нахождения точки пересечения, выясняется, виден ли каждый источник света из этой точки. Для этого нужно найти первое пересечение теневого луча (начало в точке пересечения, направлен к источнику света) и сравнить расстояние до точки пересечения с расстоянием до источника света. Для всего этого требуется [math]\displaystyle{ O(N_{lights} \cdot N_{obj}) }[/math] операций.
Для каждого луча понадобится [math]\displaystyle{ O(depth_{max}(N_{obj} + N_{obj} \cdot N_{lights})) = O(depth_{max}\cdot N_{obj} \cdot N_{lights}) }[/math] операций.
Всего лучей, выходящих из камеры [math]\displaystyle{ w \cdot h \cdot N_{rays} }[/math].
В итоге получаем [math]\displaystyle{ O(depth_{max}\cdot N_{obj} \cdot N_{lights} \cdot w \cdot h \cdot N_{rays}) }[/math] операций.
Информационный граф
Ресурс параллелизма алгоритма
Входные и выходные данные алгоритма
На вход алгоритму подается описание сцены (описание камеры, видовой плоскости, объектов и их материалов, источников света).
Выходные данные алгоритма - изображение.
Свойства алгоритма
Программная реализация алгоритма
Особенности реализации последовательного алгоритма
Локальность данных и вычислений
Возможные способы и особенности параллельной реализации алгоритма
Масштабируемость алгоритма и его реализации
Динамические характеристики и эффективность реализации алгоритма
Выводы для классов архитектур
Существующие реализации алгоритма
RenderMan
V-Ray
NVIDIA IRay
NVIDIA mental ray
Hydra Renderer
Литература
http://www.raytracegroundup.com/