\chapter{Анализ литературных источников}

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

\section{Алгоритмы на центральном процессоре}


\subsection{Излучательность}

\noindent Метод излучательности (Radiosity) в компьютерной графике был впервые представлен в работе \cite{Goral:1984:MIL:800031.808601}. Идея алгоритма Излучательности заключается в том, чтобы разбить поверхность сцены на множество маленьких кусочков - патчей; после чего производится расчет переноса энергии между отдельными патчами.

\noindent Алгоритм вводит следующие ограничения на входную сцену:

\begin{enumerate}
\item Все поверхности имеют ДФО Ламберта (отражают свет равномерно во все стороны).
\item Сцена замкнута, т.е. суммарная энергия в системе должна сохраняться.
\end{enumerate} 

\noindent На первом этапе алгоритма все поверхности сцены делятся на патчи. Патч – это элементарная единица поверхности. Дискретизация поверхности на патчи позволяет заменить интеграл в уравнении освещенности на конечную сумму интегралов специального вида. Каждый такой интеграл, называемый форм-фактором, задает взаимное влияние двух отдельных патчей (т.е. сколько энергии переходит от одного патча к другому). Если мы знаем значение каждого форм-фактора (для каждой пары патчей), то процесс синтеза изображений сводится к решению системы линейных алгебраических уравнений. Основная трудность в алгоритме – расчет форм-факторов \cite{Lebedev:Radiosity_Hist:2011}.

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{lebedevrad}           
  \caption{\label{fig:lebedevrad}
           Иллюстрация алгоритма излучательности и получаемой системы уравнений \cite{Lebedev:Radiosity_Hist:2011}. }   
\end{figure}

\noindent Форм фактор, в свою очередь, можно расчитать для каждой пары патчей при помощи Аналогии Нуссельта \cite{Cohen:1993:RRI:154731} зная только геометрию сцены. 

\noindent Основными недостатками метода излучательности являются:

\begin{enumerate}
\item Поддержка только Ламбретвовых материалов.
\item Артефакты на изображении (т.е. низкая точность) в местах перехода между патчами. 
\item Большой объем потребляемой памяти и вычислений при увеличении точности и необходимости расчета резких переходов свет-тень даже при адаптивном разбиении (разный размер патчей в зависимости от того, насколько быстро на поверхности происходит смена света и тени).
\end{enumerate}

\noindent Тем не менее излучательность достаточно часто используется для вычисления вторичного диффузного освещения, когда высокая точность не требуется. 
Например, в известной компьютерной игре Battlefield 3 для предрассчета вторичного освещения \cite{Battlefield3}.

\subsubsection{Мгновенная Излучательность и VPL}

\noindent Основная идея метода мгновенной излучительности (Instant Radiosity; рис. \ref{fig:instrad}) заключается в том, чтобы апроксимировать вторичное освещение в сцене при помощи некоторого числа виртуальных точечных источников света (VPL, Virtual Point Light) \cite{Keller:1997:IR:258734.258769}. 
При этом сами виртуальные источники могут расставляться по сцене различными способами. Освещение от точечных источников света может быть вычислено при помощи трассировки лучей или метода карт теней. Однако, заслуживает внимания вопрос эффективной реализации при большом числе источников. Одно из возможных решений основывается на построении специального дерева для источников света \cite{Walter:2005:LSA:1186822.1073318}. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{instrad}           
  \caption{\label{fig:instrad}
           Иллюстрация метода мгновенной излучательности. }   
\end{figure}

\noindent Мгновенная излучательность породила целое направление исследования решения проблемы глобальной освещенности, которое называется 'Методом Многих Источников' (Many Light Method) \cite{Krivanek:2012:ORR}. Идея вычислять освещение при помощи множества источников вторичного освещения применяется в основном если предпочтение отдается скорости, но не точности решения \cite{Krivanek:2012:ORR}. Здесь следует обратить внимание на работы \cite{Walter:2005:LSA:1186822.1073318}, \cite{Walter:2012:BL:2185520.2185555}, \cite{Ou:2011:LMS:2024156.2024213}, которые значительно улучшают точность решения проблемы глобальной освещенности при помощи виртуальных источников.

\subsection{Подход на основе облаков точек}

\noindent Данный подход (называется Point Based Global Illumination или PBGI) позволяет используется в программном решении Render Man компании Pixar \cite{Christensen:2010:PBGI}. Он хорошо зарекомендовал себя в производстве анимационным компьютерных фильмов (т.е. мультфильмов), поскольку позвозволяет достаточно быстро получить свободное от шума (что особенно важно для анимации) приближение диффузного глобального освещения.

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

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{surfels}           
  \caption{\label{fig:surfels}
           Представление освещенной поверхности в виде множества сюрфелов \cite{Christensen:2010:PBGI}. }   
\end{figure}

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

\begin{enumerate}
\item Если сюрфелы в узле мало отличаются друг от друга по нормали, вычисляется один сюрфел с усредненной позицией, нормалью, площадью и цветом.
\item Если сюрфелы сильно различаются по нормали, значение излучаемой освещенности сохраняется в представлении, основанном на сферических гармониках \cite{SHLighting}. Причем оказывается, что коэффициенты гармоник для узла могут быть вычеслены как сумма коэффициентов отдельных сюрфелов.
\end{enumerate}

\noindent Описанный выше метод применяется иерархически, строя более высокие уровни окто-дерева. При этом, инетерсным фактом является то что в родительских узлах коэффициенты гармоник могут быть вычислены точно так же простым суммирование коэффициентов дочерних узлов \cite{Christensen:2010:PBGI}.

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

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

\noindent Далее, в зависимости от расстояния до сюрфела, их растериция производилась с использованием 3 стратегий: 

\begin{enumerate}
\item Если от целевой точки расстояние до сюрфела большое, используется аппроксимация в виде узла окто-дерева, сохраняющая единственный сюрфел или гармоники.

\item Для сюрфелов, для которых аппроксимация узлом недостаточна использовалась растеризация отдельно каждого сюрфела.

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

\end{enumerate}

\noindent Для поддержки множественных диффузных отскоков света в \cite{Christensen:2010:PBGI} предлагается повторять алгоритм в течении нескольких проходов. Учитывая этот, метод PBGI можно считать в некотором роде аналогом иерархической излучательности \cite{Hanrahan:1991:RHR:127719.122740}. В сочетании с собственно алгоритмом PBGI в \cite{Christensen:2010:PBGI} использовался модифицированный кэш освещенности для ускорения вычисления глобального освещения и отмечается его высокая эффективность.

\subsection{Методы на основе Монте-Карло интегрирования}

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

\subsubsection{Методы на основе MLT}

\paragraph{Kelemen style MLT} Оригинальный алгоритм MLT, разработанный Э. Вичем и Л. Гибасом \cite{Veach:1997:MLT:258734.258775} и обсуждаемый 
в первой главе использует мутации в мировом пространстве или, иначе говоря, мутации в пространстве путей. 
Причем, одной общей стратегии мутации для всех видов путей не существует. 
Такая стратегия строится путем комбинирования большого числа различных стратегий для отдельных случаев (мутации в пространстве линз, мутации каустик, мутации на источнике света и.т.д.). Это значительно осложняет реализацию, поскольку:
 
\begin{enumerate}
\item Усложняется процесс вычисления вероятностей перехода $T(X|Y)$ и $T(Y|X)$. 
      Если вероятности будут вычислены неправильно, нарушится условие, гарантирущее сходимость к правильному решению (detailed balance).
\item Параметры алгоритма мутации необходимо настраивать в зависимости от входных данных.
\end{enumerate} 

\noindent В работе \cite{Kelemen02asimple} предложен простой и стабильный (независящий от входных данных) метод MLT, основанный на мутациях 
путей в 'первичном пространстве выборок' - единичном кубе, в котором генерируются случайные числа для метода Монте-Карло. 
При этом, в работе проводятся некоторые выкладки, позволяющие значительно упростить вычисление вклада текущей выборки (текущего пути) в 
изображение, которое фактически удалось свести к формулам \ref{eq:kmlt_weight1}, \ref{eq:kmlt_weight2}:

\begin{equation} \label{eq:kmlt_weight1}
w_{new} = \dfrac{a + large\_step}{(I/b+p_{large})M}
\end{equation}

\begin{equation} \label{eq:kmlt_weight2}
w_{old} = \dfrac{1 - a}{(I_{old}/b+p_{large})M}
\end{equation}
 
\noindent В формулах \ref{eq:kmlt_weight1}, \ref{eq:kmlt_weight2} приняты следующие обозначения:

\begin{itemize}
\item $ w_{new} $ - вес нового пути. 
\item $ w_{old} $ - вес старого пути. 
\item $ a = min(1, \dfrac{I_{old}}{I}) $. 
\item $ I = lum(color) $ - яркость текущего пути.
\item $ I_{old} = lum(colorPrev) $ - яркость предыдущего принятого пути.
\item $ large\_step $ - переменная, которая равна 1 если текущий путь сделал так называемый 'большой шаг' и 0 в противном случае. Большой шаг соответствует генерированию полностью случайного вектора параметров пути (то есть все числа перегенирируются без использования мутаций).

\item $ b $ - константа нормализации, такая что плотность вероятности под-интегральной функции $ p =\dfrac{I^*}{b}  $. Под $I^*$ в данном случае понимается проинтегрированная яркость пиксела (в отличие от обозначений $I$ и $I_{old}$ для яркости отдельных выборок).
\item $ p_{large} $  - вероятность так называемого 'большого' шага.
\item $ M $  - номер текущей выборки.

\end{itemize}

\noindent Далее алгоритм генерирует случайное число с равномерным распределением от 0 до 1 и:

\begin{itemize}
\item Если это число меньше $a$, мутация принимается. Новый путь (его вклад и вес) запоминается в качестве старого, а старый возвращается в качестве результата. 
\item Если это число больше или равно $a$, мутация отвергается и возвращаются вклад и вес нового пути.
\end{itemize}

\noindent Таким образом, если мутация принимается, вычисление вклада принятого пути откладывается до следующей итерации, а 
в результирующее изображение вносит вклад предыдущей принятый путь. 
Если же мутация отвергается, то в результирующее изображение вносится вклад отвергнутого пути. 
Следует отметить также, что за счет настраивания вероятности $p_{large}$ можно улучшить свойства алгоритма на 
темных участках (где оригинальный MLT работает, как правило, хуже чем обыкновенная Монте-Карло трассировка путей поскольку 
отводит на темные области меньшее число выборок) \cite{Kelemen02asimple}.


\paragraph{Energy Redistribution Path Tracing (ERPT)} Алгоритм ERPT \cite{Cline:2005:ERP:1073204.1073330} является гибридом MLT и идеи фильтрации, о которой пойдет речь далее. В отличие от MLT, ERPT использует кототкие последовательности мутаций и перераспределяет энергию путей на плоскости изображения. В отличие от фильтров и кэша освещенности, ERPT позволяет получить несмещенное решение. Идея, которая лежит в основе алгоритма ERPT называется 'поток энергии' (energy flow) и основывается на том факте, что интегралы освещенности для соседних пикселей, как правило, кореллируют. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{eflow}           
  \caption{\label{fig:eflow}
           Иллюстрация идеи 'потока энергии' в алгоритме ERPT. Поток энергии связывает колеррилующие интегралы и устанавливает перенос энергии между ними. }   
\end{figure}


\noindent Предположим, что мы имеем два кореллирующих интеграла $I_1$ и $I_2$ на соответствующих областях определения их под-интегральных функций  $\Omega_1$ и $\Omega_2$ (рис. \ref{fig:eflow}). Далее предположим, что в процессе интегрирования $I_1$ в $\Omega_1$ была найдена точка $\overline{x}$ c высокой значимостью. Поскольку $I_1$ и $I_2$ кореллируют, высока вероятность найти на области определения $\Omega_2$ недалеко от точки $\overline{x}$ точку $\overline{y}$, которая также будет иметь высокую значимость. В действительности эта же самая идея лежит в основе MLT, однако, в алгоритме ERPT ослаблаются ограничения на реализацию потока энергии (рис. \ref{fig:detailed_balance}). Вместо детального баланса потока энергий используется так называемый общий баланс. Общий баланс потока энергий требует чтобы суммарная входящая и выходящая энергия для точки были равны. Тогда как детальный баланс является более строгим ограничением и требует чтобы для каждой пары точек энергия перетекаемая между ними сохранялась \cite{Cline:2005:ERP:1073204.1073330}.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{detailed_balance}           
  \caption{\label{fig:detailed_balance}
           Общий (general balance) и детальный (detailed balance) баланс потока энергий. }   
\end{figure}


\noindent ERPT применят так называемый фильтр сбалансированного потока энергий (balanced energy flow filters) - набор правил, при помощи которых поток энергии может быть реализован таким образом, чтобы условие общего баланса сохранялось.


\noindent Сравнения, предоставленные в работе \cite{Cline:2005:ERP:1073204.1073330} позволяют сделать вывод о том, что методы MLT и ERPT - приблизительно одинаковы по эффективности.

\paragraph{Gradient Domain MLT} В работе \cite{Lehtinen:2013:GML:2461912.2461943} было предложено расширение оригинального алгоритма MLT, позволяющее улучшить сходимость метода. Алгоритм сначала вычисляет вертикальную и горизонтальную разницу (градиенты) изображения  а также грубое приближение самого изображения при помощи обыкновенного алгоритма MLT, после чего производится реконструкция финального изображения решая уравнение пуассона.

\subsubsection{Метод соединения вершин (Vertex Merging)}

\noindent В работе \cite{Georgiev:2011:BLT:2077378.2077412} с целью улучшения алгоритма двунаправленной трассировки путей на путях вида $SDS$ был предложен метод соединения вершин (Vertex Merging, VM). Идея основывается на том, чтобы соединять пути на диффузной поверхности даже если они попали не строго в 1 точку, а в находятся друг от друга на расстоянии некотороно небольшого радиуса $r$(рис \ref{fig:vertexmering}). Данный метод может рассматриваться как расширение финального сбора, при котором вместо фотонов храняться пути от источника света, а вместо сбора непосредственно освещенности применяется полностью аналогичное BDPT соединения вершин пути сбора и путей из источника, закончившихся на той же поверхности в некотром радиусе от точки сбора.

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{vertexmering}           
  \caption{\label{fig:vertexmering}
           Иллюстрация алгоритма VM. Соединение пути в классическом BDPT изображено слева. Cоединение в VM - справа. }   
\end{figure}

\subsubsection{Регуляризация в пространстве путей}

\noindent В работе \cite{Kaplanyan2013PSR} предложена похожая на метод соединения вершин идея о том, что можно соединять точки различных путей для переотражений вида $SDS$ даже если при зеркальном отражение ($S$) попадание в точку где произошло диффузное отражение ($D$) неточно (рис \ref{fig:regularisationkapl}). Однако, данный метод отличает тем, что используеют преобразование $S$ отражения в $D$ отражение и может быть применен ко всем видам Монте-Карло интеграторов, а не только к BDPT, как метод соединения вершин. 

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

\noindent В работе \cite{Kaplanyan2013PSR} отмечается что угол регуляризации (помеченный на рисунке \ref{fig:regularisationkapl} серым треугольником) может быть вычислен из параметра радиуса сбора в методе VM (то есть сам способ соединения похож на способ в методе VM). Однако, в отличие от VM, в данном методе будут соединяться не все пути, а только так называемые 'несэмплируемые' пути. То есть те пути, которые не могут быть получены без внесения смещения методами Монте-Карло (в любом виде - PT, BDPT, MLT, ERPT и.т.д. ). Это позволяет уменьшить смещение(bias) в получаемом решении. Помимо вышеперечисленных отличий, важным результатом метода \cite{Kaplanyan2013PSR} по сравнению с предыдущими работами является возможноть получения изображения точечных (или бесконечно малых) источников света при многократных зеркальных переотражениях.



\subsection{Методы на основе фотонных карт}

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


\subsubsection{Оценка освещенности}

\noindent Алгоритм фотонных карт для решения проблемы глобальной освещенности породил подзадачу 'оценки плотности фотонов' (Photon Density Estimation) \cite{Hachisuka:2012:SAP:2343483.2343489}. Для решения этой задачи применяются различные подходы. 

\noindent В работе \cite{Schjoth:2007:PD:1321261.1321293} был предложен метод, позволяющий значительно снизить смещение (bias) в получаемом решении и сделать изображение каустик более четкими. Идея метода состоит в том, чтобы хранить дополнительно на каждый фотон так называемый 'дифференциал фотона' (photon differential) по аналогии с дифференциалами лучей \cite{Igehy:1999:TRD:311535.311555} (рис. \ref{fig:photondiff}). Дифференциалы фотонов позволяют апроксимироавть одним лучом целый пучок фотонов. Именно этот апрроксимируемый пучок фотонов и называется дифференциалом фотона в работе \cite{Schjoth:2007:PD:1321261.1321293}.  

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{photondiff}           
  \caption{\label{fig:photondiff}
           Трассировка луча и его дифференциальных векторов от точки p к точке p'. }   
\end{figure}

\noindent Трассировка таких пучков технически не отличается от обыкновенной трассировки фотонов. При взаимодействии с поверхностью дифференцыалы фотонов пересчитываются аналогично диффернциалам лучей в работе \cite{Igehy:1999:TRD:311535.311555}. Дифференциалы используются далее во время сбора для повышения точности. Обыкновенные фотоны вносят вклад в интеграл освещенности на поверхности в виде круга. Дифференциалы фотонов вносят вклад в виде эллипса, делая, таким образом, фильтрацию анизотропной. Метод дифференциалов фотонов требуют хранения допольнительно 24 байт на фотон и увеличивает объем вычислений, выполняемых при сборе за счет умножения матрицы 3x3 на вектор разницы между позицией фотона и точкой сбора. Тем не менее, повышение точности решения значительное, что позволяет считать метод практичкным.
 
\noindent В работе \cite{spencer:prorelax} предолжен подход, называемый 'прогрессивной релаксацией фотонов' (progressive photon relaxation). Ключевым моментом метода является использование фиксированного объма памяти и фиксированного числа сохраняемых фотонов, называемых в работе 'устойчивыми фотонами' (persistent photons). Метод основан на построении диаграмм Вороного с последующей регуляризацией ячеек диаграммы на основе информации о вновь поступающих, временных (transient), фотонах. В среднем по эффективности метод не на много превосходит стохастические прогрессивные фотонные карты и является достаточно сложным в реализации. Преимущество метода релаксации заключается том, что он работает в мировом пространстве и позволяет получить решение высокой точности, не производя оценку освещенности на каждом проходе (что требуется, например, в  методах PPM и SPPM). Можно сказать, что метод позволяет построить карту освещенности высокой степени детализации и качества. В числе плюсов метода также следует отметить построение структуры с высокой степенью регулярности, что может положительно сказаться на производительности сбора освещенности при реализации на массивно-параллельных системах. 

\noindent В работе \cite{Havran2005eg05} для ускорения финального сбора был предложен подход, называемый обратными фотонными картами (reverse photon mapping). Основная идея метода заключается в том, чтобы вместо фотонов запоминать лучи сбора kd-дереве специального вида. После чего производить трассировку фотонов накапливая результат для лучей сбора. Серьезным недостатком метода является большой объем потребляемой памяти поскольку число лучей сбора, как правило, чрезвычайно велико.


\subsubsection{Комбинированные методы на основе фотонных карт}

\noindent Идея разделения освещенности на независимые компоненты, по-видимому, впервые была высказана в виде алгоритма кэша освещенности \cite{Ward:1988:RTS:378456.378490}. В дальнейшем многие исследователи приходили к тому, что вычислять отдельные компоненты освещенности можно более эффективно за счет априорных знаний о природе компоненты и применении соответствующего алгоритма \cite{BekaertSlusallekCoolsHavranSeidel}, \cite{Hachisuka:2009:SPP:1618452.1618487}, \cite{Vorba:2011}. В работе \cite{Doidge:2012:MMC:2346819.2346840} аппробирована идея комбинирования метод Монте-Карло и прогрессивных фотонных карт с целью использования преимущества обоих методов. Освещение разделялось на 2 части: освещение вызванное каустиками (пути вида $E(S|G)^+DL$ или $E(S|G)^+D(S|G)^+L$) и все остальное освещение. Первую часть освещения предлагается вычислять при помощи метода стохастических прогрессивных фотонных карт, а вторую при помощи стандартной Монте-Карло трассировки путей. Поскольку множества путей из 1 и второй части не имеют пересечений, а их объединение составляет все множество путей, метод приводит к корректному результату. В работе \cite{Hachisuka:2012:PSE:2366145.2366210} предложен улучшенный способ комбинирования результата при помощи взвешивания на основе многократной выборки по значимости.


\noindent В работе \cite{Jensen95importancedriven} был предложен иной способ комбинирования фотонных карт и Монте-Карло трассировки. Вместо того чтобы использовать фотонные карты для вычисления непосредственно освещенности, оценка освещенности на основе фотонной карты использовалась для вычисления плотности вероятности падающего освещения. После чего к методу Монте-Карло применялась выборка по значимости на основе ДФО и посчитанной плотности вероятности падающего освещения. Этот метод является альтернативой MLT, когда падающее освещение неравномерно в высокой степени. 

\subsection{Методы на основе кэша освещенности}

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

\noindent В работах J.Krivanek-а \cite{krivanek05radiance_phd}, \cite{Krivanek:2006:MRI}, \cite{Gautron:2009:PGI:1855086} предложено улучшение контроля ошибки экстраполяции, 
многопроходный адаптивный алгоритм, который, при добавлении новой точки в кэш, сравнивает относительное отклонение освещенности в соседних точках, 
что уменьшает количество артефактов (но приводит к большому числу дополнительных вычислений) и расширение метода кэширования освещенности 
на умеренно-зеркальные ДФО. Падающая освещенность сохранялась при помощи сферических гармоник \cite{SHLighting}. 
Значительным преимуществом кэша падающей освещенности над простым кэшем освещенности является способность ускорять вычисления на 
поверхностях с высокой геометрической детализацией на основе микро-релефа.
 
\noindent В \cite{Jensen:2002:HQR:1281500.1281593} описывается использование процедуры репроекции вместе с кэшем освещенности. Идея репроекции состоит в том, что точки пересечения лучей и поверхностей сцены, которые были вычислены ранее и записаны в кэш, проецируются на полусферу в новой точке поверхности, в которой необходимо вычислить освещение. Недостатком подхода описываемого в \cite{Jensen:2002:HQR:1281500.1281593} является наличие 'пустот' на полусфее вокруг точки, на которую выполняется репроекция. Репроекция получила дальнейшее развитие в работе \cite{Vostryakov:2009:Graphicon}, где устранены основные недостатки кэша падающей освещенности на основе репроекции за счет специального алгоритма расслоения репроецируемых точек (это позволяет корректно реализовать экстраполяцию освещения от нескольких записей кэша с учетом закрытия одних объектов другими) и перетрассирвоки в местах, где были обнаружены пустоты, и представлена высокоэффективная параллельная реализация на основе SSE инструкций центрального процессора.
 
\noindent В работе \cite{Vostryakov:2009:Graphicon} также высказывается идея о том что можно использовать кэш падающей освещенности для оценки плотности вероятности падающего освещения, чтобы ускорить сходимость метода Монте-Карло, если необходимо получать несмещенное решение. 

\noindent В работе \cite{Arikan:2005:FAD} предложено разделять падающеее освещения на дальнее и ближнее с тем, чтобы сохранять только 
дальнее освещение в кэше падающей освещенности. Поскольку освещение от дальних объектов меняется, как правило, значительно более плавно 
чем освещение от ближних, сохранение в кэше только 'дальней' освещенности позволяет значительно снизить число записей кэша. 
Ближнее освещение в работе  \cite{Arikan:2005:FAD} вычислялось при помощи апроксимации, основанной на методе излучательности без учета функции видимости. Это значительно ускоряет вычисление ближнего освещение, однако, приводит к некорректному результату. Тем не менее, развитие идеи разделения падающего освещения на ближнее и дальнее с корректным вычислением ближнего освещения, освещенное в \cite{Vostryakov:2009:Graphicon}, можно считать перспективным направлением.  

\subsection{Методы на основе фильтрации}

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


\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{interpvsfilter}           
  \caption{\label{fig:interpvsfilter}
           Интерполяция и фильтрация. Красные точки обозначают места, в которых значения функции $x(t)$ были вычислены точно. Зеленые точки обозначают значения точечных выборок в методе Монте-Карло при оценки функции $x(t)$ на множестве точек оси $t$. Применив к ним операцию фильтрации можно восстановить с некоторой точностью искомую функцию $x(t)$.}   
\end{figure}

\noindent Фильтрация (рис. \ref{fig:interpvsfilter}) ускоряет расчет отсвещения за счет того же эффекта - переиспользования информации. Однако, вместо того чтобы вычислять значение функции в нескольких точках с высокой точностью, этот метод принимает на вход Монте-Карло выборки. А именно: результирующую яркость путей с некоторой дополнительной информацией - позицией первого пересечения, нормалью к поверхности в этой точке, текстурными координатам, координатами на линзе (если моделируется эффект глубины резкости) и случайными параметрами самого пути. Таким образом, строится многомерное пространство где каждое измерение является некоторым параметром пути. В этом пространстве затем производится фильтрация. Одно из наиболее качественных решений было получено в работе \cite{Sen:2012:FNR:2167076.2167083}, где фильтрация производится с учетом зависимости между случайными величинами в методе трассировки лучей Монте-Карло и значениями в полученной выборке.

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

\section{Резюме по методам вычисления глобальной освещенности}

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

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

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

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


\section{Краткое введение в программирование GPU}

Когда первые видеокарты только появились в продаже, они представляли собой достаточно простые (по сравнению с центральным процессором) узкоспециализированные устройства, 
предназначенные для того чтобы снять с процессора нагрузку по визуализации двухмерных данных. С развитием игровой индустрии и распространением трехмерных игр появился спрос на ускорение 3D графики. 
Со времени создания компанией 3Dfx первых видеокарт Voodoo, (1996 г.) и вплоть до 2001 года в GPU был реализован только фиксированный набор операций над входными данными.
У программистов не было никакого выбора в алгоритме визуализации, и для повышения гибкости появились шейдеры — небольшие программы, выполняющиеся видеокартой для каждой вершины либо для каждого пиксела. 
В их задачи входили преобразования над вершинами и затенение — рассчет локального освещения в точке, например по модели Фонга. 
Начиная с этого момента видеокарты перестали называть просто ускорителями и появился термин графические процессоры (Graphics Processing Unit, GPU) и на сегодняшний день графические процессоры
представляют собой чрезвычайно сложные и высоко-производительные системы с широкими возможностями. 
Опуская исторические моменты эволюции GPU и их программно-аппаратной модели, рассмотрим далее псовременные GPU и их возможности. 

\subsection{Аппаратная модель современных GPU}

Чтобы понять причину возникновения программно-аппаратной модели GPU, необходимо обратиться к современным методам построения центральных процессоров в частности и вычислительных систем вообще.
Основы построения современных вычислительных устройств достаточно подробно изложены в \cite{Hennessy:2011:CAF:1999263}. 
Рассмотрим далее одну из наиболее значительных трудностей при создании конвейеных вычислительных устройств. 
Допустим, имеется простой конвейерный процессор со следующими стадиями в конвейера:

\begin{enumerate}
\item F - Fetch (Выборка команды из памяти).
\item D - Decode (Декодирование команды и выборка данных из регистрового файла).
\item X1,X2,X3 - eXecute; (3 стадии выполнения команды - допустим мы имеем дело с простой реализацией плавающей точки).
\item W - Write Back (запись данных обратно в регистровый файл). 
\end{enumerate}

\noindent При том что латентность (время выполнения в конвейере) каждой инструкции в такой системе 6 тактов, пропускная способность - 1 инструкция за такт.
Это хорошо известное свойство конвейерных систем. Однако, конвейер позволяет достичь производительности в 1 операцию за такт только при условии независимости команд друг от друга.
Рассмотрим ситуацию (листинг \ref{asm1}), в которой условие независимости нарушается и построим для такой ситуации диаграму прохождения по конвейеру (листинг \ref{asm2}).
Будем считать что в процессоре реализована передача операндов непосредственно на X1 с выходного регистра стадии X3 (bypassing или forwarding \cite{Hennessy:2011:CAF:1999263}).
Поскольку умножение ($fmul$) выполняется в конвейере 3 такта (по условию), следующая команда ($fadd$) вынуждена ожидать значения операнда R1 на стадии декодирования, останавливая весь конвейер.

\begin{lstlisting}[label=asm1, caption=Код вызывающий задержки в конвейере]
  fmul R1, R2, R3
  fadd R0, R1, R4
  fsub R6, R0, R5
\end{lstlisting}


\begin{lstlisting}[label=asm2, caption=Диаграмма прохождения команд CPU по конвейеру c задержками]
  такт  0  1  2  3  4  5  6  7  8  9  10 11
  fmul  F  D  X1 X2 X3 W
  fadd     F  D  D  D  X1 X2 X3 W
  fsub        F  F  F  D  D  D  X1 X2 X3 W
\end{lstlisting}


\noindent Нетрудно посчитать, что в случае 3-ёх стадийного конвейера и непрерывного потока зависимых инструкций производительность упадет в 3 раза.
Современные центральные процессоры могут иметь от 10 до 20 стадий в ковейере и, как следствие, падение производительности будет значительно более существенным. 
Проблема задержек в конвейере актуальна не только для вычислительных инструкций но также и для операций с памятю (причем даже в большей степени).
Поиск в древовидных структурах на современных центральных процессорах в действительности мало-эффективен именно по этой причине.
Каждый кэш промах может стоить сотни тактов простоя. А при большом размере дерева кэш промахи практически неизбежны.
Важно отметить, что алгоритмы вычисления глобальной освещенности рассмотренные в этой и предыдущих главах используют такой поиск чрезвычайно часто и именно он, как правило, является узким местом всей системы.

\noindent Рассмотрение способов решения данной проблемы задержки конвейера и ликвидации простоев в современных центральных процессорах выходит за рамки литературного обзора данной работы.
Однако их необходимо перечислить. Это такие механизмы как внеочередное выполнение команд (out of order execution), кэши (в том числе и неблокируюшие - out of order memory system), явный паралеллизм инструкций (EPIC) \ref{fig:cpu_vs_gppu_caches}.
Эти механизмы позволяют в некоторой степени амортизировать или решить проблемы зависимости по данным. Однако, они как правило чрезвычайно требовательны к ресурсам кристалла (кроме EPIC) и крайне сложны в реализации.

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{cpu_vs_gppu_caches}           
  \caption{\label{fig:cpu_vs_gppu_caches}
            Сокрытие латентности памяти на CPU и GPU достигается при помощи разных механизмов \cite{Aaron:Houston:2009:cource}. }    
\end{figure}

Известно, что GPU выполняют потоки в SIMD группах, называемых в термином warp \cite{Nvi13b}. Однако, эта особенность не столь важна для понимания механизма, позволяющего избавиться от задержек в GPU.
Оставим пока эту особенность в стороне и рассмотрим как проблема задержек и простоев решается в GPU. В графических процессорах мультипроцессор (вычислительное ядро) всегда выполняет множество потоков одновременно 
(вернее групп потоков, но для рассматриваемого модельного случая будем считать пока, что размер группы - 1 поток).
То есть, можно условно считать что процессор имеет достаточно активных потоков, готовых к выполнению. Будем считать, что в рассматриваемом модельном случае GPU выполняет всего 3 потока одновременно (3 группы по 1 потоку).
При этом, будем помещать в конвейер команды по очереди из каждого потока - сначала из первого, затем из второго, потом из третьго и, наконец, снова из первого. 
Тогда диаграмма из листинга \ref{asm2} превратиться в диаграмму \ref{asm4} и задержки будут устранены.  

\begin{lstlisting}[label=asm4, caption=Диаграмма прохождения команд GPU по конвейеру без задержек]
  такт     0  1  2  3  4  5  6  7  8  9  10 11 12 13
  fmul(1)  F  D  X1 X2 X3 W
  fmul(2)     F  D  X1 X2 X3 W
  fmul(3)        F  D  X1 X2 X3 W
  fadd(1)           F  D  X1 X2 X3 W
  fadd(2)              F  D  X1 X2 X3 W
  fadd(3)                 F  D  X1 X2 X3 W
  fsub(1)                    F  D  X1 X2 X3 W
  fsub(2)                       F  D  X1 X2 X3 W
  fsub(3)                          F  D  X1 X2 X3 W
\end{lstlisting}



Один мультипроцессор современных коммерческих GPU способен выполнять до 2048 потоков одновременно (64 группы warp по 32 потока) \cite{Nvi13b}. 
За счет наличия такого большого числа активных потоков, GPU способны скрывать при помощи рассмотренного механизма не только латентность вычислительных операций, но и латентность доступа к памяти (рисунок \ref{fig:cpu_vs_gppu_caches}).
При этом, хотя кэши в GPU существуют, они, во-первых не являются столь необходимыми как в случае CPU, а во-вторых, работают по иному принципу, стараясь в большей степени разделять данные между потоками, 
используя пространственную когерентность данных (2D и 3D тестур) и алгоритмов. 


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{cpu_vs_gppu_arch}           
  \caption{\label{fig:cpu_vs_gppu_arch}
            Различия в аппаратной модели вычислительных ядер центральных (CPU, слева) и графических (GPU, справа) процессоров \cite{Aaron:Houston:2009:cource}. }    
\end{figure}


Теперь рассмотрим организацию типичного мультипроцессора GPU (рисунок \ref{fig:cpu_vs_gppu_arch}). В каждый определенный момент времени мультипроцессор способен выполнять на всех ядрах только одну инструкцию. 
При этом, Арифметико Логические Устройства (АЛУ) мультипроцессоров могут быть как скалярными так и векторными. Необходимость в явной поддержке векторных операций и паралеллизма на уровне инструкций отпадает, поскольку мультипроцессор имеет достаточно 
независимых команд благодаря большому числу независимых потоков. Показательным примером является отказ от VLIW архитектуры в новейших графических процессорах компании AMD c архитектурой GNC. 

\noindent Вернемся к вопросу об организации вычислений группами. Объединение потоков в группы (warp), обозначенное ранее, позволяет дополнительно экономить ресурсы кристалла, отводя больше места на функциональные устройства.
Однако это далеко не главная причина объединения потоков в группы. Гораздо более важно, что такое объединение позволяет объединять много запросов в память в один большой синхронный запрос (coalescing \cite{Nvi13b}).
Такая модель позволяет работать подсистеме памяти значительно более эффективно \cite{Hennessy:2011:CAF:1999263}.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{gpu_branches_tf}           
  \caption{\label{fig:gpu_branches_tf}
            Выполнение кода с ветвлением на GPU \cite{Aaron:Houston:2009:cource}. }    
\end{figure}

В каждый определенный момент времени мультипроцессор способен выполнять на всех ядрах только одну инструкцию, возникает вопрос - как быть с ветвлениями? 
Ведь если в коде программы встречается ветвление, то инструкции будут уже разные. Здесь применяется следующее решение (рисунок \ref{fig:gpu_branches_tf}).
При расхождении на ветвлениях сначала маскируются потоки, неактивные внутри ветки if (рисунок \ref{fig:gpu_branches_tf}) и происходит выполнение этой ветки.
После чего, замаскированные потоки становятся активными и маскируются потоки, для которых была выполнена ветка if.
Таким образом, мультипроцессор выполнит обе ветки кода. Эта особенность является достаточно важной, поскольку потери на ветвления могут быть значительными (в среднем от 2 до 4 раз). 
При создании программ для GPU, эту особенность необходимо учитывать и по-возможности не допускать расхождения потоков из 1 группы warp по различным участкам кода.

Еще одной важной особенностью GPU является то, что большая часть локальных переменных программы хранится в регистрах, а сами регистровые файлы гораздо больше чем регистровые файлы центральных процессоров. 
Тем не менее они имеют ограниченный объем и чем меньше локальных переменных использует программа, тем больше потоков может одновременно выполнять мультипроцессор GPU \cite{Nvi13b}.
А это число, как было описано выше, напрямую влияет на способность мультипроцессора скрывать латентность вычислительных операций и, что гораздо более важно, латентность подсистемы памяти.
В связи с этой особенностью, важно избегать объемного кода с большой степенью вложенности процедур и большим числом локальных переменных, поскольку компилятор для GPU старается встраивать (inline) код, 
заменяя вызов процедуры на тело самой процедуры. Эффективность встраивания падает с ростом глубины вызовов и сложности программного кода. 


\noindent Отдельно следует сказать и о некоторых непрограммируемых блоках GPU, имеющих важное значение в программировании общего назначения. 
Это такие блоки как:

\begin{enumerate}

\item ROP - блоки графического процессора, выполняющие Z-буферизацию, полноэкранное сглаживание и запись итогового изображения (с алфа-смешиванием или без) в кадровый буфер в видеопамяти. 
      Используется в программировании общего назначения для реализации атомарных операций \cite{patent:8411103}.
      Справедливости ради следует сказать, что современные центральные процессоры, как правило, также имеют специальную логику для реализации атомарных операций.
      
\item Разделяемая память и кэш второго уровня L2. Изначально предназначались для передачи данных между стадиями графического конвейера и кэширования текстур, 
      но могут быть использованы для сообщения данными между потоками одного блока (обычно группа из 64-512 потоков) \cite{Nvi13b}.
      
\item Текстурные блоки, позволяющие осуществлять не просто чтение данных, но и фильтрацию, преобразование форматов и декомпрессию текстур на аппаратном уровне на лету. 
      К сожалению, в API программирования GPU общего назначения CUDA и OpenCL возможно использование только ограниченного подмножества их функциональности.   

\item GMU - блок, появившийся в архитектуре GPU компании Nvidia - Kepler. 
      Этот блок ответственнен за механизм создания работы на GPU (GPU Work Creation, GWC) и позволяют осуществлять вызов ядер из самих ядер \cite{NvGK110}.
      
\item RDMA (Remote Direct Memory Access) блоки - состявляющая GPU, осуществляющая возможность прямого доступа (без участия CPU) к памяти GPU из других устройств \cite{NvGK110}.
      Как правило, используются для передачи данных между несколькими GPU.   
       
\end{enumerate}

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

\subsection{Программная модель современных GPU}

Программная модель cовременных GPU основана на концепции паралеллизма по данным, называемым иначе - ОПМД (Одна Программа Множественные Данные) или SIMT (Single Program Multiple Data).
Центральный процессор (а вслучае процессоров Kepler и сам GPU) асинхронно создает очереди запросов на выполнение небольших программ, называемых ядрами (kernel \cite{Nvi13b}).
Запросы извлекаются из каждой очереди по одному (но возможно параллельное извлечение из нескольких очередей - конкурентное выполнение ядер \cite{Nvi13b}) и выполняются на GPU.

Причем, GPU имеют 2 различных режима работы, в которых они могут выполнять:

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

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

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{opengl43}           
  \caption{\label{fig:opengl43}
            Современный стандарт OpenGL 4.3; В левой части изображения структура API на основе графического конвейера, в правой - структура API для вычислений общего назначения. }    
\end{figure}


\noindent Соответственно этим 2 режимам, существуют 2 программные модели GPU и 4 распространенных API (Application Programming Interface) - (CUDA и OpenCL) для вычислений общего назначения и
(OpenGL и DirectX) для программирования графического конвейера.  
В настоящий момент графические API, за иcключением специфических для производителя GPU вещей (например создание работы непосредственно на GPU в CUDA), являются значительно более мощными и в 
смысле функциональности использования возможностей GPU могут считаться надмножеством API вычислений общего назначения. 
На рисунке \ref{fig:opengl43} обозначена структура графического API OpenGL 4.3, в котором с правой стороны находятся вычислительные шейдеры (Compute Shaders), предназначенные
для реализации вычислений общего назначения средствами API OpenGL 4.3. Аналогично, API DirectX имеет свои вычислительные шейдеры. 

\noindent Рассмотрим далее принципы использованием таких API как CUDA и OpenCL. 
Каждый запрос на выполнения вычислительного ядра запусает множество копий программы-ядра на GPU. 
При этом, потоки имеют трехмерный (в общем случае) индекс и объеденены в группы, называемые блоками потоков (в OpenCL такой блок называется рабочей группой).
Совокупность всех блоков называется сеткой (grid).

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{cudathreads}           
  \caption{\label{fig:cudathreads}
            Организация потоков в группы в двух измерениях. }    
\end{figure}


\noindent На рисунке \ref{fig:cudathreads} изображена такая организация с использованием двумерного индекса; третье измерение не используется.
Каждый блок потоков гарантированно выполняется на 1 мультипроцессоре (хотя 1 мультипроцессор может выполнять несколько блоков), за счет чего возможна передача данных между потоками
в пределах одного блока через быструю разделяемую память, физически расположенную близко к мультипроцессору. 
Двумерная или трехмерная организация, в свою очередь, позволяет сделать обработку двумерных и трехмерных данных более естественной.
Процесс вычислений можно представить в виде следующей последовательности действий:

\begin{enumerate}

\item Выделить память для входных и выходных данных на GPU.

\item Скопировать входные данные из памяти CPU в память GPU.

\item Запустить 1 или более вычислительных ядер. 

\item Дождаться окончания выполнения всех ядер путем явной или неявной синхронизации (происходит при вызове функций копирования данных на сторону CPU).

\item Скопировать выходные данные обратно на сторону CPU либо предать их в область видимости графических API для дальнейшей визуализации.
      Повторять шаги 3-5 до тех пор пока необходимо что-то вычислять.
      
\item Если работа окончена, освободить ранее выделенную память.

\end{enumerate}


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

\noindent Более сложные алгоритмы предполагают синхронизацию/взаимодействие между группой потоков и/или всеми потоками сетки и должны быть реализованы в нескольких ядрах. 
Такие алгоритмы могут быть реализованы при помощи набора более простых операций - редукции, уплотнения, сортировки. 
Эти операции допускают эффективную программную реализацию на GPU с использованием разделяемой памяти и последовательным запуском ядер, каждое из которых реализует часть алгоритма.
Достаточно простым примером здесь является реализация поиска минимума некоторой функции $f$ на некотором множестве значений $data$. 
Сначала каждый поток вычисляет свое значение функции $f$ для заранее определенных входных данных $data[i]$.
Где $i$ - идентификатор потока. $data$ - массив данных. После чего минимум среди всех значений может быть найден при помощи операции редукции.
Следует отметить, что частичная редукция (внутри одного блока потоков) при этом может совершаться в том же ядре, которое вычисляло значение функции $f(data[i])$. 
Еще один пример - построение дерева. Если алгоритм построения дерева удается свести к сортировке (причем, очень желательно - к сортировке 32 или 16 битных чисел, для которых хорошо работает поразрядная сортировка), 
тогда построение такого дерева может быть эффективно реализовано на GPU.

\noindent В соответствие с рассмотренной выше моделью, для того чтобы алгоритм мог быть эффективно реализован на GPU он должен удовлетворять следующим условиям:
 
\begin{enumerate}

\item На каждом шаге алгоритм должен иметь достаточно большое число относительно-простых независимых операций. Каждая такая операция будет назначена на отдельный поток.
      Под термином 'достаточно' здесь понимается столько операций, сколько необходимо для полной загрузки GPU работой. 
      Под 'отонсительной простотой' понимается отсутствие рекурсии, ограниченное число локальных переменных и небольшой уровень вложенности функций, 
      при котором компилятор способен эффективно распределить локальные переменные по регистрам GPU.      
      
\item Зависимости по данным между различными операциями необходимо сводить к минимому. 
      Хорошим случаем является такой, при котором различным потокам нужно обмениваться данными только локально - внутри своей группы.
      В случае когда это сделать невозможно, алгоритм должен сводиться к набору примитивных операций - редукции, уплотнения, сортировки.      

\item Если значимая часть вычислений выполняется на CPU, желательно реализовывать алгоритм в виде конвейера, чтобы операции копирования данных и 
      вычисления на CPU перекрывались по времени с вычислениями на GPU.
      

\end{enumerate}

\subsubsection{API программирования GPU}

\noindent Хотя графические API являются более мощными, трудоемкость разработки программ с их использованием значительно выше по целому ряду причин:

\begin{enumerate}

\item Графические API, как правило, работают с GPU на более низком уровне, чем API для вычислений общего назначения. 
      Например, при создания буфера или текстуры в DirectX программист обязан при помощи флагов самостоятельно указывать, как данный ресурс будет использован - 
      для записи, чтения, копирования c CPU на GPU или c GPU на CPU. Если нужный флаг (например на запись) не был установлен, при попытки записи в такой ресурс будет сгенерирована ошибка.
      Во многих случаях такой детальный контроль над ресурсами не нужен, а лишь вынуждает программиста тратить время на работу с несущественными деталями. 

\item Поскольку графические API призваны в первую очередь обеспечивать доступ ко все функциональности GPU, 
      они значительно более объемные. При этом, DirectX и OpenGL имеют различные недостатки:
      
      \begin{enumerate}
      
        \item Помимо того, что DirectX поддерживается только в операционной системе Windows, его главный недостаток - ограниченное время жизни.
              DirectX всегда создавался под определенное поколение графических процессоров. При выходе нового поколения GPU появлялась новая версия DirectX.
              При этом, каждая версия DirectX жестко фиксирует набор используемых функциональных возможностей GPU. 
              Чтобы использовать новую функциональностью GPU, необходимо использовать другую версию API. 
              Это выливается в необходимость переписывания исходного кода под новые версии API, которые могут различаться достаточно сильно (особенно DirectX 9 и 10 версии).
              При этом, если в GPU появилась какая-то новая функциональность, а новая версия DirectX еще не вышла - использовать эту возможность разработчик просто не сможет.
              
        \item Основной недостаток OpenGL - слабая поддержка на уровне драйвера, особенно касательно компиляторов шейдерных программ.
              Поскольку шейдерные программы компилируются на машине пользователя, при разработке необходмо не только строго следовать стандарту, 
              но еще и использовть минимальное подмножество языка GLSL (язык шейдеров в OpenGL) и функциональности самого OpenGL. 
              Основная причина такой необходимости заключается в том, что разные производители (и даже разные версии драйверов одного и того же производителя!) немного по-разному трактуют стандарт и не совсем строго ему следуют.
              Эта проблема знакома всем кто разрабатывал программы на C++ под различные компиляторы. 
              Однако, если в случае компиляторов C++ проблемы несовместимости могут быть устранены на этапе сборки программы, то в случае шейдеров OpenGL они возникнут на машине пользователя. 
              Справедливости ради следует сказать, что такая проблема характерна и для OpenCL. Однако, вследствие меньшей функциональности самого API, для OpenCL это не так значимо.
                                
      \end{enumerate}
      
\end{enumerate}

\noindent Корнем проблем всех API программирования GPU является противоречивость требований, предъявляемых к ним - переносимость с одной стороны и доступ к широкой 
функциональности GPU + достижение высокой эффективности - с другой. API для вычислений общего назначения CUDA является проприетарным API от компании NVIDIA. 
Среди его достоинств необходимо отметить:

\begin{enumerate}

\item Стабильную поддержку на уровне драйвера и заранее известный набор аппаратуры, на которой может быть запущена программа. 
      Это упрощает разработку и оптимизацию программного кода.

\item Статическую компиляцию ядер (шейдеров) в бинарное представление. Компиляция осуществляется на машине разработчика.

\item Гомогенность кода разрабатываемой системы. CUDA компилятор использует расщирение языка C++ и поддерживает значительную часть стандарта С++.      

\item Наличие значительного числа библиотек.
    
    \begin{enumerate}
      \item thrust - библиотека, реализующая базовые примитивы параллельного программирования (Параллельную сортировку, уплотнение, вычисление префиксной суммы).
      \item CUBLAS, MAGMA, cuSPARSE, CUSP, HiPLAR - реализация базовых операций линейной алгебры.
      \item CURAND - реализация генераторов псевдо-случайных и квази-случайных чисел.
      \item cuFFT  - реализация быстрого преобразования фурье
      \item NPP    - библиотека для обработки изображений.
      \item Другие библиотеки, предназначенные для специализированных задач \cite{CUDALibs}. 
    \end{enumerate}

\end{enumerate}

\noindent При этом, существование таких проектов как PGI CUDA-x86 \cite{PGIx86} и GPUOcelot \cite{Farooqui:2011:FDI:1964179.1964192} позволяет говорить о переносимости в перспективе.


\section{Алгоритмы на массивно-параллельных процессорах}

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

\begin{enumerate}
\item Исследования о том, как эффективно реализовать некоторый известный алгоритм на GPU, если принципиальные трудности по его 
реализации отсутствуют, однако, по каким-либо причинам, оцениваемая производительность не достигается. 
Будем условно называть данную группу исследований 'исследованиями производительности алгоритмов на GPU'.

\item Исследования о том, как реализовать на массино-параллельной системе алгоритм, изначально спроектированный для однопроцессорной системы. 
Такой алгоритм в силу своей природы не может быть напрямую реализован на массивно-параллельной системе и требует переработки в самих основах. 
Будем называть эту группу 'исследованиями по распараллеливанию алгоритмов на GPU' или 'исследованиями по созданию высоко-паралельных алгоритмов'.
\end{enumerate}

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

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

\noindent Ранние работы по реализации трассировки лучей на GPU исследовали проблему эффективной реализации без-стекового поиска в ускоряющих структурах в условиях ограниченных возможностей графических API DirectX 9 и OpenGL2 и их аппаратных реализаций \cite{Foley:2005:KAS:1071866.1071869}, \cite{Horn:2007:IKT:1230100.1230129}, \cite{guenther:07:BVHonGPU}. 

\noindent Можно легко ошибиться в предположении, что с появлением новых технологий и GPU с более широкими возможностями и ресурсами задача реализации трассировки лучей на современных GPU тривиальна, поскольку достаточно лишь назначить каждому потоку свой луч и портировать существующий CPU код. В действительности это не так. Портированный 'как есть' код будет обладать чрезвычайно низкой эффективностью вследсвие особенностей выполнения кода на GPU, а именно \cite{Frolov:2009:Graphicon}:
\begin{enumerate}

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

\item Эффективность подсистемы памяти также зависит от числа активных потоков на мультипроцессоре, поскольку механизм 
      сокрытия латентности доступа к памяти основан на выполнении большого числа потоков параллельно.
  
\end{enumerate} 

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


\noindent В работе \cite{Frolov:2009:Graphicon} выделены 2 способа решения, называемые "убер-ядро" (uber-kernel или mega-kernel, и "разделенное ядро" (separate kernel). 
          При реализации способа "убер-ядро" исходный код логически разбивается на состояниея конечного автомата, в котором каждое состояние - некоторая отдельная, достаточно тяжелая часть кода. Во время выполнения, происходит периодическая смена состояний и переход между соответствующими ветками кода, сохраняя некоторые важные данные в разделяемой или локальной памяти. 
          Суть такого подхода в том, что он позволяет использовать одни и те же регистры для разных переменных. 


\noindent Однако в \cite{Frolov:2009:Graphicon} отмечаются недостатки данного подхода, среди которых неэффективное использование ресурсов, потери производительности на ветвлениях (когда одна часть группы warp находится в одном состоянии, а другая - в другом), сложность поддержки и отладки.




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

\noindent Проблемы подхода убер-ядра обсуждаются позднее в \cite{Bogolepov:2013:phdthesis} и \cite{Laine2013hpg}. В работах \cite{Frolov:2009:Graphicon} и \cite{Bogolepov:2013:phdthesis} предлагается использовать подход 'разделенного ядра', формируя и настраивая конвейер трассировки лучей за счет замещения и/или вставки ядер на различных этапах конвейера. Этот способ предполагает реализацию различных частей в отдельных ядрах и организацию взаимодействия между ними путем сохранения данных в глобальной памяти. Работа \cite{Laine2013hpg} развивает данный подход на основе очередей на выполнение работы (каждая очередь запускает на выполнение только одно конкретное ядро), сортировки потоков по номеру ядра и механизма создания работы на GPU.  



\subsection{Монте-Карло Трассировка путей на GPU}

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

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

\noindent В работе \cite{Novak:2009:masterthesis} представлено простое и эффективное решение данной проблемы, называемое русской рулеткой на warp (per warp russian roulette). Основная идея данного метода заключается в том, чтобы применять метод русской рулетки не к каждому потоку отдельно (как было рассмотрено в 1 главе), а сразу ко всей группе потоков warp. Решение о терминировании, таким образом, принимается сразу для всей группы, что исключает возникновение групп warp с малым числом активных потоков при терменировании лучей посредством русской рулетки.

\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{perwarprr}           
  \caption{\label{fig:perwarprr}
           Применение на различных сценах русской рулетки на warp из работы \cite{Novak:2009:masterthesis}. Красные столбцы отражают число потоков, находящихся на GPU. Зеленые столбцы отражают число активных, выполняющих полезную работу, потоков. Вернхие 3 изображения сделаны для обыкновенной русской рулетки на поток. Нижние 3 - для русской рулетки на группу потоков warp.  }   
\end{figure}

\noindent Для решения проблемы неоднородной загрузки GPU из-за различной глубины терминирования потоков, в работах \cite{Novak:2009:masterthesis} и \cite{Bogolepov:2013:phdthesis} исследовался метод повышения эффективности использования GPU, называемый регенерацией путей. Сообщается об ускорении от 20 до 100\% для случаев, требующих учет большого числа переотражений. Метод был впервые предложен в \cite{vanAntwerpen:2011:ISE:2018323.2018330} и иследован на стандартной трассировке путей, двунаправленной трассировке путей и алгоритме MLT. В работе \cite{vanAntwerpen:2011:ISE:2018323.2018330} сообщается об ускорении от 20 до 50\%. Метод работает следующим образом:

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

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

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

\noindent В работе \cite{sadeghi2009cpt} описывается подход, называемый когерентной трассировкой путей. 
          Идея данного подходя заключается в том, чтобы использовать одно состояние случайного герератора для всех лучей (в группе варп, блоке или вообще всех на экране).


\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{coherentpt}           
  \caption{\label{fig:coherentpt}
           Стандартная трассировка путей (слева). Когерентная трассировка путей (справа). \cite{sadeghi2009cpt}.  }   
\end{figure}

\noindent Данный подход является эффективным но меняет характер артефактов, превращая шум в более заметные для глаза артефакты с резкими границами. Для того чтобы убрать границы в \cite{sadeghi2009cpt} применяется случайное перемешивание к схеме нумерования пикселей. 

                
\subsection{Параллельный Кэш освещенности}

\noindent Реализация параллельного кэша освещенности даже на центральном процессоре сопряжена с рядом трудностей. 
\begin{enumerate}
\item Первая проблема (проблема синхронизации) связана с обеспечением консистентности данных кэша - каждый поток должен видеть один и тот же экземпляр данных кэша.
\item Вторая проблема (последовательная природа) связана с тем, что алгоритм вычисляет расположение записей кэша и их радиусы валидности 
на основании эвристических методов, опираясь на то, что точки вставляются в кэш последовательно, одна за одной \cite{Gautron:2009:PGI:1855086}. 
\end{enumerate}
 
\noindent Один из простых способов решения проблемы синхронизации, используемый в [mitsuba, проверить] и \cite{Vostryakov:2009:Graphicon} 
заключается в том, чтобы не синхронизировать данные кэша между потоками вообще.
Первый недостаток этого подхода состоит в том, что на изображении могут появиться видимые границы между областями, 
которые были использованы разными потоками т.к. каждый поток будет иметь доступ к своему уникальному подмножеству записей кэша - и, 
как следствие, в одной той же точке трехмерного пространства для разных потоков результат интерполяции будет различен.
Второй недостаток заключается в том, что при увеличении числа потоков будет возрастать объём избыточных вычислений, 
поскольку каждый поток будет вычислять полностью независимый набор записей, а их перекрытие никак не контролируется.

\noindent Следующий способ решения проблемы снхронизации - использование разделяемой структуры данных в том или ином виде.
Здесь возможны различные варианты - синхронизация при вставке единичной записи \cite{Ward:1994:RLS:192161.192286}, 
синхронизация на основе слияния групп \cite{conf/acpc/KoholkaMG99} и более сложные подходы - \cite{conf/egpgv/DebattistaSC06}, \cite{conf/egpgv/DublaDSC09}.

\noindent В работе \cite{conf/egpgv/DebattistaSC06} предложен подход, использующий пространственную когерентность кэша освещенности для реализации
алгоритма на кластерных системах. Все узлы кластера разбивались на 2 подмножества - узлы, занимающиеся выислением записей кэша и узлы, отвечающие за построение изображение в целом.  
За счет использование такого разделения и высокой пространственной когерентности операций чтения из кэша, авторам \cite{conf/egpgv/DebattistaSC06} 
удалось сохранять частоту обновления кэша на высоком уровне бес существенных потерь в производительности.  
 
\noindent Работы \cite{Ward:1994:RLS:192161.192286}, \cite{conf/acpc/KoholkaMG99} и \cite{conf/egpgv/DebattistaSC06} используют для синхронизации критические секции.
Корректная реализация с использованием критических секций блокирует как запись в кэш, так и чтение из него. 
Это негативно сказывается на масштабируемости реализации при увеличении числа параллельно исполняющихся потоков \cite{conf/egpgv/DublaDSC09}.
В работе \cite{conf/egpgv/DublaDSC09} предлагается подход к построению разделяемой структуры кэша освещенности на 
основе атомарных операций современных CPU. За счет использования атомарных операций в работе \cite{conf/egpgv/DublaDSC09} удается избежать
блокирования при чтении из кэша освещенности. 

\noindent Наконец, еще один способ решения проблемы синхронизации заключается в том, чтобы  отделить процесс построения кэша освещенности от 
процесса его использования и разделить, таким образом, процесс построения изображения на 2 стадии:

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

\noindent Идея такого подхода обсуждается в \cite{Krivanek:2006:MRI}, \cite{Gautron:2009:PGI:1855086} и \cite{Pharr:2010:PBR:1854996}.
Причем, в указанных работах двухстадийный алгоритм используется не для параллельной реализации, а для повышения точности интерполяции (рисунок \ref{fig:ic_two_pass}).
Ранее была рассмотрена ситуация, при которой разные потоки, имея различные подмножества записей кэша в одной и той же точке пространства
будут давать различные результаты интерполяции. 

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{ic_two_pass}           
  \caption{\label{fig:ic_two_pass}
           Сравнение различных подходов построения изображения с применением кэша освещенности из работы \cite{Gautron:2009:PGI:1855086}.
           Под 'размытием' (случай d) в данном случае понимается небольшое увеличение радиусов валидности для всех записей (в 1.6 раза). }   
\end{figure}


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


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

\begin{enumerate}
\item Алгоритм создания самого множества записей кэша не предполагает одновременную вставку многих точек в кэш.
      Причем, здесь существуют 2 проблемы:
      \begin{enumerate}
      \item Проблема остановки алгоритма построения кэша. Сам процесс расстановки записей реализует ленивую схему вычислений.
            Оригинальный алгоритм (\ref{ic_alg}) вычисляет новые точки (или записи) кэша по запросу.
            Поэтому в нем проблема остановки не существовала. При необходимости новая точка (запись) всегда может быть вычислена. 
            В двух-проходной схеме это не так. Процесс создания кэша освещенности должен уметь останавливаться.  
            Реализация из \cite{Krivanek:2006:MRI} использует последователный алгоритм (называемый 'Full convergence for a set of shading points'). 
            На 1 проходе алгоритм запоминает в kd-дереве все точки, в которых потребуется осуществляь выборку из кэша (shading points).
            При вставке новой записи в кэш освещенности, из дерева удаляются точки (shading points), покрытые областью, на которую распространяется влияние
            вставляемой записи. 
            В \cite{Gautron:2009:PGI:1855086} предлагается использовать упрощение рассмотренного выше алгоритма, помечая пикселы изображения, для которых не встретилось
            высоких значений оценки ошибки интерполяции. Когда все пикселы помечены, можно остановить алгоритм.                    
            Хотя параллельная реализация обоих решений теоретически возможна, эта проблема требует исследования, поскольку
            здесь возникает все та же необходимость синхронизации структур данных, хранящих shading points или маркеры пикселей. 
            
      \item Проблема выбора оптимального расположения записей кэша. 
            В рассмотренных ранее подходах, при увеличении числа одновременно-вставляемых в кэш записей неизбежно 
            падение эффективности решения в целом. Это связано с тем, что при большом числе одновременно вставляемых точек (записей),
            некоторое число запросов на вычисление освещенности могут возникнуть в близко-расположенных областях пространства.
            Такие запросы являются избыточными, поскольку в этих областях пространства последовательный алгоритм
            обошелся бы всего одной точкой. 
      \end{enumerate}

\item Алгоритм вычисления радиусов валидности. Эвристика Neighbour Clamping \cite{Gautron:2009:PGI:1855086}, значительно повышающая 
точность решения на основе Кэша освещености, предполагает последовательную вставку записей в кэш. 
То же верно и об алгоритме контроля плотности, описываемом в разделе 5.1 работы \cite{Krivanek:2006:MRI}.
Отсутствие учета соседства точек при параллельной вставке большого числа независимых записей приведет к тому, что радиусы валидности записей будут слишком большими. 
Как будет показано в главе 3 - это, в свою очередь, ведет к высокой степени перекрытия записей и медленной интерполяции.

\end{enumerate}

\subsubsection{Кэш освещенности на GPU}

Впервый GPU был использован при реализации кэша освещенности в работе \cite{conf/rt/GautronKBP05} и в дальнейшем обсуждается в \cite{Gautron:2009:PGI:1855086}.
Рассматриваемый в этих работах подход позволяет заменить стандартный алгоритм интерполяции на растеризацию сплатов.
Это позволяет избежать поиска в дереве на GPU при финальной визуализации сцены, что было особенно для ранних GPU, не имеющих поддержку циклов.  
В работах \cite{conf/rt/GautronKBP05} и \cite{Gautron:2009:PGI:1855086} также предлагается использовать растеризацию на GPU при 
рассчете освещенности в точках кэша. Метод, обсуждаемый в этих работах имеет следующие ограничения:

\begin{enumerate}

\item Интерполяция при помощи растеризации сплатов позволяет визуализировать вторичное освещение только на первично-видимых поверхностях.
Хотя существуют методы визуализации отражений на основе растеризации (и в этом случае, растеризация сплатов также применима), 
такие методы работают для специальных случаев и/или визуализируют отражения не совсем корректно \cite{CubeMapGL}.
Поэтому, при точной визуализации и/или визуализации на основе трассировки лучей, данный метод неприменим. 
                                   
\item При использовании растеризации для вычисления непрямого освещения, значения освещенности сохраняемые в кэше, 
несут информацию только о первом переотражении. Такое решение неприемлимо, если необходимо учитывать большее число переотражений света. 
Хотя этот недостаток можно побороть при использовани фотонных карт в сочетании с картами светимости (iradiance maps), остается открытым вопрос об
эффективности реализации такого подхода на современных GPU. 
Оценим эффективность этого подхода с учетом поледних данных о производительности трассироки лучей и растеризаци на GPU.

\begin{enumerate}

\item Современные GPU могут обеспечить производительность обработку порядка 100-300 миллионов (M) лучей в секунду \cite{Aila:2009:UER:1572769.1572792}. 
С учетом необходимости трассировки лучей на большую чем 1 глубину и/или необходимости выполнения финального сбора, эти числа нужно разделить как минимум на 2. 
Однако, с учетом того что лучи из точек кэша можно выпускать когерентными группами и что большинство из них достаточно быстро находит 
пересечение с ближайшими объектами, слишком сильно занижать производительность не следует. Поэтому примем за среднюю производительность 50M монте-карло выборок в секунду. 

\item Современные компьютерные игры способны растеризовать сцену со сложным первичным освещением в несколько миллионов треугольников со скоростью 100 кадров в секунду. 
Чтобы усилить сравнение, будем считать, что мы можем растеризовать сцену с производительностью в 500 кадров в секнду.

\item Тогда за 1 секунду при помощи растеризации можно вычислить порядка 500 точек. При этом, если среднее число лучей на точку поряка 1000,
за 1 секунду при помощи Монте-Карло трассировки лучей можно вычислить освещенность в (50 000 000/1000) = 50 тысяч точках.  
 
\end{enumerate}
                            
\end{enumerate}                                

\noindent Полученный результат не должен вызывать удивления. Алгоритм растеризации для каждой точки вынужден полностью пропускать всю геометию
по графическому конвейеру и плохо подходит для вычисления освещености на множестве точек.
И хотя существуют различные способы, позволяющие в теории значительно сократить время растеризации при использовании такого подхода 
(например многомерная растеризация \cite{Nilsson:2012:DNU:2383795.2383797}), в системе рассчета освещенности на основе трассировке лучей
видится более выгодным использовать трассировку лучей для вычисления освещености в точках кэша.

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

\noindent Фотонные карты на GPU – довольно популярная тема для исследований. 
Одна из основных трудностей при реализации фотонных карт – эффективное построение ускоряющих структур для последующего поиска ближайших фотонов. 
Здесь необходим метод, который бы позволял быстро находить ближайшие фотоны в заданном радиусе вокруг заданной точки трехмерного пространства. 
Реализация такого метода на GPU встречает следующие трудности:

\begin{enumerate}

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

\end{enumerate}


\noindent В работе \cite{Purcell:2003:PMO} алгоритм фотонных карт был впервые реализован на графических процессорах. 
Аппаратные ограничения GPU того времени не позволяли производить запись в произвольную ячейку памяти, 
что создавало целый ряд проблем при реализации ускоряющих структур на GPU.
В \cite{Purcell:2003:PMO} предложено 2 метода построения ускоряющих структур на основе регулярной сетки. 
Первый метод основан на сортировке фотонов по номеру воксела, в которой он попадает. 
После чего при помощи бинарного поиска находился первый фотон в отсортированном массиве фотонов и формировался список индексов.      
Второй метод использует аппаратные возможности на основе буфера трафарета и является болеее быстрым, 
но позволяет сохранять ограниченное число фотонов на каждый воксел. 
Метод реализовывает запись индекса фотона в определенную ячейку используя аппаратный механизм на основе буфера трафарета. 
Недостатком обоих методов из работы \cite{Purcell:2003:PMO} является высокий расход памяти.


\noindent В работе \cite{Zhou:2008:RKC:1409060.1409079} был предложен эффективный алгоритм построения kd-деревьев 
c учетом эвристики VVH \cite{wald:04:VVH} и алгоритм поиска k-ближайших фотонов. 
Алгоритм из работы \cite{Zhou:2008:RKC:1409060.1409079} изначально ориентирован на построения kd-дерева 
над множеством полигонов. Он использует 2 различных подхода в зависимости от количества примитивов в узле и 
хранит данные в динамических списках на GPU. Алгоритм \cite{Zhou:2008:RKC:1409060.1409079} достаточно объемный и сложный в реализации 
и, кроме того, расходует больше памяти, чем его CPU аналог. 
Что касается сбора k-ближайших фотонов - производительность алгоритма сбора низкая вследствие итеративного алгоритма поиска оптимального радиуса сбора.

\noindent В \cite{Frolov_2009}, в целях упрощения алгоритма построения и повышения его эффективности был предложен следующий подход:
Использованое для ускорения трассировки лучей kd-дерево подразбивалось далее так, чтобы листья его
имели сопоставимый с радиусом сбора размер. После трассировки фотонов, для каждого фотона вычислялся индекс узла kd-дерева, в который этот 
фотон попадает и индексы копировались сторону CPU, где при помощи механизма быстрых списков за 2 прохода. 
строилась результирующая ускоряющая структура:

\begin{enumerate}
\item На первом проходе итерировались все индексы листьев для фотонов и для каждого листа 
      kd-дерева составлялся список фотонов попадающих в него. 
\item На втором проходе все списки объединялись в линейный массив, который копировался на сторону GPU.
\end{enumerate}                                       

\noindent Алгоритм построения, таким образом, имеет линейную сложность - $O(N)$. 
Причем, само дерево строилось 1 раз, а при обновлении порции фотонов достаточно было изменить только результирующий массив индексов фотонов. 
Среди недостатков описанного выше подхода следует отметить необходимость априорного знания о геометрии сцены и копирование данных между CPU и GPU.

\noindent В \cite{fabianowski_egsr2009} для ускорения сбора освещености было использовано BVH-дерево и фиксированный радиус сбора.
Среди алгоритмов построения BVH дерева на GPU достаточно эффективны подходы из работ \cite{Garanzha:2011:SFH:2018323.2018333, Karras12}.
Использование BVH дерева в целом является неплохим решением, и основным недостатком этого и других подходов, использующих классические деревья, является рекурсивный алгоритм обхода дерева при сборе освещенности. 
Причины падения производительности при использовании такого обхода - высокая степени дивергентности потоков на GPU и использования локальной памяти.

\noindent В работе \cite{Bogolepov:2011:phmaps} предложен эффективный подход построения ускоряющей структуры 
на основе лексикографической сортировки координат фотонов.
После окончания трассировки все фотоны сортировались в лексикографическом порядке по координатам соответствующих вокселей. 
При таком выборе оператора сравнения все фотоны, принадлежащие одному вокселю, являются эквивалентными и занимают непрерывный сегмент в
отсортированной фотонной карте. Над отсортированным массивом строилась воксельная ускоряющая структура на основе 3D текстур. 
Данный метод похож на первый алгоритм из работы \cite{Purcell:2003:PMO} и его основным недостатком является высокий расход памяти.

\noindent В работе \cite{Alcantara:2009:RPH:1618452.1618500} предложен метод пространственных хэш-таблиц на основе “хэширования кукушки”. 
Авторы не исследовали метод в применении к алгоритму фотонных карт и в работе отмечается, 
что для различных применений могут потребоваться различные хэш-функции.

\noindent Авторы работы \cite{Fleisz_2009} использовали хэш-таблицы для реализации фотонных карт, причем, для того чтобы исключить коллизии во время поиска, 
дополнительно производили сортировку пар (ключ, значение) по ключу. 
При этом поиск ближайших фотонов в точке (x,y,z) производился в ячейке hash(x,y,z) и также в 27 соседних ячейках. 
\cite{Fleisz_2009} отмечает, что для глобальной фотонной карты по сравнению с kd-деревом из работы \cite{Zhou:2008:RKC:1409060.1409079} 
этот метод не дает ускорения сбора с ростом искомого числа ближайших фотонов (т.е. с ростом радиуса сбора). 
Это может быть объяснено потерей производительности на ветвлениях из-за наличия вложенного цикла, что необходимо для поиска в 27 соседних ячейках. 
При большом радиусе сбора и, как следствие, значительном времени итерирования фотонов в листьях, время обхода узлов дерева становится пренебрежимо малым.
Однако, если в лист попадает только часть потоков группы warpы (назовем эту часть активной), 
остальные потоки будут вынуждены ждать до тех пор, пока лист не будет пройден активной частью группы warp.
Но эта проблема одинакого наблюдается как для классических деревьев, так и для хэш-таблиц из работы \cite{Fleisz_2009}.

[НУЖНАКАРТИНКА]

\noindent В работе \cite{Hachisuka:2010:PPP:1899950.1900004} предложен иной подход на основе пространственных хэш-таблиц. 
При возникновении коллизии авторы предлагают стохастически сохранять только 1 фотон на элемент таблицы. 
Данный подход замедляет трассировку фотонов (в 3-4 раза в соответствии со средним числом коллизий), что может быть нецелесообразно, т.к. 
при сложных условиях освещения трассировка фотонов будет занимать много больше времени, 
чем построение ускоряющих структур при помощи любого из известных методов. 
Позже, в работе \cite{Carlberg_2011} на основе метода из \cite{Alcantara:2009:RPH:1618452.1618500} был реализован алгоритм фотонных карт
и было продемонстрировано, что решение на основе метода из работы \cite{Hachisuka:2010:PPP:1899950.1900004} работает медленее чем 
решение на основе метода из работы \cite{Alcantara:2009:RPH:1618452.1618500}.

\noindent Среди практичных подходов следует упомянуть работу \cite{mcguire09imagespace}, авторы которой использовали растеризацию сфер с 
включенным альфа-смешиванием. Это позволяет реализовать сбор освещенности на GPU вообще без построения каких-либо укоряющих структур.
Данный подход аналогичен растеризации сплатов из работы \cite{conf/rt/GautronKBP05}, рассмотренной ранее и имеет тот же набор ограничений - 
а именно, - может быть использован только совместно с растеризацией но неприменим в трассировке лучей. В работе \cite{Mara13Photon} 
в дальнейшем исследовалась производительность различных подходов к решению задачи сбора освещенности в применении к 
интерактивной визуализации на основе растеризации сплатов.

\noindent Также стоит отметить методы, вычисляющие карты освещенности (а также карты светимости) в текстурном пространстве (texture space) \cite{neumannEurographics2005/}.
Данные методы, как правило, вместо сбора освещенности с фотонов используют обратную операцию - распределение освещености с отдельных фотонов
по карте освещенности. Значительное преимущество этого подхода - возможность рассчета освещенности независимо от положения камеры.
Недостаток подхода состоит в необходимости наложения текстурных координат для создания карты освещенности, что не во всех случаях
может быть сделано полностью автоматически с высоким качеством.

Хотя подходы на основе окто-дереьвев используются при реализации фотонных карт достаточно часто, публикаций 
по даной тематике относительно немного. Причина этого заключается в том, что значительное число работ по построению окто-деревьев на GPU
было опубликовано в смежных областях. Одна из работ, использующих окто-деревья на GPU для фотонных карт - \cite{Li:2012:KOT}.
Причем, в работе \cite{Li:2012:KOT} обсуждается эффективная реализация поиска k-ближайших фотонов над окто-деревьями. 
Сравнения, сделанные в данной работе показывают примерно одинаковое качество при использовании сбора освещенности с фиксированным радиусом и
использовании сбора с k-ближайших фотонов. По этой причине далее будем считать, что можно использовать фиксированный радиус сбора и 
рассмотрим некоторые работы, обсуждающие возможность построения окто-деревьев на GPU.

\paragraph{Построение окто-деревьев на GPU.} 
Алгоритмы построения окто-деревьев на GPU используют факт взаимосвязанности окто-дерева и трехмерной Z-кривой, 
которая получается в результате вычисления Z-индекса (или кода Мортона) \cite{morton1966computer}. 
Зафиксируем некоторую глубину дерева ‘k’. Пусть у нас имеется такой набор точек в трехмерном пространстве, что в каждый узел построенного в будущем окто-дерева попадет 
не более чем 1 точка. Если вычислить для точек коды Мортона и затем отсортировать точки по этим кодам, мы получим массив, 
эквивалентный окто-дереву, где в каждом листе хранится ровно по одному элементу (одной точке). 
Выполняя далее поиск в трехмерном пространстве в произвольной точке (x,y,z), мы вычисляем ее код Мортона, 
и применив бинарный поиск для полученного кода Мортона (ZInxed(x,y,z)), находим листовой узел окто-дерева, в который попала точка (x,y,z). 

\noindent В работе \cite{Ajmera:2008:FPG:1342250.1357022} было впервые замечено соответствие между Z-кривой и окто-деревом и использовано свойство отбрасывания младших битов Z-индекса 
для получения идентификаторов старших узлов. Окто-дерево в \cite{Ajmera:2008:FPG:1342250.1357022} представлялось набором массивов L1..Lk, размером 8i элементов для 
каждого уровня глубиной i (т.е. несуществующие узлы хранились в памяти на регулярной сетке). 
Вопрос хранения каких-либо элементов данных в окто-дереве в работе \cite{Ajmera:2008:FPG:1342250.1357022} не рассматривается. 
Таким образом, построенное окто-дерево не может быть использовано непосредственно для пространственного поиска элементов, 
а может только отвечать на вопрос имеются ли искомые элементы в некотором объеме или нет.             

\noindent В работе \cite{Zhou:2011:DOS:1963054.1963112} окто-дерево было использовано для построения поверхности из набора точек при помощи алгоритма марширующих кубов. 
Алгоритм построения использовал серию из параллельных операций уплотнения и сортировки. 
В отличие от работы \cite{Ajmera:2008:FPG:1342250.1357022}, Zhou и др. сохраняют только существенные узлы в памяти, а также сохраняют списки точек на всех уровнях дерева. 
При этом, в \cite{Zhou:2011:DOS:1963054.1963112} было замечено, что имея отсортированный массив точек по их Z-индексу в узле достаточно сохранять лишь 2 смещения, 
обозначающих начало и конец последовательности точек для данного узла в  отсортированном массиве. 
Причем, при переходе от уровня k на уровень k-1 это свойство сохраняется и точки пересортировывать не нужно. 
Достаточно пересчитать указатели соответствующие началу и концу последовательности.

\noindent Если в рассмотренных ранее работах только отдельные уровни дерева строились параллельно, то в работе \cite{Karras12} был предложен метод, 
использующий массивный параллелизм GPU при построении всего дерева. 
Основная идея алгоритма заключается в том, чтобы пронумеровать все узлы дерева при помощи индексов таким образом, 
чтобы индекс родительского узла всегда являлся префиксом по отношению к индексам его дочерних узлов, если рассматривать индексы как 
битовые строки. Тогда, отсортированный массив индексов будет эквивалентен дереву. 
Следует отметить, что свойство префиксов кодов Мортона было описано ранее в работе \cite{Ajmera:2008:FPG:1342250.1357022}.

\subsection{Приближенные методы на GPU}

\noindent Существует достаточно большое число работ, вычисляющих глобальное освещение приближенно.
Сюда можно отнести такие методы как Light Propagation Volumes (LPV) \cite{Kaplanyan:2010:CLP:1730804.1730821}, 
трассировка конусов над воксельными моделями \cite{CNSGE11b}, метод 'Radiance Hits' \cite{conf/egh/Papaioannou11},
методы на основе мгновенной излучательности \cite{Hasan:07:MRCS} и старейший из методов - 
апроксимация вторичного освещения при помощи предрассчитанных сферических гармоник \cite{G03}.

\noindent Данные методы обладают хорошими временными характеристиками. Многие из них широко используются в современных компьютерных играх.
Однако, в данной работе они не рассматриваются детально, поскольку их точность при этих временных характеристиках чрезвычайно низкая.
Если же говорить о тех из них, которые могли бы обеспечить необходимую точность в пределе - \cite{Kaplanyan:2010:CLP:1730804.1730821} и \cite{CNSGE11b},
то при такой точность они работают не быстрее, чем обыкновенная Монте-Карло трассировка путей 
(а в случае метода \cite{Kaplanyan:2010:CLP:1730804.1730821} потребуется еще и невероятно большой объем памяти).                                                               

\section{Программные решения на центральных процессорах}

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

\begin{enumerate}
\item Mental Ray, Pixar RenderMan, Houdini, Arnold, Indigo, Maxwell, VRay, HyperShot, Thea Render, Lightworks'Artisan, ZEANY, Lumion, Clarisse, 
      Twinmotion, DirectViz, keyshot, modo, PhotoView 360, fryrender, SWAP, Turtle, Brazil, final render, Cheetah3D, kerkethea, EIAS9, Team Render (CINEMA 4D).
\item Specter, Inspirer, Inspirer 2, lumicept, Tonatiuh, SolTRACE, rtt , dialux, lighttools, speos, ASAP, tracepro, zemax, DiamCalc.                                                                                                                                                
\item PBRT, Mitsuba, LuxRender, POV-Ray, Sunflow, Yafaray, Pixie, Aqsis, Radiance. 
\end{enumerate}

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

\noindent Такое большое разнообразие существующих программных решений обусловлено в первую очередь чрезвычайно широким спектром применения
рендер-систем и, в некоторых случаях, высокой специализацией (например система DiamCalc ориентирована на реалистичную визуализацию драгоценных камней).
В действительности, даже рендер-системы общего назначения имеют свою специализацию.
Например, системы RenderMan, Houdini, Arnold ориентированы в первую очередь на создание кино-эффектов.
В таких системах важна поддержка огромного спектра эффектов и, что гораздо более важно, возможность тонкой настройки всех 
этапов визуализации. Причем, настройка подразумевает не просто установку предопределенных параметров, а, как правило,
создание специализированного производственного процесса, направленного на визуализацию определенного материала в кино.
В качестве примера можно привести генерацию и визализацию лесного массива в мультфильме САВВА, 
моделирование катастрофы поезда в фильме Метро, визуализация боевых действией в фильме Сталинград (2013 года) и.т.д.          
По этой причине существует специальная высоко-оплачиваемая профессия - технический художник (techical artist).
Эти люди призваны организовывать производственный процесс целиком, превращая работу художников, программистов и актеров в конечный результат.
При этом временные затраты на создание самих моделей и анимации настолько велики, 
что время рассчета самих кадров фильма, хотя и важно, не всегда имеет решающее значения.  

\noindent Другие системы (Например MentalRay, VRay) ориентированны на работу индивидуальных дизайнеров, работающих над архитектурными 
или рекламными проектами в небольших дизайнерских студиях.
В таких системах чрезмерно широкие возможности и настройки не только не нужны, но и создают дополнительные трудности пользователям, 
которым необходимо уметь разбираться иногда не только функциональных возможностях системы, но и в параметрах используемых алгоримов.
Для этих систем важен относительно-простой и быстрый процесс проектирования и визуализации. 
Причем, скорость является критичным параметром, поскольку время визуализации может занимать значительную часть рабочего времени дизайнера. 
Более того, необходимо уметь демонстрировать потенциальному заказчику интерьер и/или экстерьер в различных условиях освещения,
а также различные варианты самого интерьера. По этой причине достаточно востребованной функциональностью рендер-систем является "переосвещение" (relighting) \cite{Relighting:1999} - 
возможность менять условия освещения в реальном времени на основе предрассчитанной информации.
Пользователи таких систем, как правило, профессиональные дизайнеры и художники, имеющие навыки в работе со сложными
пакетами трехмерного моделирования.

\noindent Такие системы как HyperShot еще в большей степени ориентированы на быстрое создание рекламных роликов и презентаций.
Основная цель - создание рекламной презентации за несколька минут, имея на руках сырую модель, спроектированную в САПР пакете.
В таких системах важна не только скорость визуализации, но и скорость создания самой модели. 
Не имея специальных навыков (как в случае предыдущей группы систем), пользователь должен иметь возможность за небольшое число действий достигать конечного результата.

\noindent Наконец, наличие самых разнообразных специализированных программных решений (например программа для создания вируальных моделей мебели - PRO100) 
и существование разнообразных сапр-систем, создает спрос на интегрированные в производственный процесс систем визуализации.

\noindent Вторая причина большого разнообразия рендер-систем - низкая производительность существующих систем и, как следствие, конкуренция систем,
интегрированных в одни и те же пакеты трехмерного моделирования. Здесь следует отметить чрезвычайно популярный пакет 3D Studio Max
компании Autodesk, который используется во многих областях промышленности и бизнеса и служит полем для соперничества
самых различных программных решений.

Что касается систем, ориентированных на высокоточное моделирование, такие системы, как правило, обладают в высокой степени
спецификой решаемой задачи. Среди таких систем необходимо отметить программные и программно-аппаратные комплексы, разработанные в Институте Прикладной Математики имени Келдыша РАН (ИПМ РАН) - Specter, Inspirer, Inspirer 2, lumicept. 
Данные комплексы прочно закрепились в современной промышленности и обеспечивают интеграцию с популярной САПР-системой CATIA. 
Обозначенные выше системы позволяют проектировать сложные оптические системы и решают такие задачи как (\cite{Galaktionov:2006:phdthesis}, \cite{Voloboi:2012:phdthesis}):

\begin{enumerate}

\item Физически-аккуратный расcчет освещенности и построение фотореалистичных изображений сцен содержащих оптически сложные материалы 
при различных условиях искусственного и естественного \cite{Valiev:2008:IMI:1921264.1921274} освещения \cite{Adinetz_physicallyaccurate}.  

\item Моделирование осветительных приборов со сложными свойствами распределения света, сложных оптических систем и устройств, включая системы, содержащие
новейшие микроструктурные рассеивающие материалы, применяемые в современном проектировании оптических светопроводящих систем \cite{Zhdanov:2013:fluorescence}, \cite{Zhdanov:2013:Design}.

\item Измерение оптических свойств плоских образцов материалов \cite{Voloboi:2006:SH}.

\item Моделирование сложных автомобильных красок, красящих покрытий с высокой концентрацией пигмента \cite{Voloboi:2010:KRASKI}, тканей \cite{Lobalzo:2010}.

\item Рассчет карт освещенности и последующая интерактивная визуализация в режиме навигации сцен \cite{Barlyadian:2010:Lightmaps}.

\item Моделирования освещенности и синтез реалистичных изображений через Интернет \cite{Barlyadian:2005:Internet}.

\item Автоматизацию задания параметров сцены с целью упрощения процесса работы пользователей \cite{Voloboi:2012:phdthesis}.

\end{enumerate}

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

Большинство рендер-систем, закрепившихся в промышленности на сегодняшний день используют центральные процессоры.
Эти системы были спроектировны в 80-ых и 90-ых годах прошлого века и в лучшем случае предполагали паралеллизм на основе кластерных систем.     
Причем, паралеллизм не закладывался в алгоритмы изначально, а реализовывался постфактум, в буквальном смысле 'как получится'.
Показательным примером является рендер-система Radiance, использующая алгоритм кэширования освещенности и
работы \cite{Ward:1994:RLS:192161.192286, conf/acpc/KoholkaMG99, conf/egpgv/DebattistaSC06}, 
пытающиеся сделать кэш освещенности параллельным, не меняя оригинального алгоритма и архитектуру системы.

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

\section{Программные решения на графических процессорах}

\noindent Среди существующих рендер-систем реалистичной визуализации,использующих GPU следует отметить тикие проекты как 
Arion, Octane, VRay RT, LuxRender GPU, IRay, OptiX (API), OpenRL (API), Gelato (проект закрыт), RenderAnts, Render Bro (проект по-видимому закрыт), Centi Leo, Cycles Render, RTRays, Hydra.

\noindent Последняя в списке система (Hydra) была создана в результате данной диссертационной работы. 
Она будет детально описана в 4 главе данной работы. 

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

\noindent В результате, несмотря на то, что современные GPU при одинаковой стоимости способоны до 10 раз 
опережать CPU на задаче трассировки лучей \cite{Aila:2009:UER:1572769.1572792}, 
время построения изображения в современных рендер-ситемах на CPU (при примерно одинаковом качестве) 
сравнимо или даже меньше времени построения изображения при помощи систем на основе GPU.   
Причина этого заключается в том, что некоторые методы, рассмотренные в 1 и 2 главах позволяют при сравнимом качестве получаемого изображения 
значительно сократить объем производимых вычислений. За счет этого, CPU решения в ряде случаев продолжают выигрывать у GPU решений.
Известная рендер система Octane, является ярким примером такого случая \cite{OctaneRender}. 
Характерным здесь является тот факт, что Octane позиционируется именно как одна из самых быстрых несмещенных рендер-систем 
(компания-разработчик сравнивает его исключительно c другими 'CPU unbiased renderers'), 
поскольку время расчета изображения в нем может быть в итоге достаточно большим.   


\noindent Среди существующих систем многие используют графичские процессоры лишь частично, для реализации отдельных алгоритмов на GPU 
(например для ускорения пересечение луча и сцены). К таким системам относятся Arion и LuxRender.
Такие решения имеют преимущество в легкости поддержки исходного кода и расширяемости. 
Однако их производительность ограничена во-первых из-за необходимости передачи данных между CPU и GPU и, во-вторых, все же, производительностью CPU кода.

\noindent Среди рендер-систем иcпользующих графические процессоры, следует отметить разработку Нижегородского Государственного Университета имени Лобачевского - 
систему RTRays. Помимо обычной Монте-Карло трассировки путей, данная система реализует собственный оригинальный алгоритм 
рассчета освещенности - усеченную двунаправленной трассировку путей \cite{Bogolepov:2013:phdthesis}. 
Этот алгоритм позволяет рассчитывать значительно более широкий спектр эффектов нежели при помощи обыкновенной Монте-Карло трассировки.
Из других положительны сторон следует отметить расширямость системы, возможность ее использования в обучающем процессе и быстрый алгоритм посроения BVH дерева, реализованный полностью на GPU.

\noindent Стоит отметить и другую отечественную разработку - систему Centi Leo. Как и в системе RTRays, 
в Centi Leo построение BVH дерева реализовано полноcтью на GPU. Однако, Centi Leo - единственная из существующих 
на сегодня систем интегрирования освещенности, способная эфективно работать с большими объемами геометрии и текстур (объемы, не помещающиеся в памяти GPU).

\noindent Следует сказать и о таких системах как OptiX и OpenRL, которые не являются законченными рендер-системами, а представляютиз себя 
специализированные API, предназначенные для создания рендер-систем на их основе. 
Причем, OptiX пренданзачен для разработке систем на основе трассировки лучей для GPU производимых компанией Nvidia, 
а OpenRL имееет в настоящий момент специализированную аппаратную реализацию от компании Caustic Graphics \cite{OpenRL}.

\noindent Хотя специализированный чип Caustic Graphics имеет определенные преимущества, по соотношению цена/производительность в настоящий момент он
проигрывает решению от Nvidia (OptiX). Из недостатков OptiX следует отметить серьезное падение производительности при реализации сложных алгоритмов. 
Это падение связано с тем, что в системе OptiX весь пользовательский код собирается компилятором в 1 большое мега-ядро, 
которое работает по принципу конечного автомата.
Такой подход позволяет достичь значительной гибкости, однако, он негативно влияет на производительность решения 
(в 2-3 раза медленнее чем подход на основе множества небольших ядер) \cite{Laine2013hpg}.
 

