Участник:Vlamakarenko/Трассировка лучей: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
 
(не показаны 244 промежуточные версии 3 участников)
Строка 1: Строка 1:
авторы: [[Участник:vlamakarenko|В.А.Макаренко]], [[Участник:Ruslixag|Р.А.Габдуллин]]
+
авторы: [[Участник:vlamakarenko|В.А.Макаренко]], [[Участник:Rgabdullin|Р.А.Габдуллин]]
  
 
= Свойства и структура алгоритма =
 
= Свойства и структура алгоритма =
Строка 5: Строка 5:
 
== Общее описание алгоритма ==
 
== Общее описание алгоритма ==
  
Трассировка лучей - это технология визуализации трехмерных сцен путем отслеживания обратной траектории лучей света (от наблюдателя к источнику).
+
'''Трассировка лучей''' - это технология визуализации трехмерных сцен путем отслеживания обратной траектории лучей света (от наблюдателя к источнику).
 
Достоинствами данного метода являются реалистичность итоговых изображений, возможность визуализации гладких объектов без аппроксимации полигональными поверхностями, простота реализации отражений, преломлений, взятия в фокус, реалистичного освещения; возможность параллельной обработки лучей.
 
Достоинствами данного метода являются реалистичность итоговых изображений, возможность визуализации гладких объектов без аппроксимации полигональными поверхностями, простота реализации отражений, преломлений, взятия в фокус, реалистичного освещения; возможность параллельной обработки лучей.
  
 
== Математическое описание алгоритма ==
 
== Математическое описание алгоритма ==
  
=== Базовые определения ===
+
=== Некоторые понятия из радиометрии ===
  
=== Входные данные ===
+
[[file: Radiance_raytracer.png|thumb|right|300px|Энергетическая яркость]]
 +
[[file: BRDF Diagram.svg|thumb|right|300px|BRDF]]
  
<math>scene = \{objects, lights, camera\}</math> - сцена, которая будет визуализирована алгоритмом, где
+
'''Энергия излучения''' - энергия, переносимая оптическим излучением. Обозначается через&nbsp; <math>Q</math>.
  
<math>objects</math> - множество объектов (например сфера, плоскость, треугольник, составной объект и т. п.).
+
'''Поток излучения''' - мощность, переносимая оптическим излучением через какую-либо поверхность.&nbsp; <math>\Phi = \frac{dQ}{dt}</math>.
  
<math>lights</math> - множество источников света (точечные источники, направленные источники, светящиеся объекты)
+
'''Плотность потока излучения''' - отношение потока излучения, проходящего через малый участок поверхности, к площади этой поверхности. &nbsp; <math>\frac{d\Phi}{dA}</math>, &nbsp; где &nbsp; <math>dA</math> - дифференциал поверхности.
  
<math>camera</math> - камера наблюдения ("виртуальный глаз").
+
'''Облученность''' - отношение потока излучения, падающего на малый участок поверхности, к площади этой поверхности. &nbsp; <math>E=\frac{d\Phi}{dA}</math>, &nbsp; где &nbsp; <math>dA</math> - дифференциал поверхности.
  
Объекты задаются произвольным образом так, чтобы можно было найти пересечение луча с этим объектом (например, для сферы достаточно задать положение в пространстве и радиус, а для плоскости - точку, принадлежащую ей и нормаль). Также, для каждого объекта необходимо задать материал - информацию о его отражающих, преломляющих, поглощающих свойствах.
+
'''Энергетическая светимость''' - отношение потока излучения, испускаемого малым участком поверхности, к площади этой поверхности.&nbsp; <math>M=\frac{d\Phi}{dA}</math>, &nbsp; где &nbsp; <math>dA</math> - дифференциал поверхности.
  
Источники света задаются по-разному, в зависимости от способа распространения света от источника. Например, точечные источники света задаются положением в пространстве и интенсивностью, направленные - интенсивностью и вектором направления.
+
'''Сила излучения''' - мощность излучения источника света в некотором направлении. Равна отношению потока излучения, исходящего от источника, к малому телесному углу. &nbsp; <math>I=\frac{d\Phi}{d\omega}</math>.
 +
 
 +
'''Энергетическая яркость''' - отношение потока излучения через малую площадку поверхности и распространяющегося в малом телесном угле, к площади проекции этой площадки на плоскость, перпендикулярную направлению распространения, и величине телесного угла.&nbsp; <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> <br>
 +
2) <math>E(\overline{p}) = \int_{\Omega_i} L_i(\overline{p}, \overline{\omega}_i)\cos \theta_i \, d\omega_i</math>,&nbsp; где интегрирование проводится по направлениям <math>\overline{\omega}_i</math> таким, что <math>(\overline{n}, \overline{\omega}_i)</math> > 0,&nbsp; то есть по полусфере.<br> (<math>\overline{n}</math> - нормаль к поверхности в точке <math>\overline{p}</math>, &nbsp; <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) > 0, \,\, (\overline{n}, \overline{\omega}_o) > 0 </math><br>
 +
Будем считать, что &nbsp; <math>dL_o(\overline{p}, \overline{\omega}_o)</math> пропорционален &nbsp; <math>dE_i(\overline{p},\overline{\omega}_i)</math>:<br>
 +
<math>dL_o(\overline{p}, \overline{\omega}_o) \sim dE_i(\overline{p},\overline{\omega}_i)</math>.<br>
 +
Функция &nbsp; <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> &nbsp; называется '''двулучевой функцией отражательной способности (BRDF)'''.<br>
 +
В точке <math>\overline{p}</math> может быть задано несколько BRDF - множество <math>\{f_{r_k}(\overline{p}, \overline{\omega}_i, \overline{\omega}_o)\}</math>.<br>
 +
Тогда итоговая 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>.<br>
 +
<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>. &nbsp; - уравнение отражений.<br>
 +
Важная часть алгоритма трассировки лучей - определить отраженную яркость по падающим яркостям, то есть посчитать этот интеграл. Обычно для вычисления интеграла применяют метод Монте-Карло.
 +
 
 +
==== Уравнение рендеринга ====
 +
 
 +
Уравнение рендеринга - это интегральное уравнение, которое определяет количество излученного в определенном направлении света как сумму отраженного и собственного излучения. Физической основой для уравнения рендеринга является закон сохранения энергии.
 +
 
 +
Уравнение отражений является частью уравнения рендеринга, поэтому остается только добавить собственное излучение объектов <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>.
 +
 
 +
По сути, алгоритмы рендеринга, в том числе и трассировка лучей, занимаются поиском приближенного решения уравнения рендеринга.
 +
 
 +
==== Поиск приближенного решения уравнения рендеринга ====
 +
 
 +
Введем оператор поиска первого пересечения.
 +
<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>.
 +
 
 +
Этот метод называется ''методом последовательных приближений''.
 +
 
 +
=== Определения ===
 +
 
 +
[[file: Ray trace diagram.svg|thumb|right|300px|Камера и видовая плоскость]]
 +
 
 +
'''Объект''' - поверхность в трехмерном пространстве.<br>
 +
'''Материал объекта''' - множество BRDF.<br>
 +
'''Сцена''' - четверка '''{камера, объекты, источники света, видовая плоскость}'''.
 +
 
 +
==== Камера ====
 +
 
 +
'''Камера''' - "виртуальный глаз", через который мы видим объекты сцены.  
  
 
Камера задается положением в пространстве и ортонормированным базисом <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>.  
 
Камера задается положением в пространстве и ортонормированным базисом <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>.  
  
=== Генерация луча ===
+
==== Видовая плоскость ====
 +
 
 +
Перед камерой на расстоянии <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> влияет на масштаб.
 +
 
 +
==== Источники света ====
 +
 
 +
Типы источников света:<br>
 +
1) '''Point light''' - точка, излучающая свет одиноково во всех направлениях.<br>
 +
2) '''Directional light''' - бесконечно удаленный источник света, задается единичным направлением на себя.<br>
 +
3) '''Area light''' - объект, излучающий свет одинаково во всех направлениях.
 +
 
 +
=== Алгоритм трассировки лучей ===
 +
 
 +
==== Генерация луча ====
  
 
Перед камерой на расстоянии <math>d</math> поставим пиксельную сетку <math>w \times h</math> с размером пикселя <math>s \times s</math> так, чтобы камера была направлена в центр сетки. Итоговое изображение будет получено из этой сетки покраской пикселей в нужный цвет.
 
Перед камерой на расстоянии <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>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}</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}.x = s (j -\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}.y = s (i -\frac{h}{2} + \frac{1}{2}) </math>
+
<math>\overline{dir}_n.z = d</math>,
  
<math>\overline{dir}.z = d</math>
+
где <math> \delta_n = (\delta_n.x, \delta_n.y) \in [0;1]^2</math> - смещения начала n-ого луча.
  
Нормируем вектор <math>\overline{dir}</math> и переведем луч в глобальную систему координат.
+
Нормируем вектора <math>\overline{dir}_n</math> и переведем лучи в глобальную систему координат.
  
=== Пересечение луча с объектами сцены ===
+
==== Пересечение луча с объектами сцены ====
  
 
Все точки луча <math>ray(\overline{p}, \overline{d})</math> можно представить в виде:
 
Все точки луча <math>ray(\overline{p}, \overline{d})</math> можно представить в виде:
Строка 53: Строка 141:
 
Если луч пересекает объект <math>O_i</math>, то <math>h_i(x,y,z) = \overline{p}+t_i \cdot \overline{d}</math> - их первое пересечение, <math>t_i</math> - расстояние от начала луча до точки пересечения.
 
Если луч пересекает объект <math>O_i</math>, то <math>h_i(x,y,z) = \overline{p}+t_i \cdot \overline{d}</math> - их первое пересечение, <math>t_i</math> - расстояние от начала луча до точки пересечения.
  
 
+
==== Вычисление энергитической яркости <math>L_o</math>, выходящей из точки пересечения ====
=== Вычисление цвета в точке пересечения ===
 
  
 
Найдем <math>T = \{t_i\}</math> - параметры пересечения с объектами (см. предыдущий пункт).
 
Найдем <math>T = \{t_i\}</math> - параметры пересечения с объектами (см. предыдущий пункт).
Строка 62: Строка 149:
 
Если <math>t_{min} = \infty</math>, то луч не пересекает ни один объект и нужно покрасить соответствующий пиксел в черный цвет (цвет фона).
 
Если <math>t_{min} = \infty</math>, то луч не пересекает ни один объект и нужно покрасить соответствующий пиксел в черный цвет (цвет фона).
  
Если <math>t_{min} < \infty</math>, то нужно найти выходящее излучение (цвет) из точки пересечения в камеру. Выходящее излучение зависит от свойств материала объекта, прямого (излучение, исходящее от источников света) и непрямого (отраженный свет) излучений, входящих в точку пересечения; угла между нормалью к поверхности в точке пересечения и направлением луча.
+
Если <math>t_{min} < \infty</math>, то нужно найти выходящее излучение (цвет) из точки пересечения в камеру. Выходящее излучение зависит от свойств материала объекта и входящего излучения в точку пересечения.
 
 
=== Вычисление прямого излучения ===
 
 
 
<math>h_{min}(x,y,z) = \overline{p}+t_{min} \cdot \overline{d}</math>.
 
 
 
<math>L_{is} = \sum_i L(S_i, h_{min})</math> - прямое входящее излучение.
 
 
 
<math>L(S_i, h_{min})</math> - излучение источника света <math>S_i</math>, приходящее в точку <math>h_{min}</math>.
 
 
 
Суммирование ведется по <math>i</math> таким, что <math>S_i</math> "виден" из точки <math>h_{min}</math>, в противном случае какой-то объект загораживает источник света, и в точку <math>h_{min}</math> будет падать тень от этого объекта.
 
 
 
=== Вычисление непрямого излучения ===
 
 
 
Рассмотрим упрощенную модель непрямого излучения. Пусть материал объекта <math>O_{min}</math> обладает способностью отражать свет.
 
 
 
Тогда нужно найти излучение, приходящее со стороны отраженного луча:
 
  
<math>ray(\overline{h}_{min}, \overline{w}_i)</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> методом Монте-Карло.
  
<math>\overline{w}_i = \overline{d} -2\cdot (\overline{n}, \overline{d})\cdot \overline{n}</math>.
+
Заметим, что каждый луч необходимо продолжить в некотором направлении(а не создавать новый пучок), так как необходимое количество лучей для численного интегрирования можно заранее сгенерировать в камере.
  
<math>L\{ray(\overline{h}_{min}, \overline{w}_i)\}</math> - излучение приходящее со стороны отраженного луча. Оно находится тем же самым алгоритмом, что и цвет пикселя (то есть рекурсивно).
+
Таким образом, алгоритм моделирует обратное движение пучка лучей света, который пришел в камеру. Пучок лучей движется из камеры и, сталкиваясь с объектами, распространяется в пространство в случайных направлениях, согласно свойствам материалов объектов.
  
 
=== Улучшение визуальных качеств изображения ===
 
=== Улучшение визуальных качеств изображения ===
  
Чтобы сгладить изображение и избавиться от шума, обычно через пиксел проводят не один луч, а достаточно много лучей, для каждого считают цвет, а потом усредняют по всем лучам (берут среднее арифметическое).
+
Заметим, что чем больше лучей выпускается через каждый пиксел изображения, тем лучше получается картинка.
  
 
== Вычислительное ядро алгоритма ==
 
== Вычислительное ядро алгоритма ==
  
Основные вычисления связаны с:
+
Основные вычисления связаны с поиском пересечения луча с объектами сцены.
  
1) поиском пересечения луча с объектами сцены.
+
== Макроструктура алгоритма ==
  
2) расчетом излучения от каждого источника света (требуется определить, "виден" ли источник света из точки пересечения с объектом).
+
Для каждого пикселя изображения генерируется пучок лучей, проходящих через этот пиксел. <br>
  
== Макроструктура алгоритма ==
+
Для каждого луча находится точка первого пересечения с объектами сцены, расчитывается входящие излучения и на основе этого - итоговое выходящее. <br>
 +
Для упрощения модели разобьем входящее излучение на ''прямое'' и ''непрямое''.<br>
 +
Прямое излучение - излучение, приходящее в точку непосредственно от источников света.<br>
 +
Непрямое излучение состоит из отраженного от других объектов сцены света.
  
Для каждого пикселя изображения генерируется пучок лучей, проходящих через этот пиксел. Для каждого луча находится первое пересечение с объектами сцены, расчитывается прямое и непрямое входящие излучения и итоговое выходящее. После этого цвет усредняется по всем лучам данного пикселя.
+
После этого цвет усредняется по всем лучам, проходящим через данный пиксел.
  
 
== Схема реализации последовательного алгоритма ==
 
== Схема реализации последовательного алгоритма ==
Строка 167: Строка 241:
  
 
//цвета
 
//цвета
direct_illumination = {0,0,0};
+
direct_illumination = {0,0,0}; //прямое излучение
indirect_illumintaion = {0,0,0};
+
indirect_illumintaion = {0,0,0}; //непрямое излучение
  
 
for (uint i = 0; i < light_count; i++)
 
for (uint i = 0; i < light_count; i++)
Строка 189: Строка 263:
  
 
== Последовательная сложность алгоритма ==
 
== Последовательная сложность алгоритма ==
 +
Вычислим последовательную сложность алгоритма в терминах пересечения луча с объектами сцены.
  
 
<math>w</math> - ширина полотна
 
<math>w</math> - ширина полотна
Строка 194: Строка 269:
 
<math>h</math> - высота полотна
 
<math>h</math> - высота полотна
  
<math>depth_{max}</math> - максимальная глубина луча в рекурсии (помимо первичного луча)
+
<math>depth_{max}</math> - максимальная глубина луча в рекурсии
  
 
<math>N_{obj}</math> - количество объектов сцены
 
<math>N_{obj}</math> - количество объектов сцены
Строка 213: Строка 288:
  
 
== Информационный граф ==
 
== Информационный граф ==
 +
 +
[[file: Raytracer_graph.png|400px]]
 +
 +
Зеленые вершины - генерация лучей для данного пикселя.<br>
 +
Желтые вершины - трассировка лучей данного пикселя.
  
 
== Ресурс параллелизма алгоритма ==
 
== Ресурс параллелизма алгоритма ==
 +
 +
Из информационного графа видно, что обработка каждого пикселя происходит независимо от других, а затем происходит соединение всех пикселей в итоговое изображение. Пусть <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>.
  
 
== Входные и выходные данные алгоритма ==
 
== Входные и выходные данные алгоритма ==
 +
 +
На вход алгоритму подается описание сцены (описание камеры, видовой плоскости, объектов и их материалов, источников света).<br>
 +
Выходные данные алгоритма - готовое изображение.
  
 
== Свойства алгоритма ==
 
== Свойства алгоритма ==
Строка 225: Строка 312:
 
== Возможные способы и особенности параллельной реализации алгоритма ==
 
== Возможные способы и особенности параллельной реализации алгоритма ==
 
== Масштабируемость алгоритма и его реализации ==
 
== Масштабируемость алгоритма и его реализации ==
 +
 +
=== Реализация с помощью MPI ===
 +
Исследование проводилось на суперкомпьютере "Ломоносов".
 +
{| class="wikitable"
 +
|-
 +
! Компилятор C++
 +
|-
 +
| Intel 13.0.1 + impi 4.1.0
 +
|-
 +
|}
 +
 +
1)<br>
 +
Размер изображения 700 на 700 пикселей, 256 лучей на пиксел, глубина рекурсии 3.<br>
 +
Сцена:<br>
 +
[[file: Raytracer_scene5.jpeg|300px]]
 +
 +
{| class="wikitable"
 +
|-
 +
! Число процессов
 +
! Время (с)
 +
|-
 +
|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)<br>
 +
Размер изображения 1340 на 1340 пикселей, 1024 лучей на пиксел, глубина рекурсии 5.<br>
 +
Сцена:<br>
 +
[[file: Raytracer_scene8.jpeg|300px]]
 +
 +
{| class="wikitable"
 +
|-
 +
! Число процессов
 +
! Время (с)
 +
|-
 +
|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-моделей).
 +
 +
=== Реализация с помощью NVIDIA CUDA ===
 +
 +
Использовалось следующее программное обеспечение:
 +
 +
{| class="wikitable"
 +
|-
 +
! OS
 +
|-
 +
| Microsoft Windows 10 Education Edition
 +
|-
 +
| Ubuntu 16.10 Yakkety Yak
 +
|}
 +
{| class="wikitable"
 +
|-
 +
! Компиляторы C++
 +
|-
 +
| nvcc из пакета NVIDIA CUDA SDK 8.0.44 + Microsoft Visual C++ 2015(Windows)
 +
|-
 +
| nvcc из пакета NVIDIA CUDA SDK 8.0.44 + gcc-4.9(Linux)
 +
|-
 +
|}
 +
 +
Результаты тестирования получены на следующих машинах:
 +
 +
{| class="wikitable"
 +
|-
 +
! 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 изображение делится на окна меньшего размера.
 +
Затем происходит отрисовка каждого окна, при этом цвета пикселей в каждом окне вычисляются параллельно.
 +
 +
[[File:CudaRaytracer_scene1.png | 500px]] [[File:CudaRaytracerPerformance.png| 500px]]
 +
 +
{| class="wikitable"
 +
|-
 +
! Лучей в окне
 +
! 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.
 +
 
== Динамические характеристики и эффективность реализации алгоритма ==
 
== Динамические характеристики и эффективность реализации алгоритма ==
 
== Выводы для классов архитектур ==
 
== Выводы для классов архитектур ==
 +
 
== Существующие реализации алгоритма ==
 
== Существующие реализации алгоритма ==
  
RenderMan
+
[https://ru.wikipedia.org/wiki/RenderMan RenderMan] — очень известный пакет программ для рендеринга 3D-анимации, который был создан компанией Pixar в 1986 году.
 +
Данный пакет использовался во множестве известных фильмов и мультфильмов, например: [https://en.wikipedia.org/wiki/Avatar_(2009_film) Аватар](2009), [https://en.wikipedia.org/wiki/Harry_Potter_and_the_Order_of_the_Phoenix_(film) Гарри Поттер и Орден Феникса](2007), [https://en.wikipedia.org/wiki/Finding_Nemo В поисках Немо](2003).
  
V-Ray
+
[https://ru.wikipedia.org/wiki/V-Ray V-Ray] — система рендеринга, разработанная компанией Chaos Group в 2000 году.
 +
VRay — это рейтрейсный рендерер, в котором присутствует несколько алгоритмов просчета глобального освещения.
  
NVIDIA IRay
+
[http://www.nvidia.com/object/nvidia-iray.html NVIDIA IRay] — технология создания фотореалистичных изображений с использованием графических ускорителей, созданная компанией NVIDIA.
  
NVIDIA mental ray
+
[https://ru.wikipedia.org/wiki/Mental_ray mental ray] — профессиональная система рендеринга и визуализации изображений, разработанная компанией mental images(является дочерней компанией NVIDIA) в 1986 году.
 +
Использован в фильмах [https://en.wikipedia.org/wiki/Terminator_3:_Rise_of_the_Machines Терминатор 3: Восстание машин](2003), [https://en.wikipedia.org/wiki/Tron:_Legacy Трон: Наследие](2010).
  
Hydra Renderer
+
[http://www.raytracing.ru/ Hydra Renderer] — это плагин для рендера реалистичных изображений в 3D Studio Max. Разработан при поддержке Лаборатории компьютерной графики и мультимедиа факультета ВМК МГУ.
  
 
= Литература =
 
= Литература =
  
http://www.scratchapixel.com/
+
1. [http://www.raytracegroundup.com/ Suffern, Kevin G. "Ray Tracing from the Ground Up / Kevin G. Suffern", 2007]
 +
 
 +
2. [http://www.scratchapixel.com/ Scratchapixel 2.0. Learn Computer Graphics From Scratch!]
 +
 
 +
3. Параллельные вычисления на GPU. Архитектура и программная модель CUDA: Учеб. пособие / А.В. Боресков и др. - М.: Издательство Московского университета, 2012. - 336 с.
  
http://www.raytracegroundup.com/
+
= Полученные изображения =
  
http://www.ray-tracing.ru/
+
<gallery>
 +
File:Raytracer_scene1.jpeg
 +
File:Raytracer_scene2.jpeg
 +
File:Raytracer_scene3.jpeg
 +
File:Raytracer_scene4.jpeg
 +
File:Raytracer_scene5.jpeg
 +
File:Raytracer_scene6.jpeg
 +
File:Raytracer_scene7.jpeg
 +
File:CudaRaytracer_scene1.png
 +
</gallery>

Текущая версия на 17:24, 7 декабря 2016

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

Содержание

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 Информационный граф

Raytracer graph.png

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

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 Входные и выходные данные алгоритма

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

1.10 Свойства алгоритма

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

2.1 Особенности реализации последовательного алгоритма

2.2 Локальность данных и вычислений

2.3 Возможные способы и особенности параллельной реализации алгоритма

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

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

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

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

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

Число процессов Время (с)
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.
Сцена:
Raytracer scene8.jpeg

Число процессов Время (с)
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 изображение делится на окна меньшего размера. Затем происходит отрисовка каждого окна, при этом цвета пикселей в каждом окне вычисляются параллельно.

CudaRaytracer scene1.png CudaRaytracerPerformance.png

Лучей в окне 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.5 Динамические характеристики и эффективность реализации алгоритма

2.6 Выводы для классов архитектур

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 Полученные изображения