\chapter{Введение в предментную область}


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


\section{Трассировка лучей}

Алгоритм трассировки лучей был впервые продемонстрирован Turner Whitted-ом в 1980-ом году \cite{Whitted:1980:IIM:358876.358882}. Этот алгоритм часто называют обратной трассировкой лучей или "Whitted-style" трассировкой. Он позволяет получить такие эффекты как тени, отражения и преломления. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{whitted}           
  \caption{\label{fig:Whitted}
           Whitted-style трассировка лучей.}
\end{figure}


\noindent Обратная трассировка лучей выглядит следующим образом: из виртуального глаза через каждый пиксел изображения испускается луч и находится точка его пересечения с поверхностью сцены (для упрощения изложения объемные эффекты вроде тумана не рассматриваются). 
Лучи, выпущенные из глаза называют первичными. Допустим, первичный луч пересекает некий объект в точке H1 (рис. ~\ref{fig:RT_Scheme1}). 

\begin{figure}[htb]
  \centering
2  \includegraphics[width=.75\linewidth]{ray_tracing_scheme}           
  \caption{\label{fig:RT_Scheme1}
           Cхема обратной трассировки лучей.}
\end{figure}


\noindent Для расчета тени необходимо определить для каждого источника освещения, видна ли из него эта точка. Предположим пока, что все источники света точечные (т.е. имеют пренебрежимо малый размер). 
Тогда для каждого точечного источника света до него испускается теневой луч из точки H1 и проверяется, пересекает ли луч какие-то объекты на своем пути. 
Если теневой луч находит пересечение с другими объектами, расположенными ближе чем источник света, значит, точка H1 находится в тени от этого источника и освещать ее не надо. Иначе, считаем освещение по некоторой локальной модели (Модель Фонга, Блина-Фонга, Кука-Торранса и.т.д.). 

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

\noindent Освещение со всех видимых (из точки H1) источников света просто складывается. Далее, если материал объекта 1 имеет отражающие свойства, из точки H1 испускается отраженный луч и для него вся процедура трассировки рекурсивно повторяется. Аналогичные действия должны быть выполнены, если материал имеет преломляющие свойства. 

\noindent В трассировке лучей для представления самого луча, как правило, используется параметрическое уравнения прямой \ref{eq:ray_param}. 

\begin{equation}\label{eq:ray_param}
\overline{X(t)} = \overline{O} + \overline{D}*t \\
\end{equation}

\noindent В выражении \ref{eq:ray_param} $ \overline{X(t)} $ обозначает некоторую точку на прямой. 
\noindent $ \overline{O} $ - точка начала (сокращение от Origin) луча, или точка вылета. $ \overline{D} $ - нормализованный (имеющий единичную длинну) вектор направления луча (сокращение от Direction). 
$t$ - скалярная величина, большая или равная нулю.

\noindent Таким образом для задания луча достаточно хранить 2 вектора - $ \overline{O} $ и $ \overline{D} $. Задания значения параметра $t$ позволит получить любую точку лежащую на этом луче.

\noindent Важным моментом при вычислении отраженного и преломленного лучей является необходимость добавлять небольшое смещение к позиции луча. 
В случае отраженного луча нужно прибавлять смещение $ \varepsilon * \overline{n} $, в случае преломленного - $(-1)*sign(dot(\overline{n},\overline{D})) * \varepsilon * \overline{n} $. 
Это делается для того, чтобы отраженный (или преломленный) луч не встретил пересечение на следующем шаге трассировки с той же самой поверхностью. В выражениях использованы следующие обозначения:

\begin{enumerate}
\item $ \overline{n}$ - вектор нормали к поверхности.
\item $ \overline{D}$ - вектор направления луча.
\item $ \varepsilon $ - некоторая малая величина, близкая к пределу машинной точности (но больше него) при вычислении точки пересечения луча и поверхности.
\item $ sign $ - функция, возвращающая 1 если ее аргумет больше или равен 0 и -1 в противном случае.
\item $ dot $ - скалярное произведение двух векторов.
\end{enumerate}

\subsection{Генерация луча}

\noindent Как уже было сказано, на начальном этапе алгоритма трассировки лучей через каждый пиксел экрана необходимо выпустить луч. 
Это можно сделать руководствуясь несложными геометрическими соображениями, рассчитав направления лучей в пространстве камеры (выражения \ref{eq:eye_ray}). 
Пространство камеры - это такое пространство, где камера находится в точке (0,0,0) и смотрит в положительном направлении оси z.

\begin{eqnarray}\label{eq:eye_ray}
fov &:=& \dfrac{\pi}{2}    \nonumber \\
D.x &:=& x + 0.5 - w/2     \nonumber \\
D.y &:=& y + 0.5 - h/2     \nonumber \\
D.z &:=& -w/tan(\dfrac{fov}{2}) \nonumber \\
D &:=& normalize(D)   
\end{eqnarray}

\noindent Где:

\begin{enumerate}
\item $ \overline{x}$ - горизонтальная координата пиксела в пространстве изображения.
\item $ \overline{x}$ - вертикальная координата пиксела в пространстве изображения.
\item $ \overline{w}$ - размер изображения по горизонтали.
\item $ \overline{h}$ - размер изображения по вертикали.
\item $ fov $ - Угол обзора камеры. Обычно равен $ \dfrac{\pi}{2} $.  
\item $ normalize $ - функция, осуществляющая нормализацию вектора (приведение к вектору с единичной длинной).
\end{enumerate}

\noindent Что касается начала луча $ \overline{O}$, то в пространстве камеры это точка (0,0,0). 

\noindent Данный подход моделирует так называемую камеру-обскура. Однако он обладает одним недостатком. 
При создании системы визуализации, в которой сцена может отрисовываться как трассировкой лучей, так и растеризацией (например при помощи библиотеки OpenGL), могут возникнуть трудности с получением геометрически-совпадающей картинки для растеризатора и трассировщика лучей. 
Причина этого заключается в том, что во всех современных системах, основанных на растеризации (и не только), используюется понятия видовой матрицы (mView) и матрицы проекции (mProj). 
Видовая матрица переводит мир из мирового пространства в пространство камеры (т.е. преобразовывает пространство так что, камера перемещается в (0,0,0) и смотрит в положительном направлении оси z). 
Эта матрица не вызывает осложнений, т.к. всегда можно преобразовать луч при помощи обратной матрицы, чтобы эмитировать перемещение камеры. 
Однако матрицу проекции необходимо пересчитать. Следующая формула используется для вычисления матрицы перспективной проекции в OpenGL \cite{Wright:2007:OSC:1406247}:

$$
 mProj  = \begin{bmatrix}
           e & 0 & 0 & 0 \\
           0 & e/a & 0 & 0 \\
           0 & 0 & -\dfrac{f+n}{f-n} & -\dfrac{2fn}{f-n} \\
           0 & 0 & -1 & 0 \\
           \end{bmatrix} 
$$


\noindent Где:

\begin{itemize}
\item e   = 1/tan(FOV/2)
\item FOV - Filed Of View; угол обзора. Обычно равен $ \dfrac{\pi}{2} $.  
\item a   - aspect ratio; a = h/w;
\item f   - far clip plane; дальняя плоскость отсечения.
\item n   - near clip plane; ближняя плоскость отсечения. 
\end{itemize}

\noindent В случае камеры-обскура, f cтремится к бесконечности, а n - к нулю. Таким образом, получаем:

$$
 mProj  = \begin{bmatrix}
           e & 0 & 0 & 0 \\
           0 & e/a & 0 & 0 \\
           0 & 0 & -1 & -2n \\
           0 & 0 & -1 & 0 \\
           \end{bmatrix}  =  \begin{bmatrix}
                             e & 0 & 0 & 0   \\
                             0 & e/a & 0 & 0 \\
                             0 & 0 & -1 & 0  \\
                             0 & 0 & -1 & 0  \\
                             \end{bmatrix}
$$


\noindent В целях обеспечения совместимости с существующими приложениями, где матрица проекции уже задана (или задается извне и менять нельзя), необходимо использовать другой алгоритм генерации луча (выражение \ref{eq:eye_ray2}):  

\begin{eqnarray}\label{eq:eye_ray2}
P1.x &:=& 2*(x + 0.5)/w - 1  \nonumber \\
P1.y &:=& -2*(y + 0.5)/h + 1 \nonumber \\
P1.z &:=& 0 \nonumber \\ 
P1.w &:=& 1 \nonumber \\
P2 &:=&(mView*mProj)^{-1} * P1 \nonumber \\
P2 &:=& P2*(1/P2.w)            \nonumber \\
P2.y &:=& -P2.y                \nonumber \\
D &:=& normalize(P2)           \nonumber \\
\end{eqnarray}


\noindent В выражение \ref{eq:eye_ray2} используется распространенная в трехмерной компьютерной графике однородная система координат \cite{Wright:2007:OSC:1406247}, 
где вводится четвертая координата $w$ в целях расширения множества преобразований, которые можно записать при помощи матриц. В соответствии с этим, матрицы $mView$ и $mProj$ имеют размер 4x4.

\subsection{Устранение ступенчастости}

\noindent Рассмотрим причину появления ступенчастости. Изображение хранится в памяти как двумерная матрица пикселей. 
Однако считается, что изображение предметов реального мира это не просто набор пикселей. 
Авторы \cite{Pharr:2010:PBR:1854996} склонны рассматривать изображение как непрерывную двумерную функцию. 
То, что мы видим на экране – приближение этой функции, ее дискретизованное в заданном разрешение представление. Ступенчатость - результат дискретизации.

\begin{figure}[htb]
  \centering
  \includegraphics[width=.75\linewidth]{antialiasing}           
  \caption{\label{fig:RT_Scheme}
           Устранение ступенчастости.}
\end{figure}

\noindent Трассировка лучей в самом общем понимании - это метод, позволяющий восстановить значения функции изображения с помощью точечных выборок, то есть значений этой функции в точках. Самый простой способ устранения ступенчатости - трассировать больше чем 1 лучей на пиксел с некоторыми смещениями внутри самого пиксела, а полученный результат усреднять. 
Однако для того чтобы уменьшить ступенчатость, достаточно часто не только трассируют больше чем 1 луч на пиксел, но и применяют фильтрацию. 
Рассмотренное выше усреднение соответствует так называемому box-фильтру. Достаточно простым, но более эффективным методом является билинейная фильтрация, при которой вклад от луча распределяется по 4 ближайшим пикселам.

\noindent Как правило процесс вычисления значения функции стараются отделить от способа хранения изображения. По этой причине во многих системах визуализации и рассчета освещения присутствуют такие классы как Sampler (отвечает за фильтрацию) и Integrator (отвечает собственно за рассчет отсвещения).

\subsection{Прозрачные объекты и тени}

\begin{figure}[htb]
  \centering
  \includegraphics[width=.75\linewidth]{blue_cow}           
  \caption{\label{fig:blue_cow}
           Цветная тень.}
\end{figure}

\noindent При расчете тени, можно использовать более сложную модель. 
Если если есть вероятность того, что один объект перекрывается другим прозрачным объектом, можно трассировть теневой луч сквозь прозрачный объект, 
расчитывая по закону Бугера-Ламберта-Бэра затухание в прозрачной среде в зависимости от пройденного расстояния (
этот же закон нужно использовать и при учете затухания внутри прозрачных объектов для обычных лучей). 
При таком алгоритме тень становится цветной. 
Разумеется, тени, рассчитанные таким образом корректны, только если прозрачный объект, отбрасывающий тень, имеет близкий к единице коэффициент преломления (считаем что коэффициент преломления воздуха равен 1). 
Если это не так, то под прозрачным объектом образуется сложная картина, называемая каустиком. 
Типичный пример каустика – солнечный зайчик от стакана воды, когда через него просвечивает солнце. 
Корректные способы расчета каустиков будут рассмотрены далее в разделе про вычисление интеграла освещенности. 

\subsection{Корректный рассчет преломлений}

\noindent В рассмотренном ранее алгоритме не обсуждалась возможность получения реалистичных изображения преломляющих объектов. Здесь необходимо отметить 2 основных момента.
Во-первых существует эффект так называемого полного внутреннего отражения при котором выходя из более плотной среды в менее плотную, 
при большом угле падения, луч полностью отражается внутрь более плотной среды и преломления не происходит. 
Во-вторых, учета только закона Снэллиуса для описания преломляющий объектов недостаточно. 
Поведение света для преломляющих объектов с достаточно-высокой степенью реалистичности описывается при помощи формул Френеля, 
которые устанавливают зависимость между долей отраженного и преломленного света в зависимости от угла падения луча и показателя преломления материала \cite{Suffern:2007:RTG:1324795}. 
Хотя следует отметить, что и этого может быть недостаточно для реалистичной визуализации некоторых кристаллов, где необходимо учитывать такие эффекты как двулучепреломление, 
поляризацию и интерференцию (на тонкой пленке оксида или другого покрытия кристалла) \cite{10.1109/TVCG.2012.304}.


\subsection{Поиск пересечений}

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

\subsubsection{Треугольник, барицентрический тест}

\begin{figure}[htb]
  \centering
  \includegraphics[width=.5\linewidth]{ray_triangle_baricentric}           
  \caption{\label{fig:ray_triangle_baricentric}
           Барицентрический тест треугольника и луча.}
\end{figure}

\noindent Введем следующие обозначения:

\begin{itemize}% начало помеченного списка
\item P - точка из который вылетает луч;
\item D - напраление луча;
\item A,B,C - вершины треугольника;
\item Тройка (u, v, 1-u-v) - барицентрические координаты. Они представляют собой отношения площадей маленьких треугольников к большому треугольнику. То есть u = S(ZCB)/S(ABC), v = S(ZAC)/S(ABC), 1-u-v = S(ZAB)/S(ABC); 
\end{itemize}% конец помеченного списка


\begin{eqnarray}\label{eq:fourierrow}
Z(u,v) &=& u*A + v*B + (1-u-v)*C \nonumber \\
Z(u,v) &=& P + t*D  \nonumber \\ 
P + t*D &=& u*A + v*B + (1-u-v)*C 
\end{eqnarray}

\noindent Имея 3 точки на плоскости, можно выразить любую другую точку через ее барицентричечкие координаты. Первое уравнение получается из определения барицентрических координат, выражая точку пересечения Z. С другой стороны, эта же точка Z лежит на прямой. Второе уравнение таким образом, это параметрическое уравнение прямой. Приравняв правые части уравнений 1 и 2 получаем третье уравнение, которое, по сути, является системой 3-х уравнений (P,D,A,B,C - векторы) с 3-мя неизвестными (u,v,t). Проведя алгебраические преобразования получим ответ в следующем виде:

\begin{align}
E1 &:= B - A \nonumber \\
E2 &:= C - A \nonumber \\
T  &:= P - A \nonumber \\
P  &:= cross(D,E2) \nonumber \\
Q  &:= cross(T, E1) 
\end{align}

$$ \begin{bmatrix}\label{eq:raytriangleformula}
t\\
u\\
v\\
\end{bmatrix} = \frac{1}{dot(P,E1)} * \begin{bmatrix}
                                      dot(Q,E2) \\
                                      dot(P,T) \\
                                      dot(Q,D) \\
                                      \end{bmatrix}
$$


\subsubsection{Треугольник, юнит-тест (Woops unit test)}

\noindent Основная идея данного алгоритма заключается в том, чтобы посчитать матрицу преобразования треугольника в некий единичный треугольник (отсюда название) с вершинами (1,0,0); (0,1,0); (0,0,0); и нормалью (0,0,1). Во время подсчета пересечения с треугольником луч преобразуется этой матрицей в пространство, где треугольник имеет единичное представление. Назовем такое проебразование $T_\Delta$. После преобразования вычислить пересечение намного проще, так как нужно считать пересечение с заранее известным треугольником - (1,0,0); (0,1,0); (0,0,0). 

\begin{figure}[htb]
  \centering
  \includegraphics[width=.5\linewidth]{unit_test}           
  \caption{\label{fig:RT_Scheme}
           Преобразование $T_\Delta$.}
\end{figure}


\noindent Пусть задан треугольник с вершинами A,B,C. Рассмотрим преобразование $T_\Delta^{-1}$:

\begin{eqnarray}\label{eq:fourierrow}
T_\Delta^{-1} = \begin{bmatrix}
                  A_x - C_x & B_x - C_x & N_x - C_x \\
                  A_y - C_y & B_y - C_y & N_y - C_y \\
                  A_z - C_z & B_z - C_z & N_z - C_z \\
                \end{bmatrix} * X + \begin{bmatrix}
                                     C_x \\
                                     C_y \\
                                     C_z \\
                                    \end{bmatrix}
\end{eqnarray}

\noindent Заметим, что применив данное преобразование к треугольнику (1,0,0); (0,1,0); (0,0,0) - получим исходный треугольник (A,B,C). То есть $T_\Delta^{-1}$ в действительности является обратным преобзазованием к $T_\Delta$


\begin{eqnarray}
T_\Delta^{-1} * \begin{bmatrix}
                  1 \\
                  0 \\
                  0 \\
                \end{bmatrix} = A & 
T_\Delta^{-1} * \begin{bmatrix}
                  0 \\
                  1 \\
                  0 \\
                \end{bmatrix} = B \\
T_\Delta^{-1} * \begin{bmatrix}
                  0 \\
                  0 \\
                  0 \\
                \end{bmatrix} = C & 
T_\Delta^{-1} * \begin{bmatrix}
                  0 \\
                  0 \\
                  1 \\
                \end{bmatrix} = N            
\end{eqnarray}

\noindent Для того чтобы найти нужное нам преобразование необходимо дополнить матрицу T$_\Delta^{-1}$ до 4x4 и (добавить еще одну строчку (0,0,0,1)) и найти обратную матрицу. Это будет соответствовать преобразованию T$_\Delta$. 

\noindent Для каждого треугольника матрицу T$_\Delta$ нужно рассчитать один раз и сохранить. Пересечение луча и треугольника в юнит-тесте вычисляется следующим образом (здесь $M$ - матрица преобразования T$_\Delta$).

\begin{align}
\overline{O} &:= mul3x4(M,\overline{O}) \nonumber \\
\overline{D} &:= mul3x3(M,\overline{D}) \nonumber \\
t &:= -o.z/d.z; \nonumber \\
u &:= o.x + t*d.x \nonumber \\
v &:= o.y + t*d.y 
\end{align}


\noindent Функция mul3x4 выполняет умножение подматрицы 3x3 на трехмерный вектор и добавляет к результату последний столбец (3 его компоненты). Функция mul3x3 просто умножает подматрицу 3x3 на трехмерный вектор. 

\subsubsection{Прямоугольный паралеллипипед (AABB)}

\noindent Прямоугольный парлеллипипед (по англ. Axis Aligned Bounding Box - или AABB) обычно представляется в виде координат 2 точек. Нижнего левого угла (третье измерение опущено) - boxMin и правого верхнего угла - boxMax. 
Все точки, координаты которыйх находятся в интервале [boxMin, boxMax] по всем 3 измерениям находятся внутри паралеллипипеда. 
Для того чтобы вычислить пересечение луча и прямоугольного паралеллипипеда сначала вычисляют точки пересечения луча со всеми 6 плоскостями паралеллипипеда (получая значения координат t в уравнении луча), 
после чего из каждой пары плоскойтей (для x,y и z) выделяют минимум и максимум (рис. ~\ref{fig:rayboxintersection}). 
Результирующие координаты ищут как максимальный из минимумов ($t_{min}$) и минимальный из максимумов ($t_{max}$) - выражение \ref{eq:rayboxintersect}.

\begin{figure}[htb]
  \centering
  \includegraphics[width=.75\linewidth]{raybox}           
  \caption{\label{fig:rayboxintersection}
           Пересечение луча и AABB. Точки m1 и m2 - минимумы пересечения луча с парой плоскостей. M1 и M2 - максимумы.}
\end{figure}


\begin{align} \label{eq:rayboxintersect}
l_x &:= (1/D.x)*(boxMin.x - O.x); \nonumber \\
h_x &:= (1/D.x)*(boxMax.x - O.x); \nonumber \\
l_y &:= (1/D.y)*(boxMin.y - O.y); \nonumber \\
h_y &:= (1/D.y)*(boxMax.y - O.y); \nonumber \\
l_z &:= (1/D.z)*(boxMin.z - O.z); \nonumber \\
h_z &:= (1/D.z)*(boxMax.z - O.z); \nonumber \\
m1 &:= min(l_x, h_x);             \nonumber \\
M1 &:= max(l_x, h_x);             \nonumber \\
m2 &:= min(l_y, h_y);             \nonumber \\
M2 &:= max(l_y, h_y);             \nonumber \\
m3 &:= min(l_z, h_z);             \nonumber \\
M3 &:= max(l_z, h_z);             \nonumber \\
t_{min}  &:= max(m1, m2, m3)      \nonumber \\
t_{max}  &:= min(M1, M2, M3)
\end{align}


\subsection{Ускорение поиска пересечений}


\subsubsection{Регулярные сетки}

\noindent Ускорение поиска пересечений в трассировке лучей - довольно обширная тема. 

\noindent Если сцена представлена набором примитивных объектов (сфера, AABB, треугольник), то ускорять процесс поиска пересечений при помощи структур пространственного разбиения имеет смысл, 
как правило, если число объектов больше одного-двух десятков. В противном случае, стоимость поиска может быть сравнима со стоимостью перебора всех объектов вслепую. 
Как правило, это происходит из-за большего количеста ошибок предсказателя ветвлений при выполнении сложного кода (особенно верно для иерархического поиска в дереве).

\noindent Данная ускоряющая структура регулярно разбивает пространство на $N^3$ вокселей, где N - разрешение сетки. 
При поиске в регулярной сетке проверяются только те объекты, которые оказались внутри вокселей на пути луча.

\begin{figure}[htb]
  \centering
  \includegraphics[width=.75\linewidth]{regular_grid}           
  \caption{\label{fig:regular_grid}
           Регулярное разбиение пространства.}
\end{figure}


\noindent Алгоритм построения сетки сводится к нахождению геометрических пересечений между вокселами и объектами сцены.
Каждый воксел должен содержать список всех указателей или индексов объектов с которыми он имеет пересечение.  

\noindent Алгоритм поиска носит фамилию своего автора Фуджимото \cite{Fujimoto:1988:AAR:95075.95111}. Функция поиска состоит из 2 частей. Инициализация и перебор вокселей. 
Инициализирующая часть весьма объемная, однако, она выполняется для каждого луча 1 раз. Алгоритм \ref{fujimoto_trav_alg} и рисунок \ref{fig:fujimoto_traverse} демонстрируют поиск Фуджимото.

\begin{figure}[htb]
  \centering
  \includegraphics[width=.75\linewidth]{fujimoto_traverse}           
  \caption{\label{fig:fujimoto_traverse}
           Поиск в регулярной сетке.}
\end{figure}

\begin{algorithm}[H]
\SetAlgoLined
\KwData{Луч - ray, ограничивающий куб - box}
\KwResult{Перебор вокселей, лежащих на пути луча }

Инициализирующая часть алгоритма:\\

$ (tMaxX,tDeltaX,stepX) \leftarrow $ ВычислитьЗначения(ray,box)\;
$ (tMaxY,tDeltaY,stepY) \leftarrow $ ВычислитьЗначения(ray,box)\;

Основная часть алгоритма:\\
\While{ (x < N) and (x >= 0) and (y < N) and (y >= 0) }{
ОбработатьВоксел(x,y)\;
\eIf{$tMaxX < tMaxY$}{
  $tMaxX \leftarrow tMaxX + tDeltaX$\;
  $x \leftarrow x + stepX$\;
}{
  $tMaxY \leftarrow tMaxY + tDeltaY$\;
  $y \leftarrow y + stepY$\;
 }
}
\caption{Обход регулярной сетки Fujimoto}
\label{fujimoto_trav_alg}
\end{algorithm}

\noindent Несмотря на простой алгоритм перебора вокселей, регулярная сетка - не очень хорошая ускоряющая структура. Причин этому несколько:

\begin{enumerate} % начало помеченного списка
\item \underline{Высокие затраты памяти}, пропорционально $N^3$.
\item \underline{Отсутствие адаптивности} или проблема чайника на стадионе. Регулярная сетка одинаково разбивает все пространство независимо от того, какая геометрия и как расположена на сцене. Размер воксела выбирается исходя из некоторого среднего размера объектов. При этом если выбрать размер сетки большим, поиск пересечений значительно замедлится для небольших высокополигональных объектов. А если размер воксела выбрать маленьким, не только значительно возрастут затраты памяти но и замедлится трассировка на всей остальной сцене, т.к. в процессе трассировки луча придется перебирать большее число вокселей.
\item \underline{Проблема повторных пересечений.} Эта проблема также усугубляется с ростом разрешения сетки. Объекты, пересекающие несколько вокселей при поиске пересечений будут встречаться несколько раз. Следовательно, при трассировке в регулярной сетке луч будет вычислять пересечение с такими примитивами более чем один раз. Проблему можно частично амортизировать сохраняя идентификатор луча в примитивах, с которым пересечение уже было посчитано. Если в данный момент времени идентификатор луча равен идентификатору сохраненному в примитиве, вычислять пересечение не нужно. Однако этот подход трудно применять при параллельной реализации трассировки лучей.
\end{enumerate}% конец помеченного списка


\subsubsection{Окто-деревья и иерархические сетки}

\noindent Перейдя к иерархическому представлению, можно в некоторой степени решить первые 2 проблемы регулярной сетки (затраты памяти и неадаптивность). Иерархическая сетка представляет из себя дерево (обычно небольшой глубины, но зато очень широкое), в узлах которого лежат относительно небольшие регулярные сетки. Если размер сетки положить равным 2x2x2 получаем окто-дерево. Таким образом, окто-дерево представляет из себя иерархическое разбиение пространства на 8 частей (рис. ~\ref{fig:octree2}).

\begin{figure}[htb]
  \centering
  \includegraphics[width=.75\linewidth]{octree2}           
  \caption{\label{fig:octree2}
           Окто-дерево.}
\end{figure}

\noindent Окто-деревья широко используются в компьютерной графике для различных задач пространственного поиска, но они довольно плохо подходят для ускорения трассировки лучей. 
Алгоритм обхода окто-деревьев достаточно сложен. При этом, адаптивность у окто-дерева низкая, так как чтобы ограничить небольшой объект (чайник на стадионе) нужно много уровней подразбиения. 
Окто-дерево, таким образом, обычно содержит множество пустых узлов и сложные регионы обходятся медленно (рис. ~\ref{fig:bunny_octree_32}). 
При этом, проблема повторных пересечений по прежнему остается актуальной для окто-деревьев, поскольку при разбиение один объект очень часто попадает сразу в несколько дочерних узлов.

\begin{figure}[htb]
  \centering
  \includegraphics[width=.75\linewidth]{bunny_octree_32}           
  \caption{\label{fig:bunny_octree_32}
           Модель стэнфордского кролика в окто-дереве. Количество белого цвета на изображении визуально позволяет оценить сложность трассировки. }
           
\end{figure}

\subsubsection{kd-деревья}

\noindent Рассмотрим структуру бинарного пространственного разбиения, называемую kd-дерево (аббревиатура kd расшифровывается как k-dimentional). Эта структура представляет собой бинарное дерево ограничивающих параллелепипедов, вложенных друг в друга. Каждый параллелепипед в kd-дереве разбивается плоскостью, перпендикулярной одной из осей координат на два дочерних параллелепипеда.

\noindent Вся сцена целиком содержится внутри корневого параллелепипеда, но, продолжая рекурсивное разбиение параллелепипедов, можно прийти к тому, что в каждом листовом параллелепипеде будет содержаться лишь небольшое число примитивов. Таким образом, kd-дерево позволяет использовать бинарный поиск для нахождения примитива, пересекаемого лучом. 

\noindent Если плоскость, разбивающую пространство выбирать каждый раз по середине паралеллипипеда (при этом текущий узел можно разбивать всегда по оси, в которой он имеет максимальный размер), kd-дерево будет эквивалентно окто-дереву и наследует все его недостатки (за исключением сложного алгоритма обхода). Однако, основное преимущество kd-дерева над окто-деревом заключается как раз в том, что плоскость разбиения можно ставить не только по середине разбиваемого узла, а в любом месте. 

\noindent Алгоритм построения kd-дерева можно представить следующим образом (будем называть прямоугольный параллелепипед сокращенно AABB).  

\begin{enumerate}% начало помеченного списка

\item 'Добавить' все примитивы в ограничивающий AABB. Т.е построить ограничивающий все примитивы AABB, который будет соответствовать корневому узлу дерева.
\item Если примитивов в узле мало или достигнут предел глубины дерева, завершить построение.
\item Выбрать плоскость разбиения, которая делит данный узел на два дочерних. Будем называть их правым и левым узлами дерева.
\item Добавить примитивы, пересекающиеся с AABB левого узла в левый узел, примитивы, пересекающиеся с AABB правого узла в правый.
\item Для каждого из узлов рекурсивно выполнить данный алгоритм начиная с шага 2. 
 
\end{enumerate}% конец помеченного списка

\noindent На рисунке ~\ref{fig:kdtreelr} изображен процесс построения kd-дерева с учетом некоего оптимального алгоритма выбора разбивающей плоскости. Самым сложным в построении kd-дерева является 3-ий шаг. От него напрямую зависит эффективность ускоряющей структуры. Существует несколько способов выбора плоскости разбиения, рассмотрим их по порядку.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{kd_tree1}           
  \caption{\label{fig:kdtreelr}
           Построение kd-дерева. L - обозначает левое поддерево. R - правое.  }
           
\end{figure}

\paragraph{Разбиение по середине.} \noindent

\noindent Самый простой способ - выбирать плоскость разбиения по центру. Сначала выбираем ось (x, y или z), в которой AABB имеет макимальный размер, затем разбиваем его по центру. 
Ранее обсуждалось, что kd-дерево в этом случае имеет плохую адаптивность и эквивалентно окто-дереву.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{kd_tree_split_center}           
  \caption{\label{fig:kd_tree_split_center}
           Разбиение по середине. }
           
\end{figure}

\paragraph{Разбиение по медиане.} \noindent

\noindent Если разбивать узел на два дочерних таким образом, чтобы в правом и левом поддереве количество примитивов было одинаково, будет построено сбалансированное дерево. 
Это не очень удачная идея. Все дело в том, что сбалансированные деревья могут помочь только если искомый элемент каждый раз находится в случайном узле, то есть если распределение лучей по узлам во время поиска будет равномерно. В действительности, это не так. Лучей в среднем пойдет больше в тот узел, который больше по своей площади поверхности, а при медианном разбиении эти площади у узлов могут быть разные. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{kd_tree_split_median}           
  \caption{\label{fig:kd_tree_split_median}
           Разбиение по медиане. }
           
\end{figure}

\paragraph{Разбиение на основе оптимизации функции стоимости.} \noindent

\noindent Каковы же критерии хорошо-построенного kd-дерева? На интуитивном уровне такой критерий можно описать фразой 'как можно больше пустого пространство должно быть отброшено как можно быстрее'. 
Для решения поставленной задачи используем формальный подход. Введем функцию стоимости поиска, которая будет отражать, насколько дорого по вычислительным ресурсам производить поиск в данном узеле случайным набором лучей. 
Будем вычислять ее по следующей формуле \cite{MacDonald:1990:HRT:87851.87856}:

\begin{align*}
  SAH(x) = CostEmpty + SA(Left)*N(Left) + SA(Right)*N(Right)
\end{align*}

\noindent В приведенной выше формуле SA(Left) и SA(Right) - площадь поверхности (SurfaceArea) соответственно левого и правого дочерних узлов, получающихся при разбиении. 
N(Left) и N(Right) - количество примитивов, попавших при разбиении соответственно в левое и правое поддерево. Аргумент функции x является одномерной координатой плоскости разбиения.  
На рисунке ~\ref{fig:kd_tree_split_sah} можно увидеть, что SAH сразу отбрасывает большие пустые пространства, плотно ограничивая геометрию.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{kd_tree_split_sah}           
  \caption{\label{fig:kd_tree_split_sah}
           Разбиение с учетом оптимизации стоимости. }
           
\end{figure}

\noindent Хорошими кандидатами на минимум SAH могут служить границы примитивов. Простой алгоритм посторения выглядит следующим образом: каждый раз при выборе плоскости нужно перебрать всевозможные границы примитивов по трем измерениям, вычислить в них значение функции стоимости и найти минимум среди всех этих значений. 
Когда мы вычисляем SAH для каждой плоскости, то нам необходимо знать N(left) и N(right) - количества примитивов справа и слева от плоскости. 
Если вычислять N простым перебором, в итоге получится квадратичный по сложности алгоритм построения.

\noindent Действенным способом ускорения поиска минимума функции SAH является использование какого-либо метода минимизации одномерной функции с несколькими начальными приближениями. 
Например, метод золотого сечения. Это избавляет от необходимости полного квадратичного перебора. Также популярен метод называемый binning \cite{popov:06:ESC, Wald:2007:FCS:1524874.1524953}, 
который ускоряет поиск за счет снижения точности поиска положения найденного минимума.


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{kd_tree_chair}           
  \caption{\label{fig:kd_tree_chair}
           kd-дерево, построенное с учетом SAH. }
           
\end{figure}


\paragraph{Критерий остановки.} \noindent

\noindent Одним из критериев остановки может служить сама функция SAH. Если оцениваемая функцией суммарная стоимость поиска в дочерних узлах больше стоимости поиска в родительском узле, разбиение стоит остановить. 
Однако, весьма неочевидным фактом является то, что такой критерий остановки далеко не оптимален. 
Дело в том, что увеличение SAH при разбиении текущего уровня дерева еще не означает замедление трассировки для всего дерева. 
Рассмотрим цилиндр, составленный из полигонов на рисунке ~\ref{fig:tricilinder}. Узел kd-дерева, обозначенный буквой A является проблемным. 
Какая-бы плоскость разбиения не была выбрана, SAH полученных дочерних узлов будет больше чем SAH родительского узла. 
Таким образом, SAH в данном случае говорит о том, что подразбиение нужно остановить. Но это не всегда так.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{tricilinder}           
  \caption{\label{fig:tricilinder}
           Цилиндр, составленный из треугольников. Какую бы плоскость разбиения не была выбирана для узла A, это не позволит уменьшить SAH. }   
\end{figure}


\noindent Если представить, что цилиндр - высокополигональный объект, дальнейшее подразбиение (даже с учетом выбора плоскости по центру) приведет к тому, что в листьях окажется не более чем 1-2 примитива. 
Однако остается открытым вопрос - будет ли при этом быстрее трассировка лучей и если да, то на сколько. И очередной возникающий вопрос формулируется так: 
Если существуют такие плоскости разбиения, которые увеличиват суммарный SAH текущего узла (который может быть оценен как SAH(левого) + SAH(правого)), но при этом все-же позволяют ускорить трассировку, когда именно следует выбирать такие плоскости (т.е. на каком уровне дерева)?

\noindent Техника, позволяющая решить описанную выше проблему и ответить на поставленные вопросы называется в англоязычной литературе 'Early Split Clipping' \cite{Ernst:2007:ESC:1524874.1524958}. 
Будем называть ее ранним подразбиение примитивов. Она говорит, что примитивы, которые плохо аппроксимируются при помощи AABB нужно подразбивать на несколько более мелких примитивов еще до построения дерева с тем, 
чтобы рассматривать такие 'сложные' плоскости на равне со всеми остальными в процессе оптимизации SAH. Подробнее данная техника будет рассмотрена при обсуждении BVH деревьев.

\paragraph{Поиск в kd-дереве.} \noindent

\noindent Классический алгоритм бинарного поиска в kd деревьях (kd-tree traversal в англоязычной литературе), состоит в следующем: 
На первом шаге алгоритма необходимо посчитать пересечение луча с ограничивающим сцену корневым параллелепипедом (AABB) и запомнить информацию о пересечении в виде двух координат (в пространстве луча) – t\_near и t\_far, обозначающих пересечение с ближней и дальней плоскостями соответственно. На каждом следующем шаге необходима информация только о текущем узле (его адрес) и этих двух координатах.

\noindent При этом нет необходимости вычислять пересечение луча и дочернего параллелепипеда, достаточно лишь узнать пересечение с разбивающей параллелепипед плоскостью (обозначим соответствующую координату как t\_split). 

\noindent В случае если t\_split >= t\_far (рис. \ref{fig:kdtree_a}) или если t\_split < t\_near (рис. \ref{fig:kdtree_b}) луч пересекает только один дочерний узел, поэтому можно просто отбросить правый (соответственно левый) узел и продолжить поиск пересечения в оставшемся узле.

\noindent В случае если луч пересекает оба дочерних узла (рис. \ref{fig:kdtree_c}), необходимо сначала поискать пересечение в ближнем узле и если оно не найдено, искать его в дальнем. 
Так как в общем случае неизвестно, сколько раз произойдет последнее событие(отсутствие пересечения в ближнем узле), необходим стек. 
Каждый раз, когда луч пересекает оба дочерних узла, адрес дальнего узла, t\_near и t\_far помещаются в стек и поиск продолжается в ближнем. 
Если в ближнем узле пересечение не найдено, из стека достаются адрес дальнего узла, t\_near, t\_far и поиск продолжается в дальнем узле.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{kdtree_a}           
  \caption{\label{fig:kdtree_a}
           Поиск в kd-дереве. (t\_split >= t\_far ). }   
\end{figure}

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{kdtree_b}           
  \caption{\label{fig:kdtree_b}
           Поиск в kd-дереве. (t\_split < t\_near ). }   
\end{figure}

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{kdtree_c}           
  \caption{\label{fig:kdtree_c}
           Поиск в kd-дереве. (t\_near <= t\_split <= t\_far ). }   
\end{figure}


\noindent Почти все рассмотренные ранее проблемы регулярных сеток kd-дерево успешно решает. 
Оно обладает простым алгоритмом поиска, хорошей адаптивностью, быстро отбрасывая пустые пространства и значительно амортизирует проблему повторных пересечений. 
Однако серьезный недостаток kd-дерева - сложный алгоритм построения. Дело здесь не только в необходимости вычисления функции стоимости SAH, но, что гораздо боле важно, перед началом построения kd-дерева невозможно сказать заранее, сколько памяти потребуется для его построения. 
Поскольку при разбиении узла все примитивы теоретически могут попасть как в левый, так и в правый дочерний узел, при разбиении узла с числом примитивов N необходимо иметь 2N свободной памяти.

\subsubsection{BSP-деревья}

\noindent В отличие от кd-деревьев, в которых плоскости разбиения имеют всего 3 возможных ориентации (XoY, XoZ, YoZ), в BSP (Binary Space Partion) деревьях плоскости разбиения могут быть ориентированны произвольным образом \ref{fig:bsp2}. Это значительно усложняет выбор плоскости с учетом оптимизации SAH (и как следствие сам алгоритм построения), но частично решает проблему обозначенную на рисунке ~\ref{fig:tricilinder}. В работе \cite{ize:08:BSP} было показано преимущество BSP-деревьев в скорости поиска над kd-деревьями. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{bsp2}           
  \caption{\label{fig:bsp2}
           kd-дерево (слева); BSP-дерево (справа). }   
\end{figure}

\noindent Тем не менее, даже BSP дерево не всегда гарантирует уменьшение числа примитивов в дочерних узлах (рис. \ref{fig:bsp2}, узел с тремя треугольниками справа вверху). На практике из этого следует, что проблема с неизвестным количеством памяти, требующимся для построения дерева, все еще актуальна для BSP деревьев. 

\subsubsection{BVH-деревья}

\noindent Этот вид ускоряющей структуры является наиболее практичным и часто используемым в современных трассировщиках лучей. BVH расшифровывается как Bounding Volume Hierarchy - Иерархия Ограничивающих Объемов. BVH-дерево состоит из вложенных друг в друга объемов, ограниченных некоторыми примитивами (рис. \ref{fig:bvh1}). Такие деревья можно классифицировать по типу ограничивающего примитива. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{bvh1}           
  \caption{\label{fig:bvh1}
           BVH-дерево. В отличие от других видов пространственных деревьев, узлы дерева могут пересекаться друг с другом. }   
\end{figure}

\noindent В дальнейшем будем рассматривать только деревья, использующие прямугольные паралеллипипеды (Axis Aligned Bounding Box, AABB) в качестве ограничивающих примитивов и под BVH деревом понимать дерево состоящее из вложенных друг в друга прямоугольных паралеллипипедов, стороны которых выровнены по осям координат. 

\paragraph{Две стратении разбиения в BVH дереве} \noindent 

\noindent При построении дерева существует небольшая, но довольно существенная путаница в том, какие примитивы включать в дочерние узлы и как именно производить разбиение. 
Вопрос, который порождает путаницу звучит так: "Что делать, если примитив пересекает оба дочерних узла?". Существует 2 стратегии разбиения.

\noindent Первая стратегия называется \underline{разбиением по объектам} (рис. \ref{fig:bvh_object_split}). 
При использовании этой стратегии необходимо найти некоторое разбиение нашего исходного множества примитивов на 2 подмножества. Затем достаточно лишь посчитать ограничивающий обьем для обоих подмножеств. 
В этом случае ответ на поставленные ранее вопрос звучит так: Такой примитив всегда идет в тот узел, в котором он содержится полностью. Такой узел всегда существует по построению указанного разбиения. 

\noindent Вторая стратегия называется \underline{пространственным разбиением} (рис. \ref{fig:bvh_object_split}). 
При использовании этой стратегии сначала выбирается некоторая плоскость, которая делит объем текущего узла на 2 части. Для полученных дочерних объемов ограничивающие AABB пересчитываются в соответствии с геометрией. 
В этом случае ответ на поставленный вопрос звучит так: "Такой примитив нужно добавлять в каждый из дочерних узлов, который он пересекает".

\noindent Рассмотренные стратегии порождают 2 принципиально разных вида BVH дерева. Стратегия разбиения по объектам позволяет строить такое дерево, в котором для каждого объекта хранится не более 1 ссылки. 
Будем называть такие дереьвя деревьями с одиночными ссылками (single reference trees). 

\noindent В противоположность этому, стратегия пространственного разбиения  может породить несколько ссылок на 1 объект. 
Будем называть такие деревья деревьями со множественными ссылками (multiple reference trees).


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{bvh_object_split}           
  \caption{\label{fig:bvh_object_split}
           Разбиение по объектам (слева) и пространственное разбиение (справа). Тонкий длинный треугольник справа при стратегии разбиения по объектам попадет только в левый узел L. А при стратегии пространственного разбиения он попадет в оба дочерних узла L и R.}   
\end{figure}

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

\paragraph{Построение BVH дерева.} \noindent

\noindent Алгоритм построения оптимального BVH дерева столь-же сложен, сколь и алгоритм построения оптимального kd-дерева. В теории даже сложнее. 
Когда мы строим kd-дерево, то на каждом шаге подразбиения минимизируем функцию стоимости SAH. 
В случае kd-дерева мы выбираем только плоскость разбиения и можем перебрать все возможные границы примитивов по трем измерениям (обычно SAH вычисляют на границах), чтобы найти такую плоскость, в которой функция SAH будет минимальна. 
То есть в худшем случае, при отсутствие оптимизаций нам нужно вычислить SAH лишь 6*N раз (где N - число примитивов) и найти минимум.

\noindent В случае BVH всё сложнее, поскольку вариантов возможных разбиений не 6*N а $S(N,2) = 2 ^{N-1} - 1 $ (где S - число Стирлинга 2 рода). 
Перебрать их все не представляется возможным. В отличие от kd-дерева, BVH на каждом уровне вообще говоря может произвести разбиение на 2 совершенно произвольных подмножества примитивов. Всего таких подмножеств $2 ^{N-1} - 1$. 

\noindent Однако можно использовать упрощение, не сильно ухудшающее качество дерева, но позволяющее значительно упростить процедуру его построения. Будем строить BVH дерево сверху вниз. 
При построении, в каждом узле дерева сортируем примитивы по минимумам (или центрам) их AABB и 3 осям координат. Таким образом, получаем 3 массива. 
Первый отсортирован по X координатам минимумов (или центров) AABB объектов, второй по их Y координатам и третий по их Z координатам. 
Затем для каждого из полученных массивов (длинной n) необходимо перебрать все разбиения на подмножества вида:

\begin{align*}
[0     ; 1..n] \\
[0..1  ; 2..n] \\
[0..3  ; 4..n] \\
...            \\
[0..n-1; n]    \\
\end{align*}

\noindent Для каждого из этих разбиений вычисляем значение SAH. Запоминая ограничивающие AABB в специальном массиве для получаемых подмножеств мы можем не пересчитывать эти AABB полностью, 
что сводит поиск минимума SAH всего к 2 проходам по массиву. С конца к началу и с начала к концу. Описанный выше подход используется в реализации трассировки лучей из работы \cite{Aila:2009:UER:1572769.1572792}.

\noindent Среди эффективных подходов, позволяющих значительно сократить время построения BVH дерева для анимированных сцен следует отметить подход из работы \cite{journals/cgf/Garanzha09}.
За счет кластеризации входных данных, в работе \cite{journals/cgf/Garanzha09} достигнуто уменьшение объема обрабатываемых построителем примитивов и, как следствие, значительное ускорение построения BVH дерева.

\subsubsection{Резюме по ускоряющим структурам}

\noindent В сочетании с ранним подразбиением, BVH дерево на практике показало себя с наилучшей стороны \cite{Aila:2009:UER:1572769.1572792}. Алгоритм построения BVH достаточно прост и занимает определенный объем памяти, пропорционально количеству входных примитивов.

\section{Проблема глобальной освещенности}

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

\noindent Один из способов получения фотореалистичных изображений - моделирование оптических процессов с целью точного вычисления освещенности каждой точки сцены и.
Как следствие, появляется необходимость в вычислении глобальной освещенности. 
Трассировка лучей, описанная выше, позволила бы рассчитать изображение корректно, если бы поверхности объектов реального мира были гладкими и отражали свет строго по закону 'угол падени равен углу отражения'. 
Однако, из-за наличия микрорельефа поверхностей это не так, и свет приходящий строго с одного направления, распределяется по полусфере вокруг нормали к поверхности в соответствии со свойствами материала. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{brdf}           
  \caption{\label{fig:brdf}
           Двунаправленная Функция Отражения.}   
\end{figure}

\noindent Для моделирования свойств материала вводят Двунаправленную Функцию Отражения (ДФО) \cite{BRDFOriginal:65}. В английской терминологии обозначается как BRDF (Bidirectional Reflectance Distribution Function). 
ДФО связывает по некоторому закону интенсивность и угол падающего света с интенсивностью и углом отраженного поверхностью света. 
Зная ДФО материала, можно вычислить освещенность точки поверхности при помощи интеграла освещенности \cite{Kajiya:1986:RE:15886.15902}:

\begin{eqnarray}\label{eq:integral}
I(\phi_{r}, \theta_{r}) = \iint\limits_{\phi_{i}\theta_{i}} L(\phi_{i},\theta_{i})R(\phi_{i},\theta_{i},\phi_{r},\theta_{r})cos(n,l_{\phi_{i},\theta_{i}}) \mathrm{d}\phi_{i}  \mathrm{d}\theta_{i}
\end{eqnarray}

\noindent В уравнении \ref{eq:integral} функция $ L(\phi_{i},\theta_{i}) $ задает яркость света, падающего с направления задаваемого вектором $ l_{\phi_{i},\theta_{i}} $. $ R(\phi_{i},\theta_{i},\phi_{r},\theta_{r}) $ - ДФО. 
$n$ - вектор нормали к поверхности. Вычисляемая интегрированием интенсивность освещения в точке $ I(\phi_{r}, \theta_{r}) $ является не числом, а функцией. 
Другими словами, она дает значения интенсивности света, отражаемой поверхностью под разными углами. 
Таким образом, выполняя интегрирование по полусфере приходящего (падающего) в точку освещения $L$ c учетом свойств поверхности $R$, получаем новую функцию распределения освещения в пространстве $I$, 
обусловленную свойствами отражения поверхности (материала) в некоторой точке x. 
Если зафиксировать конкретный луч (например выпущенный для конкретного пиксела и виртуальной камеры), направление $ \phi_{r},\theta_{r} $ было бы фиксированно и задавалось бы этим лучом. 

\noindent Формула \ref{eq:integral} берется их тех соображений, что освещенность в точке поверхности складывается из света, падающего на поверхность со всех направлений и отражающегося в заданном направлении по определенному закону (который задает ДФО). Эта формула не более чем математическая модель и она верна не всегда. Например, в случае стекла нужно учитывать свет, приходящий из под поверхности и проводить интегрирование по полной сфере. 
В случае кожи ситуация еще сложнее, так как необходимо учитывать такой эффект как подповерхностное рассеивание \cite{Nicodemus:1992:GCN:136913.136929, Jensen:2002:RHR:566654.566619}.

\noindent Размерность интеграла освещенности (формула \ref{eq:integral}) зависит от количества учитываемых переотражений света, поскольку функция L в некоторой точке x будет зависеть от функции I в некотрой другой точке y. Из-за большой размерности интеграла для его вычисления обычно используют метод Монте-Карло, поскольку его скорость сходимости не зависит от размерности интеграла.

\subsection{Метод Монте-Карло}

\noindent Пусть u - случайная величина, равномерно распределенная на отрезке [a,b]. Тогда f(u) - тоже случайная величина. Исходя из определения, ее математическое ожидание равно:

\begin{eqnarray}\label{eq:montecarlo}
Ef(u) = \int_a^b f(x)p(x)dx \nonumber \\
p(x) = \dfrac{1}{b-a}
\end{eqnarray}

\noindent Где $p(x)$ - плотность вероятности равномерного распределения. Выражая $ \int_a^b f(x)dx $ из формулы \ref{eq:montecarlo} получаем:

\begin{eqnarray}\label{eq:montecarlo2}
\int_a^b f(x)dx = (b-a)Ef(u) \approx \dfrac{b-a}{N} \sum_{i=1}^N f(u_{i})
\end{eqnarray}

\noindent В применении к интегралу освещенности (имеем право взять полусферу единичной площади):

\begin{eqnarray}\label{eq:montecarlo3}
I(\phi_{r}, \theta_{r}) \approx \dfrac{1}{\pi} \dfrac{ \sum_{i=1}^N L(\phi_{i},\theta_{i})R(\phi_{i},\theta_{i},\phi_{r},\theta_{r})cos(n,l_{\phi_{i},\theta_{i}}) }{N}
\end{eqnarray}

\noindent Итак, мы можем оценивать значения интеграла при помощи вычисления суммы достаточно большого числа значений под-интегральной функции \cite{Sobol:73}. 
Далее будут рассмотрены некоторые определения и понятия математической статистики. 

\subsubsection{Выборка по значимости (Importance Sampling).}

\noindent Эффективным способов улучшения сходимости метода Монте-Карло является выборка по значимости (метод существенной выборки, importance sampling) \cite{Sobol:73}. 
Данный метод предполагает использование случайной величины не с равномерным распределением, а с таким распределением, которое пропорционально модулю вычисляемой функции. Идея выборки по значимости базируется на том, что некоторые значения случайной величины в процессе моделирования имеют большую значимость для оцениваемой функции, чем другие. Если эти 'более вероятные' значения будут появляться в процессе выбора случайной величины чаще, дисперсия оцениваемой функции уменьшится и скорость сходимости, таким образом, увеличится \cite{Sobol:73}. При этом, для того чтобы метод давал корректный результат, необходимо  учитывать получаемые значения с весами, обратно-пропорциональными плотности вероятности случайной величины. Таким образом, если распределение случайных величин не равномерное, оценка для метода Монте-Карло запишется в более общем виде \cite{Sobol:73}:

\begin{eqnarray}\label{eq:montecarlo4}
\int_a^b f(x)dx \approx \dfrac{1}{N} \sum_{i=1}^N \dfrac{f(u_{i})}{p(u_{i})}
\end{eqnarray}

\noindent где $ p(u_{i}) $ - функция плотности вероятности случайной величины u.

\subsubsection{Смещенные и несмещенные оценки и методы.}

\noindent В данном разделе мы кратко рассмотрим понятие смещенных и несмещенных оценок в применении к вычислению интеграла освещенности чтобы пояснить такие часто встречающиеся англо-язычные термины как biased renderer и unbiased renderer.

\noindent Пусть $ X_{1} .. X_{n} $ случайная выборка из распределения, зависящего от параметра $ \theta \in \Theta $

\noindent \underline{Статистикой} называется любая измеримая функция от выборки $\theta_{i}$, не зависящая от $\theta$ (важно, что значение статистики можно вычислить при каждой реализации выборки, не зная значения $\theta$). В случае метода Монте-карло под-интегральное выражение может называться статистикой.

\noindent \underline{Точечной оценкой} называется любая статистика, принимающая некоторые знаения на области определения $ \Theta $.

\noindent Оценка называется \underline{несмещенной} (unbiased) если ее математическое ожидание равно оцениваемому параметру. В противном случае оценка называется \underline{смещенной} (biased).

\noindent Оценка называется \underline{состоятельной} (consistent) если она сходится по вероятности (сходимость почти наверное) к оцениваемому параметру. 

\noindent Будем называть метод или алгоритм вычисления интеграла освещенности \underline{несмещенным}, если он позволяет получить несмещенную оценку интеграла освещенности для индивидуальной точки или набора точек трехмерного пространства. 
В противном случае будем называть метод \underline{смещенным}. 

\noindent Будем называть метод вычисления интеграла освещенности \underline{состоятельным}, если он позволяет получить состоятельную оценку интеграла освещенности для индивидуальной точки или набора точек трехмерного пространства.

\noindent Таким образом, если программное обеспечение для вычисления освещенности и фотореалистичной визуализации (англ. renderer) использует в том или ином виде метод Монте-Карло и при этом позволяет при 
бесконечно-большом времени расчета получить абсолютно точное решение, говорят что это несмещенный рендерер (unbiased renderer). 
К таким программам, например, можно отнести все виды трассировщиков путей. 
Если программа использует метод Монте-Карло, но при этом в пределе (при бесконечно-большом времени расчета) не позволяет получить абсолютно точное решение, говорят, 
что такой рендерер является смещенным (например программы, использующие фотонные карты и/или кэш освещенности). 
К программам, не использующим метод Монте-Карло понятия смещенный/несмещенный неприменимы. 
Например, некорректно иcпользовать указанные термины по отношению к программе, использующей метод излучательности \cite{Goral:1984:MIL:800031.808601}. 

\noindent Смещенность алгоритма еще не означает, что он вычислет интеграл с низкой точностью. На практике смещение (bias) проявляется в виде цветных пятен. Если смещение небольшое, пятна незаметны для глаза. 
Как правило, смещенные методы на начальном этапе работы (при небольшом числе выборок или небольшом времени работы) дают более приемлимую для человечесткого глаза оценку, чем несмещенные. 
Однако в пределе, при длительном рассчете несмещенные методы обычно дают более качественное изображение.


\section{Стохастическая трассировка лучей}

\noindent Стохастическая трассировка лучей (также называемая Монте-Карло трассировкой) вычисляет значение интеграла освещенности (для фиксировнаного направления $\phi_{r}, \theta_{r}$) напрямую по 
формуле \ref{eq:montecarlo3}. Метод был предложен James-ом Kajiya в 1986 году \cite{Kajiya:1986:RE:15886.15902}.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{montecarlo1}           
  \caption{\label{fig:montecarlo1}
           Иллюстрация монте-карло трассировки лучей.}   
\end{figure}

\noindent Из виртуальной камеры испускается луч, который трассируется в сцену. В точке пересечения с поверхностью выпускается некоторое число лучей по полусфере и для каждого луча процедура выполняется рекурсивно (рис. \ref{fig:montecarlo1}). 
Полученные значения яркости на каждом уровне рекурсии складываются с учетом формулы \ref{eq:montecarlo3}.

\noindent В дальнейшем будем рассматривать модификацию стохастической трассировки лучей, при которой отраженный луч всегда один (рис \ref{fig:pathtracing1}). 
Такое упрощение немного замедляет сходимость метода, но позволяет упростить реализацию и при желании избавиться от рекурсии в реализации алгоритма. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{pathtracing1}           
  \caption{\label{fig:pathtracing1}
           Иллюстрация монте-карло трассировки путей.}   
\end{figure}

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

\begin{algorithm}[H]
\SetAlgoLined
\SetKwProg{Fn}{Function}{ is}{end}

\KwData{ray - Луч, depth - глубина рекурсии}
\KwResult{Монте-Карло выборка значения освещенности}

\Fn{TracePath(ray: Ray; depth : Integer) : Integer}{

\If{depth == MAX\_DEPTH}
  { return Black\;}

$ hit \leftarrow RaySceneIntersection(ray)\; $
  
\If{not hit.exists}
  { return Black\;}  
  
$ m \leftarrow hit.material \; $

\If{ m.IsLight}
  { return m.emittance\;}  
  
$ newRay.O \leftarrow hit.pos \; $
$ newRay.D \leftarrow RandomUnitVectorInHemisphereOf(hit.norm) \; $
$ cosTheta \leftarrow dot(newRay.D, hit.norm) \; $
$ BRDF \leftarrow m.reflectance * cosTheta \; $

return BRDF * TracePath(newRay, depth + 1) \; 
}
\caption{Простейший вид Монте-Карло трассировки путей. В данном псевдокоде emittance обозначает светимость и измеряется в ваттах на квадратный метр.}
\label{path_tracing_alg}
\end{algorithm}

\noindent Данный алгоритм также называют обратной трассировкой путей. 
Для генерации случайного вектора по полусфере необходимо реализовать функцию $RandomUnitVectorInHemisphereOf$, отображающую случайные числа $r_1$ и $r_2$ с равномерным распределением из интервала [0,1] на сферу.
Такое отображение можно реализовать использую формулы \ref{eq:to_hemisphere1} (результирующие координаты - $(x,y,z)$) \cite{Suffern:2007:RTG:1324795}.

\begin{eqnarray}\label{eq:to_hemisphere1}
\phi &=& 2 * \pi * r_1        \nonumber \\
h &=& 2*r_2 - 1               \nonumber \\
x &=& sin(\phi)*\sqrt{1-h^2}  \nonumber \\
y &=& cos(\phi)*\sqrt{1-h^2}  \nonumber \\
z &=& h 
\end{eqnarray}

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{pathtracing2}           
  \caption{\label{fig:pathtracing2}
           Иллюстрация алгоритма трассировки путей. Корректный перенос света моделируется путем просчета большого числа путей.}   
\end{figure}


\subsection{Прогрессивное вычисление интеграла}

\noindent Для того чтобы иметь возможность визуализироваь промежуточный результат в процессе расчета изображения, удобно использовать прогрессивную схему вычисления интеграла. 
Будем полагать, что мы используем простейший бокс-фильтр, усредняющий все пути для каждого пиксела. То есть цвет пиксела вычисляется по следующей формуле:

\begin{align}\label{eq:colsumm}
C = \sum_{i=1}^N \dfrac{P_{i}}{N}
\end{align}

Где:

\begin{enumerate}
\item $ P_{i} $ - значение Монте-Карло выборки на i-ом шаге.
\item $ N $     - число Монте-Карло выборок.
\end{enumerate}

\noindent Будем вычислять цвет пиксела $С$ в виде сходящейся к $C$ последовательности $ \{C_{i}\} $ Пусть далее $ C_{i} $ - цвет пиксела изображения на i-ом проходе. Тогда он может быть вычислен прогрессивно по следующей схеме:

\begin{align}\label{eq:progressive_scheme}
C_{0} &= P_{0} \nonumber \\
C_{1} &= \dfrac{1}{2}C_{0} + \dfrac{1}{2}P_{1} \nonumber \\ 
C_{2} &= \dfrac{2}{3}C_{1} + \dfrac{1}{3}P_{2} \nonumber \\
C_{3} &= \dfrac{3}{4}C_{2} + \dfrac{1}{4}P_{3} \nonumber \\  
&... \nonumber \\
C_{n} &= \dfrac{n}{n+1}C_{n-1} + \dfrac{1}{n+1}P_{n}
\end{align}

\noindent Данная схема позволяет пользователю наблюдать постепенно улучшающееся изображение (с постепенно снижающимся уровнем шума). 
Если пользователя устраивает полученное изображение на шаге i, он может остановить расчет.

\subsubsection{Критерий остановки}

\noindent Однако ручное управление невсегда возможно применить. Например, пользователь не может контролировать каждый отдельный пиксель изображения. 
Одним из наиболее стабильных критериев, позволяющих оценить точность полученного решения (и осуществить автоматическую остановку расчета) 
является дисперсия $D$ и получаемое на ее основе стандартное отклонение $\sigma$. Дисперсию и стандартное отклонение удобно вычислять по формуле \ref{eq:disp}.

\begin{eqnarray} \label{eq:disp}
D &=& \dfrac{ \sum_{i=1}^n X_i^2  - \dfrac{ (\sum_{i=1}^n X_i)^2 }{n} }{n-1} \nonumber \\
\sigma &=& \sqrt{ \dfrac{D}{n} }
\end{eqnarray}

\noindent Для реализации формулы \ref{eq:disp} достаточно накапливать прогрессивно (по схеме \ref{eq:progressive_scheme}) сумму цветов и сумму квадратов цветов.

\subsection{Применение выборки по значимости}

\noindent Известно, что для Ламбретовых поверностей наибольший вклад дают лучи близкие к направлению нормали, 
а наименьший - лучи почти перпендекулярные нормали. 
Это следует того, что значение падающей освещенности умножается на косинус между нормальню и направлением падающей освещенности ($cosTheta$ в алгоритме \ref{path_tracing_alg}). 

\noindent Чтобы применить выборку по значимости, необходимо генерировать такой вектор отражения, который соотвествует косинусоидальному 
распределению вокруг направления нормали (рис. \ref{fig:cosine_distribution}). 
Это можно сделать эффективно используя соотношения \ref{eq:to_hemisphere2}, которые позволяют отобразить 2 равномерно-распределенных числа 
на сферу с косинусоидальным распределением \cite{Suffern:2007:RTG:1324795}.
В соотношениниях \ref{eq:to_hemisphere2} переменная $e$ задает степень косинуса в косинусоидальном распределении.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{cosine_distribution}           
  \caption{\label{fig:cosine_distribution}
           Выборка по значимости для Ламбертовского отражения.}   
\end{figure}

\begin{eqnarray}\label{eq:to_hemisphere2}
e &=& 1                          \nonumber \\
\phi &=& 2 * \pi * r_1           \nonumber \\
cosTheta &=& (1-r_2)^{1/(e+1)}   \nonumber \\
sinTheta &=& \sqrt{1-cosTheta^2} \nonumber \\
x &=& sinTheta*cos(\phi)         \nonumber \\
y &=& sinTheta*sin(\phi)         \nonumber \\
z &=& cosTheta
\end{eqnarray}


\noindent Тогда для того, чтобы результат не изменился, необходимо делить значение BRDF (или приходящей яркости умноженной на BRDF) на плотность вероятности выбранного распределения, то есть $cosTheta$. 
В результате, при применении выборки по значимости для Ламбертовских материалов $cosTheta$ уходит, а функция $RandomUnitVectorInHemisphereOf$ 
заменяется на функцию $RandomCosineVectorOf$ и алгоритм \ref{path_tracing_alg} переходит в алгоритм \ref{path_tracing_alg2}. 

\begin{algorithm}[H]
\SetAlgoLined
\SetKwProg{Fn}{Function}{ is}{end}

\KwData{ray - Луч, depth - глубина рекурсии}
\KwResult{Монте-Карло выборка значения освещенности}

\Fn{TracePath(ray: Ray; depth : Integer) : Integer}{

\If{depth == MAX\_DEPTH}
  { return Black\;}

$ hit \leftarrow RaySceneIntersection(ray)\; $
  
\If{not hit.exists}
  { return Black\;}  
  
$ m \leftarrow hit.material \; $

\If{ m.IsLight}
  { return m.emittance\;}  
  
$ newRay.O \leftarrow hit.pos \; $
$ newRay.D \leftarrow RandomCosineVectorOf(hit.norm) \; $
$ BRDF \leftarrow m.reflectance\; $

return BRDF * TracePath(newRay, depth + 1) \; 
}
\caption{Трассировка путей с учетом косинусоидального распределения отраженных лучей.}
\label{path_tracing_alg2}
\end{algorithm}

\subsection{Учет сложных материалов}

\noindent Ранее был рассмотрен алгоритм трассировки путей (алгоритмы \ref{path_tracing_alg} и \ref{path_tracing_alg2}), и делалось предположение, что все материалы в сцене Ламбертовы. 
Такие материалы отражают свет равномерно во все стороны. Однако, большинство материалов реального мира обладают намного более сложными ДФО, 
а также могут не только отражать но и пропускать свет (рис. \ref{fig:dgs}). 

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{dgs}           
  \caption{\label{fig:dgs}
           Различные виды взаимодесйствий света с материалом. Diffuse(D) - Ламбертовское отражения. Specular(S) - зеркальное. Glossy(G) - матовое.}   
\end{figure}

\noindent Идеально-зеркальное(S) и матовое(G) отражения моделируются при помощи единого механизма (который будет рассмотрен ниже), 
поэтому в дальнейшем будем рассматривать взаимодействия G и S совместно и называть такой тип отражения зеркальным (S). 

\noindent Сложные ДФО часто представляются в виде композиции более простых компонент. 
Классической композицией является композиция ДФО Ламберта и Фонга (или любой другой ДФО, моделирующей зеркальное и/или матовое (glossy) отражение). 
В этом случае ДФО получается как линейная комбинация (с коэффициентами kd и ks) отдельных компонент. 
Помимо этого, для прозрачных объектов необходимо учитывать пропускание (по аналогии с ДФО вводят функцию ДПФ или BTDF) с коэффициентом kt. 
Для того чтобы учитывать такие композитные свойства материала, используют следующий механизм для выбора отраженного/преломленного луча:

\begin{eqnarray*}\label{eq:kdkskt}
\Sigma = kd + ks + kt \\
\xi = rnd(0,1) \\
\xi \in [0, \dfrac{kd}{\Sigma}] \Rightarrow diffuse \\
\xi \in [\dfrac{kd}{\Sigma}, \dfrac{kd+ks}{\Sigma}] \Rightarrow specular \\
\xi \in [\dfrac{kd+ks}{\Sigma}, 1] \Rightarrow transmittion \\
\end{eqnarray*}

\noindent kd - коэффициент диффузного отражения. ks - коэффииент зеркального отражения. kt - коэффициент пропускания. 
Важно отметить, что при использовании модели выбора компоненты ДФО указанной выше, полученное значение яркости 
не нужно умножать на значения коэффициентов (kd, ks или kt) т.к. эти коэфициенты получаются стохастически из самой схемы. 

\noindent Если необходимо учитывать цвет, схема усложняется:

\begin{eqnarray*}\label{eq:kdkskt2}
Pd = max(kd.r, kd.g, kd.b) \\ 
Ps = max(ks.r, ks.g, ks.b) \\
Pt = max(kt.r, kt.g, kt.b) \\
\Sigma = Pd + Ps + Pt \\
\xi = rnd(0,1) \\
\xi \in [0, \dfrac{Pd}{\Sigma}] \Rightarrow diffuse \\
\xi \in [\dfrac{Pd}{\Sigma}, \dfrac{Pd+Ps}{\Sigma}] \Rightarrow specular \\
\xi \in [\dfrac{Pd+Ps}{\Sigma}, 1] \Rightarrow transmittion \\
\end{eqnarray*}

\noindent Однако это еще не все. Теперь полученную яркость необходимо делить на значение плотности вероятности распределения выбранной компоненты. 
Например, если было выбрано зеркальное отражение, яркость луча необходимо пересчитать по следующей формуле:

\begin{eqnarray*}
color.r = color.r / (\dfrac{Ps}{\Sigma}) \\
color.g = color.g / (\dfrac{Ps}{\Sigma}) \\
color.b = color.b / (\dfrac{Ps}{\Sigma}) \\
\end{eqnarray*}

\noindent Для того чтобы применять выборку по значмости с учетом матовых (glossy) ДФО, необходимо генерировать лучи со степенным косинусоидальным 
распределением в направлении отраженного вектора. 
Такое распределение представляет собой лепесток (cosine lobe) вокруг вектора отражения \cite{Suffern:2007:RTG:1324795}.
Чтобы его реализовать, необходимо изменить значение переменной $e$ в соотношениях \ref{eq:to_hemisphere2} присвоив ему соответствующую степень. 



\noindent Результат работы алгоритма представлен на рисунке \ref{fig:cornell_box_stupid}. Алгоритм позволяет корректно вычислить интеграл освещенности и получить весь спектр эффектов геометричекой оптики (рис. \ref{fig:cornell_box_stupid}).

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{cornell_box_stupid}           
  \caption{\label{fig:cornell_box_stupid}
           Пример работы алгоритма трассировки путей на классической сцене Cornell Box.}   
\end{figure}

\subsection{Теневые лучи}

\noindent Представленный ранее алгоритм работает чрезвычайно долго, причем его скорость сильно зависит от размера источника освещения; 
          алгоритм сходитя тем дольше, чем меньше размер источника. 
		  Причина такой зависимости заключается в том, что вероятность случайного луча попасть в источник света зависит от его размера и 
		  расстояния до источника. Из этого следует, что трассировка путей в том виде, в котором она была рассмотрена ранее, не применима 
		  для сцен с точечными источниками. Чтобы устранить зависимость времени расчета от размера источника и добавить поддержку точечных 
		  источников света, вводят так называемые теневые лучи. 
		  Идея заключается в том, чтобы генерировать выборки в соответствие с формой ДФО, создавая случайные отраженные лучи по полусфере, 
		  но также пускать лучи напрямую к источникам света, (рис. \ref{fig:lightsampling}). 

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{lightsampling}           
  \caption{\label{fig:lightsampling}
           Сэмплирование (генерация выборок) по ДФО (справа) и по источнику света (слева). 
		   Первую стратегию генерации выборок принято называть неявной. Вторую - явной. }   
\end{figure}

\noindent При реализации теневых лучей следует иметь ввиду 2 важных момента: 
\begin{enumerate}% начало помеченного списка
\item Чтобы не учитывать источник света 2 раза, случайные лучи, полученные в результате сэмплирования ДФО при попадании в источник света должны терминироваться и возвращать черный цвет. 
\item Если для случайного луча вероятность попасть в источник света зависит от его площади поверхности и расстояния до него, то для теневых лучей вероятность попадания луча в источник света 
равна единице по определению (поскольку лучи всегда испускаются точно в источник). 
Для того чтобы корректно применить теневые лучи, необходимо вычислять такое же значение яркости теневого луча, как если бы луч попал в источник 
случайно, сэмплируя ДФО. 
Для этого вводится переменная $implicitP$ (алгоритм \ref{path_tracing_alg3}). Она отражает вероятность попадания в источник света при использовании неявной стратегии. 
$R$ - расстояние от точки, из которой испускается теневой луч до источника света; 
$light.area$ - площадь поверхности источника; 
emittance - ранее использованная светимость поверхности источника света (измеряется в ваттах на квадратный метр).

\end{enumerate}% конец помеченного списка


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{shadowrays}           
  \caption{\label{fig:shadowrays}
           Теневые лучи обозначены пунктиром. Отраженный луч, попавший в источник зачеркнут пеерекрестной красной линией. 
		   При возникновении такие лучи терминируются с тем чтобы не учитывать источник 2 раза.}   
\end{figure}


\begin{algorithm}[H]
\SetAlgoLined
\SetKwProg{Fn}{Function}{ is}{end}

\KwData{ray - Луч, depth - глубина рекурсии}
\KwResult{Монте-Карло выборка значения освещенности}

\Fn{TracePath(ray: Ray; depth : Integer) : Integer}{

\If{depth == MAX\_DEPTH}
  { return Black\;}

$ hit \leftarrow RaySceneIntersection(ray)\; $
  
\If{not hit.exists}
  { return Black\;}  
  
$ m \leftarrow hit.material \; $

\If{ m.IsLight}
  { return Black\;}  
  
$ lpos   \leftarrow LightSample(light) \; $

$ shadow \leftarrow Visibility(hit.pos, lpos) \; $  

$ R      \leftarrow dist(hit.pos, lpos) \; $  

$ sdir       \leftarrow normalize(lpos - hit.pos) \; $ 
$ cosTheta1  \leftarrow -dot(sdir, light.norm)    \; $ 
$ cosTheta2  \leftarrow dot(sdir, hit.norm)       \; $ 

$ implicitP  \leftarrow light.area*cosTheta1*cosTheta2/(\pi*R^2) \; $ 
  
  
$ newRay.O \leftarrow hit.pos \; $
$ newRay.D \leftarrow RandomCosineVectorOf(hit.norm) \; $
$ BRDF    \leftarrow m.reflectance\; $

$ explicitColor \leftarrow shadow*BRDF*light.emittance*implicitP  $

return explicitColor + BRDF * TracePath(newRay, depth + 1) \; 
}
\caption{Трассировка путей с теневыми лучами.}
\label{path_tracing_alg3}
\end{algorithm}

\noindent Коэффициент $implicitP$ получается при переходе от рассмотренной ранее сферической формы уравнения освещенности (уравнение \ref{eq:integral}) к площадной форме этого уравнения (рис \ref{fig:areaequation}, соотношение \ref{eq:areaintegral}).

\begin{equation}\label{eq:areaintegral}
I(x,\psi) = \int_ {\psi^{'}}{} R(x,\psi^{'})L(x, \psi^{'}) \dfrac{dA^{'}cos\theta cos\theta^{'}}{||x - x^{'}||}
\end{equation}


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{areaequation}           
  \caption{\label{fig:areaequation}
           Иллюстрация к площадной  форме интеграла освещенности \cite{Shirley92distributionray}.}   
\end{figure}


\noindent На рисунке \ref{fig:shadowrays} изображена ситуация, при которой случайный отраженный луч попадает в источник света. Такой луч необходимо терминировать и вернуть 0 в качестве яркости от источника света.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{cornell_box_shadow}           
  \caption{\label{fig:cornell_box_shadow}
           Пример работы алгоритма трассировки путей c теневыми лучами (алгоритм \ref{path_tracing_alg3}) на классической сцене Cornell Box. }   
\end{figure}

\noindent Реализация трассировки путей с теневыми лучами значительно быстрее сходится (рис. \ref{fig:comparisonshadow2}) и позволяет 
обрабатывать источники любых размеров. Однако, как нетрудно заметить (рис. \ref{fig:cornell_box_shadow}), источник света стал черным и на изображении 
исчезают каустики. Первую проблему (черный источник света) легко исправить если возвращать цвет источника вместо черного при отсутствии диффузных и матовых (glossy) переотражений. 
Это всегда можно делать поскольку теневые лучи применяются только для учета матовой и диффузной компонент ДФО.

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{comparisonshadow2}           
  \caption{\label{fig:comparisonshadow2}
           Сравнение простого алгоритма трассировки путей (64 путя на пиксел, слева на изображении) и алгоритма трассировки путей 
		   использующего теневые лучи (32 путя на пиксел, справа на изображении). Оба изображения получены за примерно одинаковое время.}   
\end{figure}

\noindent Вторая проблема - отсутствие кауcтик. Чтобы вернуть каустики, определим для начала в каких условиях они возникают. 
Каустик - яркое пятно на диффузной поверхности, вызванное непрямым освещением посредством зеркальных отражений (или преломлений) света, излученного из источника. 
Таким образом, если в процессе трассироки путей после диффузного переотражения встретилось зеркальное, на месте упомянутого диффузного отражения может возникнуть каустик. 
Добавить каустики в трассировку путей с теневыми лучами, таким образом, достаточно просто:


\begin{enumerate}
\item Трассируем путь из глаза с применением теневых лучей до тех пор пока не встретилось зеркальное переотражение после диффузного.
\item Как только указанное условие выполнилось, превращаем оставшуюся часть пути в 'каустический путь'. 
     Такой путь не использует теневые лучи и строит оставшуюся его часть руководствуется только формой ДФО. 
\end{enumerate}
\noindent В результате получем изображение на рисунке \ref{fig:cornell_box_stupid}.

\noindent Рисунок \ref{fig:comparisonshadow2} показывает преимущество трассировки путей с теневыми лучами над простой трассировкой путей, 
рассмотренной ранее. Поскольку трассировка путей без теневых лучей обрабатывает пути быстрее (в описанной выше реализации на основе алгоритмов \ref{path_tracing_alg2} и \ref{path_tracing_alg3} примерно в 2 раза), 
корректно сравнивать 64 путя на пиксел для трассировки путей без теневых лучей и 32 путя на пиксел с теневыми лучами.


\subsection{Две стратегии сэмплирования}

\noindent Ранее были рассмотрены стратегими сэмплирования (стратегиями генерации выборок) - сэмплирование по ДФО и сэмплирование по источнику света (рис. \ref{fig:lightsampling}). Будем называть сэмплирование по ДФО \underline{неявной} стратегией, а сэмплироание по источнику -  \underline{явной}.

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

\subsection{Многократная выборка по значимости}

Для того чтобы корректно комбинировать явную и неявную стратегии, Эриком Вичем был предложен метод многократной выборки по значимости (англ. Multiple Importance Sampling, MIS) \cite{Veach:1998:RMC:927297}, который позволяет скомбинировать результаты нескольких способов выборки (нескольких стратегий сэмплирования) в одну несмещенную оценку. На практике реализация многократной выборки по значимости имеет некоторые неочевидные, но существенные детали. Мы рассмотрим их далее.


Пусть мы имеем сцену с ламбертовскими материалами и площадным источника света (c площадью $A$) с равномерным распределением. Пусть мы находимся в некоторой точке $X$ на поверхности. И пусть для явной стратегии мы выбираем на источники точку $L_{pos}$. Пусть $L_{e}$ - яркость пришедшая от источника света в направлении теневого луча, а $ L_{i} $ яркость в точке источника, куда попал случайно отраженный луч.  

Тогда для ламбертовского материала плотность вероятности случайной величины, распределенной косинусоидально, будет $p_{i} = cos(\theta) / \pi$ где  $ \theta $ - угол между нормалью и отраженным лучом. А плотность вероятности случайной величины при сэмплировании источника cвета будет $p_{e} = \frac{d^2}{cos(\phi)*A}$. Где $d$ - расстояние от точки на поверхности до точки пересечения луча и поверхности источника лучом при сэмплировании явной стратегией ($d = dist(L_{pos},X)$), а $cos(\phi)$ - угол между нормалью источника и направлением луча (при сэмплировании неявной стратегией) умноженного на $-1$. Следует отметить, что поскольку мы выбираем на источнике точку $L_{pos}$ равномерно, изначальная плотность вероятности при сэмплировании площадного источника равна $ \frac{1}{A}$. Множитель $\frac{d^2}{cos(\phi)} $ переводит плотность вероятности в радиальную систему координат - в ту же самую систему, в которой задана плотность вероятности $p_{i} = cos(\theta) / \pi$. Итак, имеем:

\begin{eqnarray*}
p_{i} &= cos(\theta) / \pi \\
p_{e} &= \dfrac{d^2}{cos(\phi)*A} \\
w_{e} &= \dfrac{p^{\beta}_{e}}{p^{\beta}_{i} + p^{\beta}_{e}} \\
w_{i} &= \dfrac{p^{\beta}_{i}}{p^{\beta}_{i} + p^{\beta}_{e}} \\
\beta &= 2
\end{eqnarray*}


$ w_{i} $ - вес для учета яркости полученной по неявной стратегии. $ w_{e} $ - вес для учета яркости полученной по явной стратегии. Схема многократной выборки по значимости записывается по формуле \ref{eq:mis_pt_diffuse}.

\begin{eqnarray}\label{eq:mis_pt_diffuse}
  L_{k} =  w_{e}\dfrac{L_{e}}{ p_{e}} + w_{i}\dfrac{L_{i}}{p_{i}}
\end{eqnarray}

\noindent Оптимальность такого взвешивания показана в \cite{Veach:1998:RMC:927297}.

\subsubsection{Детали реализации}

Казалось бы, все переменные известны и мы готовы к применению многократной выборки по значимости. Достаточно лишь хранить раздельно яркости $L_{e}$, $L_{i}$ и соответствующие плотности вероятности. Если случайно-выстреленный луч попал в источник, используем многократную выборку по значимости. Если не попал - не используем и просто складываем  $L_{e}$ и $L_{i}$. Однако, это приведет к \underline{неправильному} результату. Причина этого заключается в невыполнении некоторого важного условия MIS. Сумма весов для обоих стратегий должна быть \underline{стохастически} равна 1. 


Рассмотрим, что это означает. Пусть мы имеем серию выборок (или сэмплов, экспериментов) $f(x_1)$ с плотностью вероятности случайной величины $x_1$ равной $p_1$. Пусть также мы имеем другую серию экспериментов $f(x_2)$ с плотностью вероятности случайной величины $x_2$ равной $p_2$. Имея две серии экспериментов мы можем оценивать значения интеграла при помощи метода Монте-Карло (уравнение \ref{eq:mis_example}):

\begin{eqnarray}\label{eq:mis_example}
  I_{1} = \frac{1}{N1} \sum_{k=0}^{N1}{ \dfrac{f(x_{1,k})}{p_{1,k}} } \\
  I_{2} = \frac{1}{N2} \sum_{k=0}^{N2}{ \dfrac{f(x_{2,k})}{p_{2,k}} }
\end{eqnarray}
		  
Поскольку обе оценки в пределе сходятся к одному и тому же значению, имеем право смешивать их любыми коэффициентами, сумма которых равна 1 - например так $I_{new} = 0.25*I_1 + 0.75*I_2$. Обратите внимание, что весовые коэффициенты вышли за знак суммы. Это и есть условие, обсуждаемое ранее - сумма весов для обоих стратегий должна быть стохастически равна 1. То есть:

\begin{eqnarray}\label{eq:mis_wsumm}
  \frac{1}{N1}\sum_{k=0}^{N1}{ w_{i,k} } + \frac{1}{N2}\sum_{k=0}^{N2}{ w_{e,k} } = 1  \\
\end{eqnarray}

При этом следует отметить, что сумма весов на каждой итерации $k$ вовсе не обязана быть равна единицы, более того - число итераций для разных стратегий может быть разным ($N1$ и $N2$). Теперь рассмотрим рисунок \ref{fig:mis_sampling_sr} и алгоритм, правильно вычисляющий веса $w_i$ и $w_e$.


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{mis_sampling_sr}           
  \caption{\label{fig:mis_sampling_sr}
           Сэмплирование источника двумя различными стратегиями. Теневой луч (синий, суффикс s) отвечает за явную стратегию. Случайно-отраженный (зеленый, суффикс r) - за неявную стратегию. }   
\end{figure}


Для того чтобы правильно вычислить вес $w_i$ вам необходимо знать не только плотность вероятности явной стратегии для данной точки ($L_{pos [s]}$), но и плотность вероятности для неявной стратегии - как если бы вы попали в точноcти в ту же самую точку ($L_{pos [s]}$). Аналогично и для неявной стратегии. Когда вы попадаете в источник света, необходимо вычислять плотность вероятности попадания в конкретную точку $L_{pos [r]}$ обеими стратегиями. Сопоставим обозначения на рисунке \ref{fig:mis_sampling_sr} и в алгоритме \ref{path_tracing_alg4}:

\begin{itemize}
\item $D$ - $d_r$
\item $R$ - $d_s$
\item $cosTheta1$ - косинус угла $\phi_{s}$
\item $dot(light.norm, -ray.direction)$ - косинус угла $\phi_{r}$
\item $lpos$ - $L_{pos [s]}$
\item $prevHit$ - позиция и нормаль к поверхности в точке $X$.
\end{itemize}


\begin{algorithm}[H]
\SetAlgoLined
\SetKwProg{Fn}{Function}{ is}{end}

\KwData{ray - Луч, depth - глубина рекурсии}
\KwResult{Монте-Карло выборка значения освещенности}

\Fn{TracePath(ray: Ray; depth : Integer; prevHit : Hit) : Integer}{

$ ...  \; $

$ ...  \; $

$ m \leftarrow hit.material \; $

\If{ m.IsLight}
 {
     $ D     \leftarrow dist(prevHit.position, hit.position) \; $
     
     $ p_{e} \leftarrow D^2/(light.surfaceArea*dot(light.norm, -ray.direction))  \; $ 
     
     $ p_{i} \leftarrow  dot(prevHit.norm, ray.direction)/\pi \; $ 
     
     $ w_{i} \leftarrow  p_{i}^2/(p_{e}^2 + p_{i}^2) \; $ 
    
    return light.emittance*$w_{i}$ \;
    
 }  
  
$ lpos   \leftarrow LightSample(light) \; $

$ shadow \leftarrow Visibility(hit.pos, lpos) \; $  

$ R      \leftarrow dist(hit.pos, lpos) \; $  

$ sdir       \leftarrow normalize(lpos - hit.pos) \; $ 
$ cosTheta1  \leftarrow -dot(sdir, light.norm)    \; $ 
$ cosTheta2  \leftarrow dot(sdir, hit.norm)       \; $ 

$ p_{e} \leftarrow R^2/(light.surfaceArea*cosTheta1)  \; $ 

$ p_{i} \leftarrow  cosTheta2/\pi \; $ 

$ w_{e}   \leftarrow  p_{e}^2/(p_{e}^2 + p_{i}^2) \; $ 
    
$ ...  \; $


return explicitColor*$w_{e}$ + BRDF*TracePath(newRay,depth+1,hit) \; 
}
\caption{Трассировка путей с MIS.}
\label{path_tracing_alg4}
\end{algorithm}


Хорошим примером реализации техники многократной выборки по значимости может служить соттветствующи пример в \cite{OptiX}. Ситуация, в которой многократная выборка по значимости действительно дает существенный выигрыш изображена на рис \ref{fig:misexample1}.


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{misexample1}           
  \caption{\label{fig:misexample1}
           Сравнение различных стратегий генерации выборки для сцены с матовым (glossy) отражением. 
		   Явная стратегия дает лучший результат, если источник находится далеко от повехности или имеет малый размер. 
		   Неявная стратегия дает лучший результат, если источник имеет большой размер или находится близко к поверхности. 
		   Многократная выборка по значимости использует преимущества обоих пододов. }   
\end{figure}

\subsection{Русская рулетка и глубина трассировки}

\noindent Одной из причин низкой эффективности трассировки путей путей может являться завышенное значение максимальной глубины трассировки. 
          Стеклянные объекты могут требовать значительного числа переотражений для корректной визуализации. 
		  В то же время, лучи попавшие на диффузные объекты, как правило, уже после нескольких переотражений теряют значимость. 
		  Одним из простых решений является терменирование пути при достижении определенного числа диффузных отскоков и/или низкого 
		  'коэффициента пропускания' (произведение коэффициентов отражения во всех точках пути). 
		  Параметр терменирования при этом требует ручной настройки. 
		  Однако, такой метод вносит смещение в получаемое решение и может дать неверный результат на сцене с очень яркими источниками и 
		  сложным перераспределением световой энергии \cite{Bogolepov:2013:phdthesis}.

\noindent Одно из решений, позволяющих автоматически вбирать глубину терминирования пути - стохастическое терминирования пути на освнове механизма русской рулетки (алгоритм \ref{path_tracing_alg_rr}). 
          Фактически, русская рулетка позволяет трассировать меньше путей на большой глубине. Однако, она присваивает им большую значимость, за счет чего в пределе получается верный результат. 


\begin{algorithm}[H]
\SetAlgoLined
\SetKwProg{Fn}{Function}{ is}{end}

\KwData{ray - Луч, depth - глубина рекурсии}
\KwResult{Монте-Карло выборка значения освещенности}

\Fn{TracePath(ray: Ray; depth : Integer) : Integer}{

$ pabsorb \leftarrow  terminationProbability(depth, ...)\; $
  
\If{random(0.0, 1.0) < pabsorb}
  { return Black\;}  
  
...\;

return result*(1/(1-pabsorb)) \; 
}
\caption{Трассировка путей с русской рулеткой.}
\label{path_tracing_alg_rr}
\end{algorithm}
		  
		  

\noindent Механихм русской рулетки, вообще говоря, увеличивает дисперсию в оцениваемом решении \cite{SzecsiSK03}, поэтому должен 
          применяться осторожно - только тогда когда ожидаемый вклад от текущего пути достаточно мал. В этом случае:

\begin{enumerate}
\item Дисперсия повысится незначительно, так как в среднем вклад от путей, для которых применяется русская рулетка относительно мал.
\item По сравнению с простым терминированием время трассировки не увеличится, поскольку число лучей трассируемых на большую глубину невелико.
\item Решение станет более точным за счет учета большего числа переотражений.
\end{enumerate}		  
		  
		  

\subsection{Резюме по обратной трассировке путей}

\noindent Рассмотренный алгоритм обратной трассировки путей с применением теневых лучей и многократной выборки по значимости является несмещенным и достаточно эффективным методом для вычисления интеграла освещенности; за исключением каустик, которые вычисляются этим алгоритмом довольно медленно. Для эффективного вычисления каустик далее будет рассмотрен алгоритмы прямой трассировки путей (Light Tracing) и алгоритм фотонных карт.  

\subsection{Грамматика путей}

\noindent Для классификации различных ситуаций, возникающих в процессе трассировки путей часто используют понятие грамматики путей \cite{Heckbert:1990:ART:97880.97895}. Этот способ классифицирует пути при помощи строк и регулярных выражений. Каждый символ в строке обозначает определенное событие, произошеднее с лучом (путем) в процессе трассировки.

\begin{enumerate}
\item S - (Specular); зеркальное отражение/преломление
\item D - (Diffuse); диффузное отражение
\item G - (Glossy); матовое отражение/преломление
\item V - (Volume); объемное рассеивание
\item L - (Light); источник света 
\item E - (Eye); глаз
\end{enumerate}

\noindent При этом под символами S и G понимается не только отражение но и преломление света. Грамматика удобна не только для описания различных ситуаций, возникающих в трассировке путей, но она позволяет также классифицировать видимые эффекты, однозначно ставя им в соответствие классы путей которые эти эффекты вызывают (для упрощения изложения атериала объемное рассеивание V мы далее не рассматриваем). Например, $EDL$ - прямое освещение диффузной поверхности. $E(S|G)^{+}L$ - яркий блик. $EDSL$ - каустик, видимый из камеры напрямую, обусловленный одним зеркальным переотражением света (например солнечный зайчик от зеркала). $EDS^{+}L$ - каустик, видимый из камеры напрямую, вызванный одни или более зеркальным переотражением. $ES^{+}DS^{+}L$ - каустик, видимый через стекло или зеркало, вызванный одним или более зеркальным переотражением.

\noindent При помощи грамматики путей удобно классифицировать алгоритмы по типу освещения, которые они способны расчитывать. Например, Трассировка лучей Уиттеда, рассмотренная нами в самом начале, строит пути $ES^{*}L$. Обратная трассировка путей строит пути вида $ E(S|D|G)^{*}L $.

\section{Прямая трассировка (Light Tracing)}

\noindent В противоположность обратной трассировки путей, рассмотренной выше, алгоритм прямой трассировки путей начинает свою работу из источника света и строит пути вида $ L(S|D|G)^{*}E $. При каждом соударении с поверхностью создается аналог теневого луча, направленный в камеру. Если луч успешно доходит до камеры, значение яркости записывается в пиксель, соответствующий проекции начала луча на экранную плоскость. Для реализации операции проекции вы можете использовать матрицу $ mProj $, которая была рассмотрена нами в разделе "Генерация луча".

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{emitphotons}           
  \caption{\label{fig:emitlightrays}
           Иллюстрация различных видов распределения световой энергии для различных источников света.}   
\end{figure}

\noindent Пути создаваемые на источнике должны быть сгенерированы в соответствии с распределением световой энергии для источника света (рис. \ref{fig:emitlightrays}). Обратите внимание, что если каждая точка поверхности излучает свет равномерно во все стороны (например площадный источник из рис. \ref{fig:cornell_box_stupid}), то исходя из геометрических соображений и применяя выборку по значимости, необходимо использовать косинусоидальное распределение для расчета направления пути точно так же, 
         как мы вычисляли ранее диффузное отражение. В обратной трассировке путей данный косинус получался стохастически, за счет того что число лучей попавших в источник зависит от направления, с которого они были выпущены.

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

\noindent Метод, реализующий эту идею существует и называется 'усеченная двунаправленная трассировка путей' \cite{Bogolepov:2013:phdthesis}. 
          Идея этого метода в упрощенном представлении изложена ниже: 

\begin{enumerate}
\item Будем производить обратную трассировку путей с теневыми лучами. При этом отключим расчет каустиков. Результат сохраняем в изображение с номером 1. На изображении 1 мы полчим картинку без каустиков.
\item Будем производить прямую трассировку из источника света, но на этот раз будем сохранять вклад только от тех путей, для которых первое переотражение было зеркальным. Результат сохраняем в изображение с номером 2. На изображении 2 мы получим только каустики. 
\item Финальное изображение можно получить простым сложением изображений 1 и 2.
\end{enumerate}

\noindent Усеченная двунаправленная трассировка путей дает почти корректный результат, 
          поскольку изображения 1 и 2 будут нести в себе освещение, обусловленное различными типами переотражений. 		  
		  Слово 'почти' означает, что при помощи такого метода не удасться эффективно рассчитать каустики, не видимые камерой напрямую. 
		  То есть не удасться эффективно вычислять освещение, вызванное путями вида $ES^{+}DS^{+}L$.
          Время расчета таких эфектов будет сопостовимо со временем их рассчета обыкновенной трассировкой путей. 		  
		  




\section{Двунаправленная трассировка путей}

\noindent Идея объеденить преимущества прямой и обратной трассировки привела к созданию алгоритма двунправленной трассировки путей (Bidirectaional Path Tracing, BDPT) \cite{Veach:1998:RMC:927297}. Основная идея двунаправленной трассировки путей заключается в том, чтобы соединять теневыми лучами не только источник с точками путя (как в обратной трассировки) или камеру с точками путя (как в прямой трассировки), но и различные точки путей друг с другом (рис. \ref{fig:bdpt1}). 

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{bdpt1}           
  \caption{\label{fig:bdpt1}
           Идея двунаправленной трассировки путей.}   
\end{figure}
\noindent Двунаправленная трассировка путей обладает лучшей сходимостью (по сравнению с прямой или обратной трассирокой) на сценах со сложными условиями освещения (преобладающее вторичное освещение) только при реализации многократной выборки по значимости. В противном случае, изображение получается сильно зашумленным. Однако даже при реализации многократной выборки по значимости, двунаправленная трассировка путей по прежнему не позволяет эффективно рассчитывать пути вида $ES^{+}DS^{+}L$ (каустики, видимые через зеркало или стекло). Такие пути могут быть эффективно вычислены при помощи фотонных карт, переноса света Метрополиса и метода соединения вершин (vertex merging)\cite{Georgiev:2011:BLT:2077378.2077412}.

\noindent Поскольку теневые лучи в двунаправленной трасировки могут соединять произвольные точки различных путей, реализация многократной выборки по значимости значительно усложняется. Если ранее для вычисления весов достаточно было рассмотреть только 2 стратегии, то теперь число рассматриваемых вариантов значительно увеличивается. Рассмотрим путь на рисунке \ref{fig:bdpt1}.

\subsection{Многократная выборка по значимости в BDPT}

\noindent Расмотрим некоторый путь $ x_{0} .. x_{5} $. Для такой постедовательности точек существует несколько способов соединения их теневым лучом (обозначено вертикальной чертой):

\begin{eqnarray*}
x_{0} &|& x_{1} \: x_{2} \: x_{3} \: x_{4} \: x_{5} \\
x_{0} \: x_{1} &|& x_{2} \: x_{3} \: x_{4} \: x_{5} \\
x_{0} \: x_{1} \: x_{2} &|& x_{3} \: x_{4} \: x_{5} \\
x_{0} \: x_{1} \: x_{2} \: x_{3} &|& x_{4} \: x_{5} \\
x_{0} \: x_{1} \: x_{2} \: x_{3} \: x_{4} &|& x_{5} 
\end{eqnarray*}

\noindent При вычислении весов для какого-либо из случаев, многократная выборка по значимости должна учитывать все остальные случаи. Формула \ref{eq:misbdpt} применяется для вычисления веса конкретного пути. 

\begin{align}\label{eq:misbdpt}
 w_i = \dfrac{p_i^\beta}{  \sum_{j=1}^N p_j^\beta }
\end{align}

\noindent Например, для случая $ x_{0} \: x_{1} \: x_{2} | x_{3} \: x_{4} \: x_{5} $ i равно 2; N равняется 5. $ p_i $ - плотность вероятности распределения отраженного луча для точки с номером 2. Таким образом, для указанного случая из 6 точек $ x_{0} .. x_{5} $ необходимо просчитать всевозможные пути (всего 5) и яркость k-ого сэмпла может быть вычислена по формуле \ref{eq:mis_bdpt} (N = 5):

\begin{eqnarray}\label{eq:mis_bdpt}
  L_{k} =  \sum_{j=1}^N w_j \dfrac{L_j}{p_j}
\end{eqnarray}

\section{Перенос света Метрополиса (MLT)}

\noindent Трассировка путей обладает достаточно низкой скоростью на сценах со сложными условиями освещения. Сценами со сложными условиями освещения будем называть 2 вида сцен:

\begin{enumerate}
\item Сцены, на которых преобладает вторичное освещение вызванное узкими и яркими световыми пятнами (рис. \ref{fig:complexlightscene}).
\item Сцены, на которых $ES^{+}DS^{+}L$ пути вносят значительный вклад (рис. \ref{fig:cornell_water}). 
\end{enumerate}

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{complexlightscene}           
  \caption{\label{fig:complexlightscene}
           Пример сцены со сложными условиями освещения.}   
\end{figure}

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{cornell_water}           
  \caption{\label{fig:cornell_water}
           Пример сцены с $ES^{+}DS^{+}L$ путями. Каустики, под водой (рис. \ref{fig:complexlightscene}).}   
\end{figure}

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

\noindent Для улучшения скорости сходимости на таких сценах Э. Вичем и Л. Гибсом был разработан алгоритм переноса света Метрополиса (Metropolis Light Transport, MLT) \cite{Veach:1997:MLT:258734.258775}. Идея этого метода аналогична применению выборки по значимости к под-интегральному выражению \ref{eq:subint}, однако имеет совершенно иную природу. 

\begin{equation}\label{eq:subint}
L(\phi_{i},\theta_{i})R(\phi_{i},\theta_{i},\phi_{r},\theta_{r})cos(n,l_{\phi_{i},\theta_{i}})
\end{equation}  

\noindent Рассмотрим Метод Монте-Карло в общем виде (формула \ref{eq:montecarlo5}). 

\begin{eqnarray}\label{eq:montecarlo5}
\int_a^b f(x)dx \approx \dfrac{1}{N} \sum_{i=1}^N \dfrac{f(u_{i})}{p(u_{i})}
\end{eqnarray}

\noindent Если бы мы имели функцию распределения $ p(x)$ пропорциональную $f(x)$, тогда для некоторой константы $c$ можно записать $ p(x) = cf(x) $. Подставив константу в выражение для интеграла \ref{eq:montecarlo5}, получим оценку для константы: 

\begin{equation}\label{eq:mltconst}
c = \dfrac{1}{\int_a^b f(x)dx}
\end{equation}

\noindent Таким образом, мы имеем возможность оценивать константу $c$ взяв, например, выборки с равномерным распределение и применив метод монте-карло. Алгоритм Метрополиса отвечает за оставшуюся часть задачи - создание $p(x)$ пропорционально $f(x)$. Это делается следующим образом:

\begin{enumerate}
\item Пусть $X$ - начальный вектор случайных параметров, которые были использованы при вычислении $P_i$. 
\item Сгенерируем путь $P_{i+1}$ путем мутаций случайных параметров $X$. Т.e. $Y = mutate(X)$.
\item Вычислим некоторым специальным образом вероятность мутации $ T(X|Y) = p(X \rightarrow Y) $ и вероятность обратной мутации $ T(Y|X) = p(Y \rightarrow X) $.  
\item Вычислим вероятность приятия мутации $X \rightarrow Y$ по формуле \ref{eq:mltaccept}.  

\end{enumerate}

\begin{eqnarray}\label{eq:mltaccept}
  a(Y|X) = min\{1,\dfrac{f(y)T(X|Y)}{f(x)T(Y|X)}\}
\end{eqnarray}

\noindent Обоснование того почему выше-описанный алгоритм работает можно найти в \cite{Veach:1997:MLT:258734.258775}. 
Важным свойством алгоритма Метрополиса является то, что он строит $p(x)$ даже если значения $f(x)$ могут быть вычислены только стохастически, 
при помощи выборок. Это именно та ситуация, которая присутствует в трассировке путей. 
Наиболее сложная часть в алгоритме переноса света Метрополиса - выбор хорошей стратегии мутации и вычисление вероятности мутаций. 
За подробностями реализации алгоритма MLT рекомендуется обратиться к \cite{Cline05apractical} и \cite{Veach:1997:MLT:258734.258775}.



\section{Фотонные карты}

\noindent Рассмотренные ранее несмещенные методы на основе трассировки путей работают в терминах яркости. Каждый путь в них отвечал за перенос яркости от источника света к камере. Метод фотонных карт работает в терминах потока. Идея этого метода заключается в том, чтобы вычислить распределение световой энергии по сцене при помощи трассировки множества частиц, переносящих порцию световой энергии. После чего можно оценивать интеграл освещенности, выполняя на поверхности поиск ближайших фотонов. Рассмотрим интеграл освещенности (формула \ref{eq:ph_integral}). Выразим функцию падающей освещенности в терминах потока. Подставим $ L(\phi_{i},\theta_{i}) $ в интеграл и получим: 

\begin{eqnarray}\label{eq:ph_integral}
I(\phi_{r}, \theta_{r}) &=& \iint\limits_{\phi_{i}\theta_{i}} L(\phi_{i},\theta_{i})R(\phi_{i},\theta_{i},\phi_{r},\theta_{r})cos(n,l_{\phi_{i},\theta_{i}}) \mathrm{d}\phi_{i}  \mathrm{d}\theta_{i} \\
L(\phi_{i},\theta_{i}) &=& \dfrac{d^{2}\Phi}{cos(n,l_{\phi_{i},\theta_{i}})d\phi_{i}d\theta_{i}dA } \\
\end{eqnarray}

\begin{eqnarray}\label{eq:ph_integral2}
I(\phi_{r}, \theta_{r}) = \iint\limits_{\phi_{i}\theta_{i}} \dfrac{d^{2}\Phi}{cos(n,l_{\phi_{i},\theta_{i}})d\phi_{i}d\theta_{i}dA } R(\phi_{i},\theta_{i},\phi_{r},\theta_{r})cos(n,l_{\phi_{i},\theta_{i}}) \mathrm{d}\phi_{i}  \mathrm{d}\theta_{i} \\
I(\phi_{r}, \theta_{r}) = \iint\limits_{\phi_{i}\theta_{i}} \dfrac{d^{2}\Phi}{dA } R(\phi_{i},\theta_{i},\phi_{r},\theta_{r})
\end{eqnarray} 

\noindent Дифференциальный поток можно аппроксимировать суммой энергий фотонов в окрестности. Таким образом, в результате преобразований получаем сумму K ближайших фотонов. 

\begin{eqnarray}\label{eq:ph_integral3}
I(\phi_{r}, \theta_{r}) = \iint\limits_{\phi_{i}\theta_{i}} \dfrac{d^{2}\Phi}{dA } R(\phi_{i},\theta_{i},\phi_{r},\theta_{r}) \\
I(\phi_{r}, \theta_{r}) \approx  \sum_{i=1}^K R(\phi_{i},\theta_{i},\phi_{r},\theta_{r}) \dfrac{\Delta\Phi(\phi_{i},\theta_{i})}{\Delta A}  =  \dfrac{1}{\pi r^{2}} \sum_{i=1}^K R(\phi_{i},\theta_{i},\phi_{r},\theta_{r}) \Delta \Phi(\phi_{i},\theta_{i})
\end{eqnarray}

\noindent В выражении \ref{eq:ph_integral3} переменная $r$ отвечает за радиус диска на поверхности, в пределах которого выполняется поиск ближайших фотонов. Метод фотонных карт, таким образом, состоит из 4 шагов:

\begin{enumerate}
\item Испускание фотонов из источников света
\item Трассировка фотонов
\item Построение фотонной карты
\item Сбор освещенности
\end{enumerate}

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{photontrace}           
  \caption{\label{fig:photontrace}
           Процесс трассировки и сохранения фотонов на поверхностях.}   
\end{figure}

\paragraph{Испускание фотонов.} Фотоны в данном методе – это частицы, переносящие некоторую небольшую порцию световой энергии. 
На начальном этапе фотоны испускаются из источника света в соответствие с распределением световой энергии для данного источника. 
Например известно, что точечный или сферический источник света (такой как солнце) испускают свет равномерно во всех направлениях. 
Площадные источники света имеют косинусоидальное распределение, имеющее максимум по направлению, совпадающему с нормалью к плоскости источника. 
Стадия испускания фотонов не отличается от стадии испускания путей вида $ L(S|D|G)^{*}E $ в прямой трассировки путей. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{emitphotons}           
  \caption{\label{fig:emitphotons}
           Иллюстрация различных видов распределения световой энергии для различных источников света.}   
\end{figure}

\paragraph{Трассировка фотонов. Русская рулетка.} Трассировка фотонов происходит похожим на прямую трассировку путей образом за исключением одного важного отличия - механизма русской рулетки (этот механизм можно также применять и в трассировке путей \cite{Bogolepov:2013:phdthesis}, но исторически он появился именно в методе фотонны карт) \cite{Jensen:2007:HQR:1281500.1281593}. 
Легче всего продемонстрировать механизм русской рулетки на примере. Допустим мы имеем ламбертовский материал с коэффициентом отражения 0.25 и фотон падающий на эту поверхность с энергиейм равной 1. Очевидным решением является модифицировать энергию фотона, умножив ее на 0.25 и продолжить трассировку фотона. Механизм русской рулетки поступает по другому. Вместо того чтобы модифицировать энергию фотона, его энергия сохраняется, но с вероятностью 0.25 фотон отражается и продолжает трассировку, а с вероятностью 0.75 - погибает (рис. \ref{fig:rus_roulette}). Таким образом, если на поверхность падает 1000 фотонов, 250 отразится, сохранив энергию, а остальные 750 - погибнут. За счет такого механизма повышается точность решения, покольку в местах с высокой освещенностью будет сохраняться больше фотонов, чем в местах с низкой освещенностью; в неосвещенных областях фотоны будут отсутствовать. 
\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{rus_roulette}           
  \caption{\label{fig:rus_roulette}
           Русская рулетка.}   
\end{figure}
\noindent Нетрудно заметить, что при использовании русской рулетки в чистом виде все фотоны будут иметь единичную энергию. То есть, энергию фотона хранить не обязательно, а для того чтобы  оценить значение освещенности, достаточно либо посчитать число фотонов в заданном радиусе сбора, либо найти радиус, в котором лежит K ближайших фотонов. Мы вернемся к этому вопросу позже, когда будем рассматривать процесс сбора освещености. При этом для того чтобы выполнять правильно сбор освещенности, нам потребуется хранить позицию фотона и нормаль к поверхности в точке удара о поверхность. 

\paragraph{Сохранение фотонов в фотонной карте.} Вопрос сохранения фотонов на поверхностях (в фотонной карте) заслуживает отдельного рассмотрения. 
Фотоны сохраняются на всех диффузных (ламбертовых) поверхностях, которые они ударяют. 
Как правило фотоны используются только при вычислении части интеграла, обусловленной диффузной компонентой. 
Зеркальная компонента не может быть вычислена при помощи фотонов по той же причине, по которой прямая трассировка путей (Ligh Tracing) 
отображает зеркала и стекла черными: вероятность встретить фотон (или путь от источника), отразившийся в строго заданном направлении стремится к нулю. 
Из этого следует, что компонента, обусловленная матовым отражением (glossy reflection) все же может быть вычислена при помощи фотонных карт, однако эффективность такого расчета будет тем ниже, чем более зеркальный характер приобретает ДФО. Следует подчеркнуть, что фотон в течении трассировки может быть сохранен несколько раз - на всех диффузных поверхностях, которые он ударяет.   

\noindent Итак, Фотоны сохраняются на всех диффузных поверхностях, которые они ударяют. 
Однако, одним из часто-применяемых трюков являетcя в буквальном смысле 'пропуск' первого отскока для процедуры сохранения фотона. 
Фотоны в этом случае сохраняются только после первого отскока а освещение, получаемое в результате интегрирования фотонной карты - целиком вторичное. 
Этот метод значительно снижает число фотонов, необходимое для расчета вторичного освещения. 
Первичное освещение гораздо быстрее вычисляется при помощи трассировки лучей/путей, поэтому данный прием имеет важное значение для ускорения процесса рендеринга.




\paragraph{Построение фотонной карты.}

\noindent В самом простом случае можно хранить фотонну карту в виде обыкновенного массива. Однако это не позволит в последствие выполнять эффективный поиск ближайших фотонов. Поэтому под построением фотонной карты мы будем понимать построение ускоряющей структуры, позволяющей выполнять эффективный поиск ближайших фотонов. Существует довольно большое число способов, как делать это эффективно. Способы могут различаться для поиска k-ближайших или всех фотонов в заданном радиусе. Среди наиболее эффективны структур следует отметить пространственные хэш-таблицы и деревья, построенные при помощи эвристики VVH \cite{wald:04:VVH} аналогично эвристике SAH, рассмотренной нами ранее.

\paragraph{Сбор освещенности.} После того как фотонная карта построена, наступает стадия сбора освещенности. Из виртуальной камеры испускаются лучи и выполняется обратная или стохастическая трассировка лучей (или трассировка путей). В местах соударений луча и поверхности производится сбор освещенности. Причем, для стохастической трассировки лучей (или путей), диффузные переотражения не учитываются. Компонента, обусловленная диффузными переоражениями получается при интегрироании фотонной карты (фомула \ref{eq:phintegral4}). Диффузные переотражения в процессе обратной трассировки лучей/путeй вычисляются при использовании метода Финального Сбора, который мы рассмотрим позже. 

\begin{eqnarray}\label{eq:phintegral4}
I(\phi_{r}, \theta_{r}) \approx \dfrac{1}{\pi r^{2}} \sum_{i=1}^K R(\phi_{i},\theta_{i},\phi_{r},\theta_{r}) \Delta \Phi(\phi_{i},\theta_{i})
\end{eqnarray}

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{phradianceestimate}           
  \caption{\label{fig:phradianceestimate}
           Иллюстрация к формуле \ref{eq:phintegral4}. Cбор освещенности или интегрирование фотонной карты. В данной формуле $r$ - радиус диска, полученного при проекции сферы на поверхность. $\pi r^{2}$ - площадь диска. }   
\end{figure}

\noindent Cуществует 2 основные стратегии сбора освещенности:
\begin{enumerate}
\item Сбор k-ближаших фотонов (динамический радиус сбора), в соостветствии с формулой \ref{eq:phintegral4}. Рисунок \ref{fig:gatherphotons}.
\item Сбор всех фотонов в заданном радиусе (фиксированный радиус сбора). Данная стратегия упрощает процесс сбора и дает корректный результат в пределе (при стремлении общего числа фотонов к бесконечности), однако не обеспечивает постоянное значение $K$ в формуле \ref{eq:phintegral4}. 
\end{enumerate}

\noindent Поиск ближайших фотонов, как правило, использует специальную структуру 'max-heap' для упорядочивания списка фотонов \cite{Jensen:2007:HQR:1281500.1281593}. Однако, поиск k-ближайших заключает в своей основе более простой поиск всех фотонов в заданном радиусе, поскольку даже при поиске k-ближайших необходимо начинать с некоторого радиуса и уметь находить все фотоны в нем.
   
\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{gatherphotons}           
  \caption{\label{fig:gatherphotons}
           Cбор k-ближайших фотонов. Радиус сбора выбирается динамически (чтобы ограничить ровно k фотонов) в зависимости от плотности фотонной карты. }   
\end{figure}

\subsection{Финальный сбор}

\noindent Один из артефактов, возникающих при реализации сбора освещенности - темные края поверхностей(рис \ref{fig:fg1}, \ref{fig:fg2}). Данный артефакт возникает по причине того, что на границах поверхностей освещение собирается только с половины диска. При этом площадь, вычисляемая как $\pi r^{2}$ для целого диска остается прежней. Один из способов, позволяющий убрать темные края заключается в том, чтобы найти реальную площать элемента поверхности на котором происодил сбор при помощи вычисления площади пересечение сферы сбора и всех треугольников с нормалью, близкой к нормали в точке сбора (то есть вычислить честную проекцию сферы на поверхность). Второй метод называется Финальный Cбор (Final Gathering, FG). Вместо непосредственного сбора освещенности с фотонов, в методе финального сбора из заданной точки испускается некоторое число лучей по полусфере и освещенность собирается уже в тех местах, куда попали лучи (рисунок \ref{fig:fg5}).
 
\begin{figure}[htb]
  \centering
  \includegraphics[width=0.5\linewidth]{fg5}           
  \caption{\label{fig:fg5}
           Иллюстрация метода финального сбора. }   
\end{figure}

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

\paragraph{Карты светимости (radiosity maps).}  Процесс поиска ближайших фотонов является, как правило, узким местом финального сбора. Поскольку в данном случае при помощи самих фотонных карт достаточно получить лишь грубую оценку освещенности, существуют различные методы ускорения финального сбора, использующие этот факт. К таким методам относятся карты светимости \cite{Tabellion:2004:AGI:1186562.1015748} или октанные текстуры \cite{Vostryakov:2006:Graphicon}. Их смысл заключается в том, чтобы предрасчитать светимость в точках поверхности \cite{Tabellion:2004:AGI:1186562.1015748} (или в объеме \cite{Vostryakov:2006:Graphicon}) и запомнить ее в виде карты светимости - двумерной или трехмерной (возможно разряженной) текстуры. При попадании луча финального сбора в определенную точку поверхности поиск ближайших фотонов не выполняется. Вместо этого производится быстрая выборка значения из карты светимости.


\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{fg1}           
  \caption{\label{fig:fg1}
           Сбор освещенности с 10M фотонов (слева). Финальный сбор с 1M фотонов (справа). }   
\end{figure}

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{fg2}           
  \caption{\label{fig:fg2}
           Сбор освещенности с 10M фотонов (слева). Финальный сбор с 1M фотонов (справа). }   
\end{figure}

\subsection{Прогрессивные фотонные карты}

\noindent Для фотонных карт в классической реализации серьезной проблемой является объем потребляемой памяти. Для того, чтобы получить изображение высокого качества может потребоваться до нескольких миллиардов фотонов. Хранить все фотоны в памяти в этом случае не представляется разумным, поскольку даже если хранить фотоны в соответствии с компактным представлением в 16 байт, один миллиард фотонов займет $ 16*10^{9} \approx 15Gb $ .

\noindent Для того чтобы иметь возможноть обрабатывать произвольное число фотонов имея фиксированный объем памяти был разработан метод прогрессивных фотонных карт \cite{Hachisuka:2008:PPM:1409060.1409083}. Основная идея метода прогрессивных фотонных карт заключается в том, чтобы повторять все 4 шага (испускание фотонов из источника света, трассировка, построение фотонной карты и сбор освещенности) несколько раз для порций фотонов фиксированного размера (например,по 1 миллиону фотонв). При этом, при обработке каждой следующей порции фотонов радиус сбора и аккумулируемый поток уменьшаются в соответствии с формулой \ref{eq:ppm}. Такой метод называется Прогрессивными фотонными картами (Progressive Photon Mapping, PPM \cite{Hachisuka:2008:PPM:1409060.1409083}). Если выбрана стратегия с фиксированым радиусом сбора, то для каждой следующей порции фотнов радиус сбора уменьшается в соответствие с соотношением \ref{eq:sppm}. Такой алгоритм называется Стохастическими Прогрессивными Фотонными Картами (Stochastic Progressive Photon Mapping, SPPM \cite{Hachisuka:2009:SPP:1618452.1618487}, \cite{Hachisuka:2012:SAP:2343483.2343489}). Параметр $ \alpha $ варьируется в пределах от 0 до 1. Рекомендуемое значение - больше или равно 2/3.


\begin{eqnarray}\label{eq:ppm}
\dfrac{r_{i+1}^{2}}{r_{i}^{2}} = \dfrac{N_{i} + \alpha M_{i}}{ N_{i} + M_{i} } 
\end{eqnarray}

\begin{eqnarray}\label{eq:sppm}
\dfrac{r_{i+1}^{2}}{r_{i}^{2}} = \dfrac{i + \alpha}{i+1}
\end{eqnarray}

\noindent В соотношении \ref{eq:ppm} $N_i$ - число уже собранных фотонов. $M_i$ - число вновь добавленных фотонов (на очередном проходе). 

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{sppm_comparison1}           
  \caption{\label{fig:sppm_comparison1}
           Стохастические прогрессивные фотонные карты. 1M фотонов слева. 10M (10 проходов) - справа. }   
\end{figure}

\noindent Метод стохастических прогрессивных карт значительно проще в реализации, чем обыкновенные прогрессивные фотонные карты, поскольку он не фиксирует точки сбора (что необходимо делать в методе прогрессивных фотонных карт)\cite{Hachisuka:2009:SPP:1618452.1618487}, \cite{Hachisuka:2012:SAP:2343483.2343489}. Причем, в \cite{Hachisuka:2009:SPP:1618452.1618487} было показано, что стохастические прогрессивные фотонные карты обладают лучшей сходимостью. Результаты для различных порций фотонов могут быть скомбинированы простым усреднением, 
как и в обыкновенной трассировке путей (с.м. 'прогрессивное вычсление интеграла'). 
Алгоритм стохастических прогрессивных фотонных карт выглядит следующим образом:

\begin{enumerate}
\item Заранее выделяем некоторый объем памяти для хранения достаточно большой порции фотонов (например 1 миллион фотонов).

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

\item Строим ускоряющие структуры для быстрого поиска ближайших фотонов.

\item Трассируем пути из виртуальной камеры. В точках пересечения с диффузными поверхностями собираем освещенность по формуле \ref{eq:phintegral4}. При этом используем фиксированный радиус сбора $r_i$. Цвет, как и при обыкновенной трассировке путей накапливается в пикселах.

\item Уменьшаем радиус в соответствии с формулой \ref{eq:sppm}. Переходим к шагу 2.

\end{enumerate}

\subsection{Некоторые важные оптимизации}

\paragraph{Три фотонные карты} В целях эффективности (повышения точности) в \cite{Jensen:2007:HQR:1281500.1281593} предлагается хранить 3 раздельные фотонные карты:

\begin{enumerate}
\item Глобальная фотонная карта. Данная карта содержит фотоны, прошедшие путями вида $ L(S|D|G|V)^{+}D $ и используется для вычисления диффузного освещения.
\item Каустическая фотонная карта. Данная карта содержит фотоны, прошедшие путями вида $ LS^{+}D $ и используется для вычисления каустиков.
\item Объемная фотонная карта. Данная карта содержит фотоны, прошедшие путями вида $ L(S|D|G|V)^{+}V $ и используется для вычисления объемных эффектов.
\end{enumerate}

\noindent Разделение фотонных карт позволяет обрабатывать каустики и диффузное освещение (а также объемное рассеивание) при помощи различных порций фотонов. Это особенно важно при использовании финального сбора, поскольку в этом случае для диффузного освещения с финальным сбором нужно гораздо меньше фотонов чем для каустиков.  

\paragraph{Определение видимости.} Рассмотрим изображение сцены классной комнаты на рис. \ref{fig:class_fg}. Для рассчета освещения на этой сцене фотоны испускались снаружи. При этом, лишь небольшая часть фотонов проникла внутрь комнаты через окно. Остальные фотоны сохраняются с обратной стороны стен и потолка, приводя к нэффективному расходованию вычислитльных ресурсов. Для того чтобы не сохранть фотоны на снаружи комнаты, необходимо перед тем, как производить трасровку фотонов, пометить некоторым образом поверхности, для которых фотоны сохраняются в фотонной карте. 

\noindent Простейший способ заключается в трассирове некоторого небольшого числа путей для каждого пиксела из виртуальной камеры и прямой пометки пересекаемых треугольников. При этом, каждую сторону треугольника необходимо помечать отдельно. Чтобы идентифицировать сторону треугольнка, можно использовать векторное произведение для вычисление нормали. Поскольку вершины треугольника хранятся в памяти в заданном порядке, такой способ однозначно позволяет идентифицироать сторону по нормали. Недостаок данного метода заключаетс в том, что некоторые тонкие примитивы (треугольники) могут быть пропущены алгоритмом. 
Хотя, если используется финальный сбор, это не имеет значения посколько такие примитивы вносят почти нулевой вклад в результирующее значение интеграла. Т.е. чем тоньше треугольник, тем больше вероятность не попасть в него в процессе трассировки путей из виртуальной камеры. Но также уменьшается вероятность попасть в данный треугольник в процессе сбора освещенности.

\noindent Еще один способ заключается в трассировке так называемых частиц видимости \cite{Suykens:2000:DCP:647652.732120}. 
          Частицы видимости представлят из себя фотоны, испускаемые из камеры. 
		  Таким образом, сначала строится фотонная карта видимости, в которой все фотоны трассиртся из виртуальной камеры, а затем, используя карту видимости, строится фотонная карта освещенности.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{class_fg}           
  \caption{\label{fig:class_fg}
           Комната, освещенная солнцем через окно. Лишь $\frac{1}{7}$ часть выпущенных из источника света фотонов попадает внутрь комнаты.  }   
\end{figure}

\paragraph{Карты проекций.} При оптимизированной реализации построения ускоряющей структуры и сбора освещенности, одним из наиболее ресурсоемких этапов в фотонных картах становится трассировка фотонов. Причина низкой скорости трассировки заключается не том, что медленно происходит сам процесс трассировки одиночного фотона, а в том, на сложных сценах значительная часть испускаемых из источника света фотонов погибает в процесе трассировки поскольку фотоны не достигают видимых областей. Это особенно верно при реализации оределения видимости, рассмотренном выше и при трассировки каустических фотонов. 
На сцене из рисунка \ref{fig:class_fg} только $\frac{1}{7} $ часть фотонов выпущенных из источника света попала внутрь комнаты. На классической сцене 'cornell box' (рис. \ref{fig:montecarlo1}) только $ \frac{1}{10} $ часть каустических фотонов, выпущенных из источника света, попала в зеркальный или стеклянный шары, образуя каустики. Русская рулетка дополнительно усугубляет ситуацию, стохастически убивая фотоны.
 
\noindent Карты проекций \cite{Jensen:2007:HQR:1281500.1281593} - это метод, позволяющий не выпускать фотоны в те области, где они гарантированно погибают, не сохраняясь ни на одной поверхности. Идея этого метода заключается в том, чтобы спроецировать сцену на источник света (например при помощи растеризации в кубическую текстурную карту) и пометить на этой проекции 'мертвые области'. После чего в эти мертвые области фотоны не трассируются совсем, погибая еще при рождении на источнике света.

\paragraph{Приминение фильтрации.} Довольно простым и эффективным способом, позволяющим улучшить четкость каустиков при помощи метода фотонных карт является применение фильтра на основе ядра Епанечникова \cite{Schjoth:2007:PD:1321261.1321293}. При таком подходе освещенность при сборе домножается на определенный вес в зависимости от расстояния до фотона (формула \ref{eq:ph_filter}, нормирующий множитель '2' вынесен за сумму).
 
\begin{eqnarray}\label{eq:ph_filter}
w_{i} = 1 - d^2/r^2 \\
I(\phi_{r}, \theta_{r}) \approx \dfrac{2}{\pi r^{2}} \sum_{i=1}^K R(\phi_{i},\theta_{i},\phi_{r},\theta_{r}) \Delta \Phi(\phi_{i},\theta_{i})w_{i}
\end{eqnarray}

\noindent Возможно использование и других типов фильтров. Например в \cite{Jensen:2007:HQR:1281500.1281593} предлагается использовать 'cone filter' и фильтра Гаусса. Однако, фильтрация на основе ядра Епанечникова эффективнее в вычислительно плане, поскольку не требует вычисления квадратного корня и оперирует отношениями квадратов расстояний. 

\subsection{Объемные фотонные карты}

\noindent При помощи алгоритма фотонных карт можно достаточно эффективно вычислять эффекты основанные на объемном рассеянии света. В этом случае фотоны сохраняются прямо в объеме, а сбор освещенности происходит при помощи 'марширования по лучу' (ray marching).

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{volume_phmap1}           
  \caption{\label{fig:volume_phmap1}
           Трассировка фотонов внутри объема \cite{Hachisuka:2012:SAP:2343483.2343489}. }   
\end{figure}


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{volume_phmap2}           
  \caption{\label{fig:volume_phmap2}
           Сбор освещенности при помощи марширования вдоль луча \cite{Hachisuka:2012:SAP:2343483.2343489}. }   
\end{figure}


\begin{eqnarray}\label{eq:sppm2}
\dfrac{r_{i+1}^{3}}{r_{i}^{3}} = \dfrac{i + \alpha}{i+1}
\end{eqnarray}

\noindent Для стохастических прогрессивных объемных фотонных карт соотношение \ref{eq:sppm} модифицируется в \ref{eq:sppm2} с учетом того что процесс происходит в объеме.

\subsection{Резюме по фотонным картам}

\noindent Фотонные карты - метод дающий смещенную оценку (смещенный метод). 
          Прогрессивные фотонные карты и стохастические прогрессивные фотонные карты -  также методы смещенные, но при этом состоятельные (сходимость по вероятности) \cite{Hachisuka:2008:PPM:1409060.1409083}, \cite{Hachisuka:2008:PPM:1409060.1409083}). 
		  На практике смещение проявляется на изображении в виде цветных пятен (как и в кэше освещенности). 
		  В фотонных картах цветные пятна достаточно долго не удаляются с изображения в процессе рассчета и изображение лишенное пятен может потребовать 
		  обработки нескольких миллиардов фотонов. При таком большом числе фотонов гораздо более простая в реализации прямая и/или обратная трассировка 
		  путей даст более точное и приемлимое для человеческого глаза решение. 

\noindent Финальный сбор позволяет снизить число необходимых фотонов примерно на порядок. 
          Однако, он может быть использован только для вычисления диффузного вторичного освещения и его корость сравнима со скоростью обыкновеной трассировки путей. Для эффективной реализации финального сбора рекоммендуется использовать карты светимости \cite{Tabellion:2004:AGI:1186562.1015748} или октанные текстуры \cite{Vostryakov:2006:Graphicon}. Последний метод предпочтительнее, поскольку для него не нужно реализовывать отображение поверхностей трехмерных объектов на двумерную текстуру (что в общем случае является нетривиальной задачей).

\noindent Тем не менее, фотонные карты наряду с методом соединения вершин \cite{Georgiev:2011:BLT:2077378.2077412} являются одним из самых эффективных алгоритмов для вычисления путей вида $S^{+}DS^{+}$ (каустики, видимые через стекло или зеркало). Такие пути, как правило, значительно хуже вычисляются трассировкой путей.



\section{Кэш Освещенности (Irradiace Cache)}

\noindent Кэш Освещенности (Iradiance Cache, IC) – это не способ вычисления интеграла освещенности. Это лишь способ ускорения вычисления этого интеграла на множестве точек. Алгоритм был впервые представлен в работе \cite{Ward:1988:RTS:378456.378490}. 
Основная идея заключается в том, что вторичное освещение разделяется на две компоненты – низкочастотную (диффузную) и высокочастотную (отражающую). 
Низкочастотная компонента на изображении и в трехмерном пространстве меняется плавно, поэтому ее можно вычислить каким-либо из методов лишь в очень небольшом числе точек, а в остальных - интерполировать \cite{Gautron:2009:PGI:1855086}. Рисунки 1 и 2 демонстрируют пример использования кэша освещенности. 

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

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{icrecords}           
  \caption{\label{fig:icrecords}
           Освещение вычисленное при помощи кэша освещенности (слева). Визуализированные точки кэша (справа). }   
\end{figure}

\noindent Различают кэш освещенности в трехмерном пространстве сцены (world space) и в экранной плоскости (screen space).
          Реализация кэша в экранной плоскости проще и эффективнее. Интерполяция также производится в пространстве экрана. 
		  Однако, кэш освещенности в пространстве экрана может быть использован по очевидным причинам только для первично видимых поверхностей.

\noindent Как правило, кэш освещености вычисляется на лету. Каждый раз при вычислении вторичного диффузного освещения делается попытка выборки из кэша. 
Если значение может быть выбрано так, что ошибка меньше допустимой величины, оно получается как результат интерполяции ближайших точек кэша. 
Если ошибка больше допустимой величины, значение честно вычисляется с помощью трассировки лучей по всем направлениям (или финального сбора из фотонной карты) 
и полученная точка сохраняется в КЭШе \cite{Gautron:2009:PGI:1855086}.

\begin{algorithm}[H]
\SetAlgoLined
\KwData{p - точка на поверхности, n - нормаль к поверхности в точке p}
\KwResult{Значение освещенности}

\eIf{InterpolationIsPossibe(p,n)}{
   return IrradiacneCacheLookUp(p,n) \;
}{
  $I \leftarrow EvalIntegral(p,n) $ \;
  
  $R \leftarrow CalcValidityRadius(p,n) $\;
  
  $ InsertRecordInCache(I,R,p,n) $ \;
  
  return I\;
 }


\caption{Алгоритм кэширования освещенности}
\label{ic_alg}
\end{algorithm}


\noindent Алгоритм кэширования вторичного освещения содержит 2 узких места - расчет собственно вторичного освещения в точках кэша и финальный алгоритм интерполяции освещения во всех остальных точках.

\subsection{Вычисление вторичного освещения}

\noindent Стандартным способом вычисления вторичного освещения в точках кэша является монте-карло трассировка лучей или финальный сбор в сочетании с фотонными картами. Однако в \cite{Gautron:2009:PGI:1855086} предлагается вычислять вторичную освещенность в точках с помощью растеризации в кубическую текстурную карту. Данная операция может производиться очень быстро т.к. имеет аппаратную поддержку даже в довольно старых видеокартах (обычно используется для симуляции отражений на объектах).

\noindent Несмотря на преимущество по скорости, недостаток этого алгоритма в том, что он не может учитывать зеркальных отскоков света и может быть в ряде случаев медленнее трассировки лучей на больших сценах. Причина этого заключается в том, что при растеризации вся геометрия обрабатывается графическим процессором. Она проходить весь графический конвейер: преобразование вершин, возможно тесселяцию и геометрические преобразования над примитивами, интерполяция атрибутов вершин и вычисление значений освещенности в пикселах. Допустим, геометрия – это некоторое здание, разделенное на комнаты. Всё здание содержит 2 миллиона вершин. Пусть виртуальный наблюдатель находится в одной из комнат. Несмотря на то, что комната может составлять 10-20\% всей геометрии, при вычислении освещения в точке каждый раз потребуется пропускать через графический конвейер все 2 миллиона вершин. Тот факт, что 90\% сцены может быть закрыто стеной или каким-то близко-расположенным объектом никак не используется. В то же время, скорость трассировки лучей при правильной реализации зависит от размера входных данных сублинейно (обычно логарифмически), а не линейно, как в случае растеризации. 

\noindent Отдельного внимания заслуживает вопрос точности вычисления освещенности в точках кэша. Допустим для вычисления освещенности при помощи трассировки путей в пикселах на некоторой сцене будет достаточно около 1000 путей на пиксел. Тогда для вычисления освещенности в точках кэша потребуется может потребоваться большая точность (порядка 4000 путей). Недостаточная точность вычисления для попиксельной трассировски путей приводит к шуму на изображении, который в некоторой степени фильтруется глазом и мало заметен. В случае кэша освещенности, недостаточная точность приводит к появлени цветных пятен, которые гораздо более заметны для человеческого глаза, чем шум.

\subsection{Вычисление радиуса валидности}

\noindent Одним весьма неочевидным моментом при реализации кэша освещенности является вычисление радиуса валидности точки кэша 
(переменная R, алгоритм \ref{ic_alg}). Подразумевается, что каждая точка иррадианс кэша может быть использована для интерполяции только в 
некоторой определенной области. Есть по крайней мере 2 причины, по которым радиус валидности точки нужно ограничивать. 
Первая причина - стремление избежать артефактов. Вторая - скорость выборки из кэша. 
Если слишком много точек будут затрагиваться при каждой выборке, интерполяция станет медленной.

\noindent Обычно область валидности ограничивают сферой с центром в точке иррадианс кэша и некоторым радиусом, называемым радиусом валидности (хотя есть работы где для этой цели используются эллипсоиды). В [3] для вычислениия радиуса валидности было использовано 5 различных эвристик. Каждая из них имеет недостатки, но все вместе они дают удовлетворительный результат.

\noindent Из самых простых - используется эвристика расстояния до ближайших поверхностей. Эта эвристика вычисляется исключительно из геометрических соображений. При оценке освещенности, для всех трассируемых лучей запоминается расстояние, на котором они ударились о поверхность и из этих расстояний берется наименьшее. Однако при использовании этого подхода, в процессе нахождения минимума важно не принимать в расчет лучи, имеющие маленький угол с тангент-плоскостью (рис. \ref{fig:icradius1}).

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{icradius1}           
  \caption{\label{fig:icradius1}
           Точки расположенные на вогнутой поверхности. Если принимать в расчет лучи, идущие перпендикулярно нормали, радиус валидности получится слишком маленьким. }   
\end{figure}

\noindent В противном случае на вогнутых криволинейных поверхностях радиус валидности точки всегда будет равным нулю. Альтернативный метод реализации эвристики расстояния до поверхностей заключается в том, чтобы не искать минимальное расстояние а вычислить так называемое "среднее грамоническое расстояние" (формула \ref{eq:rhmd}).

\begin{equation} \label{eq:rhmd}
R^{HMD}_{i} = \dfrac{N}{ \sum_{i=1}^N \dfrac{1}{r_{i}} }
\end{equation}

\noindent Число N в формуле \ref{eq:rhmd} - это количество лучей используемых при сэмплировалнии полусферы. $ r_{i} $ - расстояние от точки кэша до поверхности в которую ударился соответствующий луч. Формула \ref{eq:rhmd} для среднего гармонического расстояния характерна тем, что чем меньше расстояние $ r_{i} $ тем больше оно влияет на финальный результат. Дополнительно размер радиуса валидности рекомендуется ограничивать снизу размером 2-4 пикселей (спроецированных на сцену в world space) и сверху размером 20-64 пикселей. Конкретные числа могут варьироваться.

\subsection{Структуры данных и интерполяция}

\noindent Для хранения точек кэша освещенности обычно используется специальное октодерево со множественными ссылками (multiple reference octree) \cite{Gautron:2009:PGI:1855086}. Точки кэша излучения в таком октодереве хранятся только в листьях, но каждый лист хранит список всех сфер, которые он пересекает. Основной плюс такого октодерева в том, что возможно осуществить простой нерекурсивный поиск от корня к листу с последующим перебором всего лишь нескольких точек кэша освещенности.

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{mroctree2}           
  \caption{\label{fig:mroctree2}
           Октодерево с одиночными ссылками (слева) и окто-дерево с множественными ссылками (справа). Cтрелки обозначают ссылки. }   
\end{figure}

\paragraph{Интерполяция.} В заданной точке трехмерного пространства, при помощи поиска в октодереве необходимо найти все сферы 
(с центрами в точках иррадианс кэша и радиусами равными радиусам валидности записей кэша), которые пересекают данную точку. 
Затем, пройдя в цикле по всем точкам, нужно собрать с них освещение. 
Одна из практичных формул была предложена Tabellion-ом и Lamorlette в \cite{Tabellion:2004:AGI:1186562.1015748}. 

\subsection{Резюме по кэшу освещенности}

\noindent Кэширование вторичного освещения - весьма эффективная техника ускорения вычисления глобально освещения, позволяющая снизитиь количество необходимых вычислений практически на 3 порядка (для больших разрешений). Допустим имеется изображение размером 1024x1024 пиксела. Если вычислять вторичное освещение в каждой точке, то потребуется суммарно протрассировать порядка $10^9$ лучей только для вычисления вторичной освещенности. При использовании кэша освещенности количество точек, в которых необходимо вычислить вторичную освещенность сокращается с миллиона до нескольких тысяч. Соответственно, число лучей падает с $10^9$ до $10^6$ - $10^7$.

\noindent Однако, кэширование освещенности имеет определенные недостатки:

\begin{enumerate}
\item Кэш освещенности хорошо работает только на сценах с преимущественно гладкими поверхносями. На сценах с большим количеством геометрических деталей (например карты нормалей для имитации микрорельефа) кэш освещености не даст преимущества и может даже замедлить процесс редеринга.
\item Кэш освещенности достаточно сложно сделать параллельным. 
\item Кэш освещенности дает мерцание при анимации. Если стохастическая трассировка путей дает шум, который эффективно удаляется при помощи размытия в движении, то пятна вызванные кэшированием освещенности удалить гораздо труднее. 
\item Кэш освещенности все же снижает точность решения и является смещенным методом. Чем ниже точность, тем больше выигрыш в скорости, но тем более заметны артефакты.
\end{enumerate} 

















