\documentclass[a4paper,12pt]{report}

\usepackage[utf8]{inputenc}
\usepackage[T2A]{fontenc}
\usepackage[english,russian]{babel}

\usepackage{amssymb,amsfonts,amsmath,mathtext,cite,enumerate,float} %подключаем нужные пакеты расширений
\usepackage{graphicx} 
\graphicspath{{images/}}%путь к рисункам


%\usepackage{courier} 

\usepackage{listings}

\lstset{language=C,
        basicstyle       = \footnotesize\ttfamily,
        keywordstyle     = \bfseries,
        showstringspaces = false,
        captionpos       = b,
        morekeywords     = {include}
}

\usepackage{amsmath}

\begin{document}


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

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

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


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

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


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

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

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

\begin{lstlisting}[label=raytracing1, caption=Трассировка лучей]
float3 RayTrace(Ray ray)
{
  float3 color(0,0,0);
 
  Hit hit = RaySceneIntersection(ray);
  if (!hit.exist)
    return color;
 
  float3 hit_point = ray.pos + ray.dir*hit.t;
 
  for(int i=0;i < lights.size();i++)
  {
   if(Visible(hit_point, lights[i]))
   {
     float R = dist(lights[i].pos, hit_point);
     float I = lights[i].Intensity;
     color += Shade(Ray, I, hit_point, hit.normal)/(R*R);  
   }
  }
   
  if (hit.material.reflection > 0)
  {
    Ray reflRay = reflect(ray, hit); 
    color += hit.material.reflection*RayTrace(reflRay);
  }
 
  if (hit.material.transparency > 0)
  {
    Ray refrRay = refract(ray, hit); 
    color += hit.material.transparency*RayTrace(refrRay);
  }
 
  return color;
} 
\end{lstlisting}


\noindent Поясним фрагмент программы (листинг 1). Луч представлен двумя векторами. Первый вектор – pos – точка испускание луча. Второй – dir – нормализованное направление луча. Цвет – вектор из трех чисел – красный, зеленый, синий. В самом начале функции RayTrace мы считаем пересечение луча со сценой (представленной просто списком объектов пока что) и сохраняем некоторую информацию о пересечении в переменной hit и расстояние до пересечения в переменной hit.t. Далее, если луч промахнулся и пересечения нет, нужно вернуть фоновый цвет (в нашем случае черный). Если пересечение найдено, мы вычисляем точку пересечения hit\_point используя уравнение луча (эквивалентное уравнению прямой с условием t>0). Листинг \ref{raystruct} отображает представление луча в виде структуры. pos - точка начала луча. dir - нормализованное направление луча. Согласно уравнению прямой, точка на прямой вычисляется как p = ray.pos + t*ray.dir. Где t - одномерная координата.

\begin{lstlisting}[label=raystruct, caption=Структура луча на основе уравнения прямой. ]
strucr Ray 
{
  float3 pos;
  float3 dir;
} 
\end{lstlisting}

\noindent Когда мы вычислили точку пересечения в мировых координатах, приступаем к расчету теней и затенения (shading). Пусть источники лежат в массиве lights. Тогда проходим в цикле по всему массиву и для каждого источника света проверяем (той же трассировкой луча), виден ли источник света из данной точки hit\_point. Если виден, прибавляем освещение от данного источника, вычисленное по некоторой локальной модели (например модели Фонга). Этим занимается функция Shade. При расчете затенения важно не забыть поделить интенсивность света источника на квадрат расстояния до него. 

\noindent Если у материала объекта, о который ударился луч, есть отражающие или преломляющие свойства, трассируем лучи рекурсивно, умножаем полученный цвет на соответствующий коэффициент отражения или преломления и прибавляем к результирующему цвету. Коэффициенты reflection и transparency могут быть как монохромными так и цветными. Всё зависит от того, какая используется математическая модель для представления материалов. Для расчета преломленного луча обычно используется закон Снелла. 

\noindent Важным моментом при вычислении отраженного и преломленного лучей является необходимость добавлять небольшое смещение к позиции луча. В случае отраженного луча нужно прибавлять eps*hit.normal, в случае преломленного - (-1)*sign(dot(normal, ray.dir))*eps*hit.normal. Это делается для того чтобы отраженный (или преломленный) луч не ударился на следующем шаге трассировки о ту же самую поверхность. sign(x) = if (x > 0) 1 else -1. dot(x,y) - скалярное произведение векторов x и y.

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

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

\begin{lstlisting}[label=eyeray1, caption=Генерация 'eye' луча для камеры-обскура.]
float3 EyeRayDir(float x, float y, float w, float h)
{
  float fov = 3.141592654f/(2.0f); 
  float3 ray_dir;
  
  ray_dir.x = x+0.5f - (w/2.0f);
  ray_dir.y = y+0.5f - (h/2.0f);
  ray_dir.z = -(w)/tan(fov/2.0f);
	
  return normalize(ray_dir);
}
\end{lstlisting}

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

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


\noindent Где:

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

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


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


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

\begin{lstlisting}[label=eyeray2, caption=Генерация 'eye' луча для камеры с заданной извне матрицей проекции.]
float3 EyeRayDir2(float x, float y, float w, float h)
{
  float4 pos( 2.0f * (x + 0.5f) / w - 1.0f, 
             -2.0f * (y + 0.5f) / h + 1.0f, 0.0f, 1.0f );

  pos = g_mViewProjInv*pos; // inverse(mView*mProj)*pos
  pos /= pos.w;             // homogeneous coordinates

  pos.y *= (-1.0f);

  return normalize(pos.xyz);
}
\end{lstlisting}

\noindent Здесь g\_mViewProjInv - обратная матрица от произведения видовой матрицы (mView) и матрицы проекций (mProj).

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

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

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

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

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

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

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

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

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

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

\begin{lstlisting}[label=raytracing2, caption=Трассировка лучей с корректным учетом преломлений.]
float3 RayTrace(Ray ray)
{
  float3 color(0,0,0);
 
  Hit hit = RaySceneIntersection(ray);
  if (!hit.exist)
    return color;
 
  float3 hit_point = ray.pos + ray.dir*hit.t;
 
  for(int i=0;i < lights.size();i++)
  {
    if(Visible(hit_point, lights[i]))
    {
      float R = dist(lights[i].pos, hit_point);
      float I = lights[i].Intensity;
      color += Shade(Ray, I, hit_point, hit.normal)/(R*R);  
    }
  }
 
  float3 reflection   = hit.material.reflection;
  float3 transparency = hit.material.transparency; 
 
  if(hit.material.type == FRESNEL)
    ApplyFresnel(ray, hit.norm, hit.material.IOR, 
                 &reflection, &transparency); 
 
  if (reflection > 0)
  {
    Ray reflRay = reflect(ray, hit);
    color += reflection*RayTrace(reflRay);
  }
 
  if (transparency > 0)
  {
    Ray refrRay;
    if(refract(ray, hit, &refrRay))
      color += transparency*RayTrace(refrRay);
  }
 
  return color;
} 
\end{lstlisting}

\noindent В приведенном выше листинге считается, что функция 'refract' возвращает false при возникновении полного внутреннего отражения. В соотетствии с \cite{Mitsuba} в листинге \ref{fresnel} мы приводим формулы Френеля для стекла в виде исходных кодов. cosTheta1 - угол между падающим лучом и нормалью. Эта величина может быть вычислена как dot(ray.dir, hit.normal). etaExt - коэфициент преломления внешней среды; etaInt -  внутренней. Обратите внимание, что в листинге \ref{fresnel} ситуация, при которой луч ударяет поверхность снизу (cosTheta1) моделируется при помощи смены местами коэффициентов преломления для etaInt и etaExt. Полученное при помощи функции 'fresnel(...)' значение, таким образом, есть коэффициент отражения; '1.0-fresnel(...)'  - коэффициент пропускания.

\begin{lstlisting}[label=fresnel, caption=Формулы Френеля для диэлектрика \cite{Mitsuba}.]
float fresnelDielectric(float cosTheta1, float cosTheta2, 
                        float etaExt, float etaInt) 
{
  float Rs = (etaExt * cosTheta1 - etaInt * cosTheta2) 
           / (etaExt * cosTheta1 + etaInt * cosTheta2);
  
  float Rp = (etaInt * cosTheta1 - etaExt * cosTheta2)
           / (etaInt * cosTheta1 + etaExt * cosTheta2);

  return (Rs * Rs + Rp * Rp) / 2.0f;
}

float fresnel(float cosTheta1, float etaExt, float etaInt) 
{
  // Swap the indices of refraction if the interaction starts
  // at the inside of the object
  //
  if (cosTheta1 < 0.0f)
    swap(etaInt, etaExt);

  // Using Snell's law, calculate the sine of the angle
  // between the transmitted ray and the surface normal
  //
  float sinTheta2 = etaExt/etaInt * 
                    sqrt(max(0.0f, 1.0f - cosTheta1*cosTheta1));
		
  // Total internal reflection!
  //
  if (sinTheta2 > 1.0f)
    return 1.0f;

  // Use the sin^2+cos^2=1 identity - max() guards against
  // numerical imprecision
  //
  float cosTheta2 = sqrt(max(0.0f, 1.0f - sinTheta2*sinTheta2));

  // Finally compute the reflection coefficient
  //
  return fresnelDielectric(abs(cosTheta1), cosTheta2, etaInt, etaExt);
}
\end{lstlisting}


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

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

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

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

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

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


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

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

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

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


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

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

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


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

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

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


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

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

\noindent Для каждого треугольника матрицу T$_\Delta$ нужно рассчитать один раз и сохранить. Код, вычисляющий пересечение луча и треугольника приведен ниже (здесь m - матрица для преобразования T$_\Delta$).
\begin{lstlisting}[label=rayunittest, caption=Код пересечения луча и треугольника (woops unit test).]
float3 o = mul3x4(m, ray_origin);
float3 d = mul3x3(m, ray_direction);

float t = -o.z/d.z;
float u = o.x + t*d.x;
float v = o.y + t*d.y;
\end{lstlisting}

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

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

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

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

\noindent Нетрудно заметить, что алгоритм работает корректно, рассмотрев различные случаи расположения луча и AABB.

\begin{lstlisting}[label=rayaabb, caption=Код пересечения луча и AABB.]
bool RayBoxIntersection(float3 ray_pos, float3 ray_dir_inv, 
                        float3 boxMin, float3 boxMax, 
                        float& tmin, float& tmax)
{
  float lo = ray_dir_inv.x*(boxMin.x - ray_pos.x);
  float hi = ray_dir_inv.x*(boxMax.x - ray_pos.x);
  
  tmin = fminf(lo, hi);
  tmax = fmaxf(lo, hi);

  float lo1 = ray_dir_inv.y*(boxMin.y - ray_pos.y);
  float hi1 = ray_dir_inv.y*(boxMax.y - ray_pos.y);

  tmin = fmaxf(tmin, fminf(lo1, hi1));
  tmax = fminf(tmax, fmaxf(lo1, hi1));

  float lo2 = ray_dir_inv.z*(boxMin.z - ray_pos.z);
  float hi2 = ray_dir_inv.z*(boxMax.z - ray_pos.z);

  tmin = fmaxf(tmin, fminf(lo2, hi2));
  tmax = fminf(tmax, fmaxf(lo2, hi2));
  
  return (tmin <= tmax);
}

\end{lstlisting}

\noindent В приведенном листинге ray\_dir\_inv - инвертированное направление луча (1.0/ray\_dir), которое для каждого луча можно вычислять 1 раз. Функции fmaxf и fminf вычисляют соттветственно максимальное и минимальное значение из 2 чисел и обычно входят в пакет <math.h> компилятора GCC.


\subsubsection{Функции расстояний (distance-aided ray marching)}

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

\noindent Пусть существует функция f(x,y,z), задающая поверхность уравнением вида f(x,y,z) = 0. Для краткости можно будем писать f(p) = 0. Если в функцию f вместо p подставить, некоторую точку, не лежащую на поверхности, мы получим некое не равное нулю число. Нетрудно заметить, что это число будет представлять собой знаковое (хотя это может зависеть от формулы) расстояние до поверхности.

\begin{figure}[htb]
  \centering
  \includegraphics[width=.5\linewidth]{distance_aided_ray_march}           
  \caption{\label{fig:distance_aided_ray_march}
           Иллюстрация итеративного поиска пересечения с поверхностями, заданным как функции расстояний \cite{DistAidedRayMarch}.}
\end{figure}

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

\noindent Для того чтобы не продолжать марширование по лучу до бесконечности, необходимо в самом начале найти пересечение луча и ограничивающего всю сцену bounding box-а. Как только вы выходите за этот бокс (t > tmax) – остановить реймарчинг. Либо вы просто можете ограничить максимальное число шагов по лучу какой-либо достаточно большой константой.

\begin{lstlisting}[label=raysurf, caption=Поиск пересечения луча и поверхности заданной при помощи функции расстояния.]

Hit DistanceAidedRayMarch(float3 ray_pos, float3 ray_dir, \
                          float tmin, float tmax)
{
  float epsilon = 0.001f;
  
  float minStep = epsilon;
  float maxStep = 10.0f;

  float t = tmin;
  float dt = 1.0f;

  float dist = 1.0f; // signed distance
	
  while(dist > epsilon && t < tmax)
  {
    float3 p = ray_pos + t*ray_dir;
    dist = DistanceEvaluation(p);
    dt = clamp(dist, minStep, maxStep); 
    t += dt;
  }
	
  float3 hitPoint = ray_pos + (t-dt)*ray_dir;

  Hit hit;
  
  if(dist<=epsilon)
  {
    hit.norm = EstimateNormal(hitPoint, epsilon);
    hit.pos  = hitPoint;
  }
	
  return hit;
}
\end{lstlisting}


\noindent В приведенном выше листинге функция DistanceEvaluation как раз и является той f(x,y,z) которая задает поверхность. Функция EstimateNormal вычисляет нормаль на основе того факта что она совпадает с градиентом f(x,y,z):


\begin{lstlisting}[label=raysurfnorm, caption=Вычилени нормали для поверхности заданной при помощи функции расстояния \cite{DistAidedRayMarch}.]
// Calculate the gradient in each dimension 
// from the intersection point
//
float3 EstimateNormal(float3 z, float eps)
{
  float3 z1 = z + float3(eps, 0, 0);
  float3 z2 = z - float3(eps, 0, 0);
  float3 z3 = z + float3(0, eps, 0);
  float3 z4 = z - float3(0, eps, 0);
  float3 z5 = z + float3(0, 0, eps);
  float3 z6 = z - float3(0, 0, eps);
  
  float dx = DistanceEvaluation(z1) - DistanceEvaluation(z2);
  float dy = DistanceEvaluation(z3) - DistanceEvaluation(z4);
  float dz = DistanceEvaluation(z5) - DistanceEvaluation(z6);

  return normalize(float3(dx, dy, dz) / (2.0*eps));
}

\end{lstlisting}


\begin{figure}[htb]
  \centering
  \includegraphics[width=1.0\linewidth]{julia}           
  \caption{\label{fig:julia_fract}
           Фрактал Julia, рассчитанный при помощи трассировки лучей и функций расстояния \cite{DistAidedRayMarch}.}
\end{figure}

\noindent Важное замечание состоит в том, что для некоторых видов поверхностей, часто используемых в компьютерной графике (например NURBS) не только не существует аналитической формулы для вычисления пересечения, но также не существует и аналитического представления в виде f(x,y,z) = 0. Такие поверхности не удасться визуализировать при помощи рассмотренного метода напрямую. 


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

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

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

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

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


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


\noindent Алгоритм построения сетки выглядит следующим образом:

\begin{lstlisting}[label=regulargrid1, caption=Построение регулярной сетки.]

for(int objId=0;objId<objects.size();objId++)
{
  AABB3f box = roundedBoundingBox(objects[objId]);
  for(int x = box.min.x; x <= box.max.x; x++)
  {
    for(int y = box.min.y; y <= box.max.y; y++)
    {
      for(int z = box.min.z; z <= box.max.z; z++)
      {
        if(objects[objId].overlaps(voxel(x,y,z)))  ]
          voxel(x,y,z).appendRef(&objects[objId]);
      }    
    }
  }
  
}

\end{lstlisting}

\noindent В данном листинге для каждого объекта находится его ограничивающий AABB и координаты найденного AABB округляются функцией roundedBoundingBox в нижнюю сторону. Поcле чего для всех вокселей внутри этого AABB проверяется их пересечение с объектами.

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

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

\begin{lstlisting}[label=regulargrid2, caption=Поиск в регулярной сетке.]

float FujimotoTraverse2D(Ray2f ray, AABB2f box)
{
  float t_min;

  if(!IntersectRay2fBox2f(ray, box, t_min))
    return;

  if(t_min < 0)
    t_min = 0;

  // initializing part
  //
  float startX = ray.pos.x + t_min*ray.dir.x;
  float startY = ray.pos.y + t_min*ray.dir.y;

  int x = int( (startX-box.vmin.x)/(box.vmax.x-box.vmin.x)*CELL_NUMBER );
  int y = int( (startY-box.vmin.y)/(box.vmax.y-box.vmin.y)*CELL_NUMBER );

  if(x == CELL_NUMBER) x--;
  if(y == CELL_NUMBER) y--;

  float2 boxSize = box.vmax - box.vmin;

  float tVoxelX, tVoxelY; 
  int stepX, stepY;
  
  if(ray.dir.x > 0)
  {
    tVoxelX  = float(x+1)/float(CELL_NUMBER);
    stepX = 1;
  }
  else
  {
    tVoxelX  = float(x)/float(CELL_NUMBER);
    stepX = -1;
  }

  if(ray.dir.y > 0)
  {
    tVoxelY = float(y+1)/float(CELL_NUMBER);
    stepY = 1;
  }
  else
  {
    tVoxelY = float(y)/float(CELL_NUMBER);
    stepY = -1;
  }

  float voxelMaxX = box.vmin.x + tVoxelX*boxSize.x;
  float voxelMaxY = box.vmin.y + tVoxelY*boxSize.y;

  float tMaxX = t_min + (voxelMaxX - startX)/ray.dir.x;
  float tMaxY = t_min + (voxelMaxY - startY)/ray.dir.y;

  const float voxelSize = (box.vmax.x - box.vmin.x)/CELL_NUMBER;
  const float tDeltaX = voxelSize/fabs(ray.dir.x);
  const float tDeltaY = voxelSize/fabs(ray.dir.y);

  // core of the traversal
  //
  while(x < CELL_NUMBER && x >= 0 && y < CELL_NUMBER && y >= 0)
  {
    float tHit = IntersectAllObjectsInVoxel(x, y, box, ray);

    if(tHit < fminf(tMaxX,tMaxY)) 
      return tHit;   
    
    if(tMaxX < tMaxY)
    {
      tMaxX += tDeltaX;
      x += stepX;
    }
    else
    {
      tMaxY += tDeltaY;
      y += stepY;
    }
  };
  
  return 1e38f; // no hit were found
}

\end{lstlisting}


\noindent Константа CELL\_NUMBER задает разрешение сетки. Переменным VoxelMaxX и VoxelMaxY - соответственно x и y координаты ближайших плоскостей-границ сетки; voxelSize - размер воксела, а tMaxX и tMaxY - расстояние от текущей точки на луче до пересечения с плоскостями VoxelMaxX и VoxelMaxY (рис.~\ref{fig:fujimoto_traverse} ). Модификация части тела цикла перебора вокселей для 3D случая представлена ниже:

\begin{lstlisting}[label=regulargrid3, caption=Модификация основной части поиска в регулярной сетке для трехмерного пространства.]
if(tMaxX <= tMaxY && tMaxX <= tMaxZ)
{
  tMaxX += tDeltaX;
  x += stepX;
}
else if(tMaxY <= tMaxZ && tMaxY <= tMaxX)
{
  tMaxY += tDeltaY;
  y += stepY;
}
else
{
  tMaxZ += tDeltaZ;
  z += stepZ;
}
\end{lstlisting}


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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


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

\noindent Одним из критериев остановки может служить сама функция SAH. Если оцениваемая функцией суммарная стоимость поиска в дочерних узлах больше стоимости поиска в родительском узле, разбиение стоит остановить. Однако, весьма неочевидным фактом является то, что такой критерий остановки далеко не оптимален. Дело в том, что увеличение SAH при разбиении текущего уровня дерева еще не означает замедление трассировки для всего дерева. Рассмотрим цилиндр, составленный из полигонов на рисунке ~\ref{fig:tricilinder}. Узел kd-дерева, обозначенный буквой A является проблеммным. Какую-бы плоскость разбиения мы не выбрали, SAH полученных дочерних узлов будет больше чем SAH родительского узла. Таким образом, SAH в данном случае говорит о том, что подразбиение нужно остановить. Но это не всегда так.
\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{tricilinder}           
  \caption{\label{fig:tricilinder}
           Цилиндр, составленный из треугольников. Какую бы плоскость разбиения мы не выбирали для узла A, это не позволит уменьшить SAH. }   
\end{figure}


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

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

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

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

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

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

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

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

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

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

\begin{lstlisting}[label=kdtree1, caption=Поиск в kd-древе.]
bool KdTreeTraverse(Ray ray)
{
  (hit, t_near, t_far) = RayBoxIntersection(ray, scene_box);
  
  if (!hit || t_near > t_far)
    return false;
  
  if (t_near < 0)
    t_near = 0;
  
  while (true)
  {
  
    while (!node.Leaf())
    {
      axis = node.GetAxis(); // axis = 0..2
      t_split = (node.GetSplitPos() - ray.pos[axis])/ray.dir[axis];
      
      if (t_split < t_near)
      {
        node = node.NearChild();
      }
      else if (t_split > t_far)
      {
        node = node.FarChild();
      }
      else
      {
        stack.push(node.FarChild(), t_far)
        t_far = t_split;
      }
    }

    t_hit = IntersectAllPrimitivesInLeaf(ray,node)
    if (t_hit <= t_far)
	  return true;
    
    if (stack.empty())
	  return false;
    
    t_near = t_far;
    (node, t_far) = stack.pop();
  }
}

\end{lstlisting}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

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

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

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

\begin{lstlisting}[label=bvh1, caption=Поиск разбиения по объектам при построении BVH дерева.]
SplitData FindObjectSplit(PrimitiveList a_plist,AABB3f a_box)
{
  struct BoxLessMax
  {
    BoxLessMax(int a) { axis=a; }
    bool operator()(const PrimitiveRef& p1, const PrimitiveRef& p2) const
    { return (p1.Box().vmax[axis] < p2.Box().vmax[axis]); }
    int axis;
  };
 

  PrimitiveList plists[3] = {a_plist, a_plist, a_plist};

  BVHTree::SplitData res;

  res.sah = SA(a_box)*a_plist.size(); // SA is SurfaceArea
  res.subdivideNext = false;
  res.axis = -1;
  std::vector<AABB3f> rightBounds(a_plist.size());
  
  for (int dim=0;dim<3;dim++)
  {

    // sort data according to axis
    //
    PrimitiveList& plist = *plists[dim];
    std::sort(plist.begin(), plist.end(), BoxLessMax(dim));

    // Sweep right to left and determine bounds.
    //
    AABB3f rightBounds1;
    for (uint i = plist.size() - 1; i > 0; i--)
    {
      rightBounds1.include(plist[i].Box());
      rightBounds[i-1] = rightBounds1;
    }

    // Sweep left to right and select lowest SAH.
    //
    AABB3f leftBounds;

    for (uint i = 1; i < plist.size(); i++)
    {
      leftBounds.include(plist[i-1].Box());
      float CE  = CostEmpty
      float sah = CE+SA(leftBounds)*i+SA(rightBounds[i-1])*(plist.size()-i);
      
      if (sah < res.sah)
      {
        res.sah  = sah;
        res.axis = dim;
        res.primsOnLeftSide  = i;
        res.primsOnRightSide = plist.size() - i;
        res.leftBounds  = leftBounds;
        res.rightBounds = rightBounds[i-1];
      }
   }

  } // for (int dim=0;dim<3;dim++)
  
  res.subdivideNext = (res.axis!=-1) && (res.sah<SA(a_box)*a_plist.size());

  return res;
}
\end{lstlisting}

\noindent Описанный пособ выбора разбиения позволяет найти близкое к оптимуму решение за линейное время. Алгоритм построения самого дерева приведен ниже:

\begin{lstlisting}[label=bvh2, caption=Построение BVH дерева.]

void Subdivide(BVHNode* a_currNode, const AABB3f& a_box, 
               int a_currDeep, PrimitiveList& plist)
{
  if(a_currDeep == 0 || plist.size() <= recomendedPrimInLeaf)
  {
    InsertListInLeaf(a_currNode, plist);
    return;
  }

  SplitData split = FindObjectSplit(plist, a_box);

  if(!split.subdivideNext)
  {
    InsertListInLeaf(a_currNode, plist);
    return;
  }

  // split primitives and construct subtrees
  //
  BVHNode (*leftNode, *rightNode) = NewNodePair(a_currNode);

  iterator middle = plist.begin() + split.primsOnLeftSide;
  
  PrimitiveList leftList(plist.begin(), middle);
  PrimitiveList rightList(middle, plist.end());

  Subdivide(leftNode,  split.leftBounds,  a_currDeep-1, leftList);
  Subdivide(rightNode, split.rightBounds, a_currDeep-1, rightList);
}
\end{lstlisting}


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

\noindent У BVH дерева использующего AABB в качестве ограничивающихпримитивов есть один существенный недостаток, свойственный также и kd-деревьям - треугольники (и некоторые другие объекты) плохо аппроксимируются при помощи прямоугольных паралеллипипедов, оставляя внутри ограничивающего примитива большое количество пустого пространства (рис \ref{fig:bvh_early1}). Даже треугольники, лежащие в плоскостях XOY, XOZ или YOZ обладают площадью поверхности, вдвое меньшей чем ограничивающий их AABB. Для всех остальных треугольников это соотношение еще выше.

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

\noindent На рисунке \ref{fig:bvh_early1} черным цветом обозначены лучи которые не пересекли ни один треугольник (представьте эту картинку в 3D).  Если в геометрии присутствует много длинных и тонких треугольников, получается BVH с огромным количеством пересекающихся друг с другом узлов, содержащих много пустого пространства.   

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{bvh_early2}           
  \caption{\label{fig:bvh_early2}
           Ранее подразбиение. Полученное множество AABB подается на вход построителя BVH дерева.}   
\end{figure}

\noindent Чтобы побороть этот недостаток, можно подразбивать сами треугольники на более мелкие, но существует более элегеантное и экономное в плане памяти решение, называемое ранним подразбиением примитивов (Early Split Clipping) \cite{Ernst:2007:ESC:1524874.1524958}. Идея раннего подразбиения изображена на рисунке \ref{fig:bvh_early2}. Она предлагает представлять треугольник не одним AABB а несколькими. То есть один длинный треугольник подается на вход BVH построителю в виде набора небольших AABB, содержащих ссылки на этот треугольник.

\paragraph{Поск в BVH дереве.} \noindent

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

\begin{lstlisting}[label=bvh3, caption=Узел BVH дерева.]

struct BVHNode
{
  int  leftChildOffset;
  bool isLeaf;
  AABB3f box;
};

\end{lstlisting}

\noindent Рассмотрим представление BVH узла на литинге \ref{bvh3}. Если узел не является листом, поле leftChildOffset следует трактовать как смещение до левого потомка. Смещение до правого потомка при этом вычисляется как leftChildOffset+1. Если узел является листом, leftChildOffset следует трактовать как смещение в массиве списков примитивов, некий указатель передаваемый в функцию IntersectAllPrimitivesInLeaf. Для сокращения записи будем считать присваивание вида '(leftNodeOffset, isLeaf) = node' автоматической распаковкой структуры BVHNode в кортеж из 2 переменных в соответствии с порядком объявления полей в структуре BVHNode. То есть leftNodeOffset = node.leftChildOffset и isLeaf = node.isLeaf. Поле 'box' при этом не используется.

\begin{lstlisting}[label=bvh4, caption=Поиск пересечений в BVH дереве.]
Hit BVHTraversal(Ray ray)
{
  Hit  hit;
  
  int  top=0;
  int  leftNodeOffset = 1;
  bool isLeaf = false;

  bool searchingForLeaf = true;

  while(top >= 0)
  {
    while(searchingForLeaf)
    {       
      float t_minLeft   = 0.0f;
      float t_maxLeft   = 0.0f;
      float t_minRight  = 0.0f;
      float t_maxRight  = 0.0f;
     
      BVHNode leftNode  = GetBVHNode(leftNodeOffset); 
      BVHNode rightNode = GetBVHNode(leftNodeOffset+1);
       
      bool traverseChild0 = RayBoxTest(ray,leftNode,&t_minLeft, &t_maxLeft);
      bool traverseChild1 = RayBoxTest(ray,rightNode,&t_minRight,&t_maxRight);

      traverseChild0 = traverseChild0 && (t_maxLeft >= t_rayMin) \ 
                                      && (t_minLeft <= hit.t);
                                       
      traverseChild1 = traverseChild1 && (t_maxRight >= t_rayMin) \
                                      && (t_minRight <=hit.t);

      // traversal decision
      //
      (leftNodeOffset,isLeaf) = traverseChild0 ? leftNode : rightNode;

      if(traverseChild0 && traverseChild1)
      {
        bool l = (t_minLeft <= t_minRight);
        (leftNodeOffset, isLeaf) = l ? leftNode  : rightNode;
        stack[top]               = l ? rightNode : leftNode;
        top++;
      }

      if(!traverseChild0 && !traverseChild1) // both miss, stack.pop()
      {
        top--;
        (leftNodeOffset, isLeaf) = stack[top];
      }

      searchingForLeaf = !isLeaf && (top >=0);
    }

    if(isLeaf && top >= 0)
    {
      IntersectAllPrimitivesInLeaf(ray, leftNodeOffset, t_rayMin, &hit);
    }
    
    top--;
    (leftNodeOffset, isLeaf) = stack[top];
    
    searchingForLeaf = !isLeaf;
  }
  
  return hit;
}

\end{lstlisting}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

\noindent При этом, если распределение случайных величин не равномерное, оценка для метода Монте-Карло запишется в более общем виде \cite{Sobol:73}:

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

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

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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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

\section{Трассировка путей}

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

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

\begin{lstlisting}[label=pathtracing1, caption=Простейший вид Монте-Карло трассировки путей. В данном псевдокоде emittance обозначает светимость и измеряется в ваттах на квадратный метр.]

Color TracePath(Ray r, int depth) 
{ 
  if (depth == MaxDepth) 
    return Black; // Bounced enough times. 

  Hit hit = RaySceneIntersection(ray);
 
  if (!hit.exist) 
    return Black; // Nothing was hit. 

  Material m = hit.material; 
  if (m.isLight())
    return m.emittance;

  // Pick a random direction from here and keep going. 
  // This is NOT a cosine-weighted distribution! 
  //
  Ray newRay; 
  newRay.origin    = hit.pos; // + hit.norm*epsilon
  newRay.direction = RandomUnitVectorInHemisphereOf(hit.norm); 
  
  // Compute the BRDF for this ray (assuming Lambertian reflection)
  // 
  float cos_theta = dot(newRay.dir, hit.norm); 

  Color BRDF = m.reflectance * cos_theta; 
  return BRDF * TracePath(newRay, depth + 1); 
}

\end{lstlisting}

\noindent Данный алгоритм также называют обратной трассировкой путей. Для генерации случайного вектора по полусфере Вам может быть полезна функция, отображающая случайные числа с равномерным распределением из интервала [0,1] на сферу \cite{Suffern:2007:RTG:1324795}:

\begin{lstlisting} [label=spheresample1, caption=Генерация случайного вектора по сфере из 2 равномерно-распределенных чисел в интервале от 0 до 1.]
float3 MapSampleToSphere(float r1, float r2) // [0..1]
{   
  float phi = r1*3.141592654f*2; //[0..2*PI]
  float h   = r2*2 - 1;          //[-1..1]

  float sin_phi, cos_phi;
  sincosf(phi, &sin_phi, &cos_phi);

  float x = sin_phi*sqrt(1-h*h); 
  float y = cos_phi*sqrt(1-h*h); 
  float z = h;  
  
  return float3(x,y,z);
}
\end{lstlisting}

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

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

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

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

\begin{align*}
C_{0} &= P_{0} \\
C_{1} &= \dfrac{1}{2}C_{0} + \dfrac{1}{2}P_{1} \\ 
C_{2} &= \dfrac{2}{3}C_{1} + \dfrac{1}{3}P_{2} \\
C_{3} &= \dfrac{3}{4}C_{2} + \dfrac{1}{4}P_{3} \\  
&... \\
C_{n} &= \dfrac{n}{n+1}C_{n-1} + \dfrac{1}{n+1}P_{n}
\end{align*}

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

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

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

\begin{align*}
C_{i}     &= \dfrac{ (C_{0} + C_{1} + C_{2} + C_{3} + C_{4} + C_{5}  + ...) }{N}   \\
C_{i}even &= \dfrac{ (C_{0} + C_{2} + C_{4} + C_{6} + C_{8} + C_{10} + ...)}{N/2} \\
C_{i}odd  &= \dfrac{ (C_{1} + C_{3} + C_{5} + C_{7} + C_{9} + C_{11} + ...)}{N/2} \\
Error     &= |C_{i}even -  C_{i}odd|                              \\
C_{i}     &= \dfrac{1}{2}(C_{i}even + C_{i}odd)
\end{align*}

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

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

\begin{lstlisting}
  float cos_theta = dot(newRay.dir, hit.norm);
  BDRF = m.reflectance * cos_theta; 
\end{lstlisting}

\noindent Применим выборку по значимости. Будем генерировать такой вектор отражения, который соотвествует косинусоидальному распределению вокруг направления нормали (рис. \ref{fig:cosine_distribution}).

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

\begin{lstlisting} [label=cosinedist, caption=Генерация вектора с косинусоидальным распределением по нормали \cite{Suffern:2007:RTG:1324795}.]
float3 MapSampleToCosineDistribution(float r1, float r2)
{
  float e = 1.0;
  float sin_phi, cos_phi;

  sincosf(2*r1*3.141592654f, &sin_phi, &cos_phi);
  
  float cos_theta = powf(1.0f-r2,1.0f/(e+1.0f));
  float sin_theta = sqrtf(1.0f-cos_theta*cos_theta);

  float x = sin_theta*cos_phi;
  float y = sin_theta*sin_phi;
  float z = cos_theta;

  return float3(x,y,z);
}
\end{lstlisting}

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

\begin{lstlisting} [label=patracing2, caption=Трассировка путей с учетом косинусоидального распределения отраженных лучей.]
Color TracePath(Ray r, int depth) 
{ 
  if (depth == MaxDepth) 
    return Black; // Bounced enough times. 

  Hit hit = RaySceneIntersection(ray);
 
  if (!hit.exist) 
    return Black; // Nothing was hit. 

  Material m = hit.material; 
  if (m.isLight())
    return m.emittance;

  // Pick a random direction from here and keep going. 
  // This is cosine-weighted distribution 
  //
  Ray newRay; 
  newRay.origin    = hit.pos; // + hit.norm*epsilon
  newRay.direction = RandomCosineVectorOf(hit.norm);

  Color BRDF = m.reflectance; // cos_theta/cos_theta = 1. 
  return BRDF * TracePath(newRay, depth + 1); 
}
\end{lstlisting}

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

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

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

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

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

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

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

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

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

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

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

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

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

\begin{lstlisting}[label=cosinedist2, caption=Генерация вектора со степенным косинусоидальным распределением вокруг вектора direction \cite{Suffern:2007:RTG:1324795}.]
float3 MapSampleToCosineDistribution(float r1, float r2, \ 
                                     float3 direction,   \
                                     float3 hit_norm,    \ 
                                     float power) 
{
  float e = power;
  float sin_phi, cos_phi;

  sincosf(2*r1*3.141592654f, &sin_phi, &cos_phi);
  
  float cos_theta = powf(1.0f-r2,1.0f/(e+1.0f));
  float sin_theta = sqrtf(1.0f-cos_theta*cos_theta);

  float3 deviation(sin_theta*cos_phi, sin_theta*sin_phi, cos_theta);

  // simple orthogonalisation
  //
  float3 ny = direction;
  float3 nx = normalize(cross(ny, float3(1.04f,2.93f,-0.6234f)));
  float3 nz = normalize(cross(nx, ny));
  swap(ny,nz);  // depends on the coordinate system

  float3 res = nx*deviation.x + ny*deviation.y + nz*deviation.z;

  float invSign = dot(direction, hit_norm) > 0.0f ? 1.0f : -1.0f;
  
  // resulting vector is below surface ?
  //
  if(invSign*dot(res, hit_norm) < 0.0f)
    res = -nx*deviation.x + ny*deviation.y - nz*deviation.z;

  return res;
}
\end{lstlisting}

\noindent В приведенном выше листинге использовался некоторый произвольный вектор (1.04f,2.93f,-0.6234f) для реализации ортогонализации Грамма-Шмидта. Более правильный вариант вычисления перпендикуляра преставлен ниже: 

\begin{lstlisting}[label=grammshhmidt, caption=Устойчивая ортогонализация Грамма-Шмидта.]
float3 GetPerpendicular(float3 a_vec) const 
{   
  float3 leastPerpendicular; 
  float  bestProjection; 
  float  xProjection = fabs(a_vec.x); 
  float  yProjection = fabs(a_vec.y); 
  float  zProjection = fabs(a_vec.z); 

  if((xProjection <= yProjection + 1e-5) && \
     (xProjection <= zProjection + 1e-5)) 
  { 
    leastPerpendicular = float3(1, 0, 0); 
    bestProjection = xProjection; 
  } 
  else if((yProjection < xProjection + 1e-5) &&  \
          (yProjection <= zProjection + 1e-5)) 
  { 
    leastPerpendicular = float3(0, 1, 0); 
    bestProjection = yProjection; 
  } 
  else 
  { 
    leastPerpendicular = float3(0, 0, 1); 
    bestProjection = zProjection; 
  } 
  
  return normalize(cross(a_vec, leastPerpendicular));
}
\end{lstlisting}

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

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

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

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

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

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


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

\begin{lstlisting}[label=pathtracing3, caption=Трассировка путей с теневыми лучами.]
Color TracePath(Ray r, int depth) 
{ 
  if (depth == MaxDepth) 
    return Black; // Bounced enough times. 

  Hit hit = RaySceneIntersection(ray);
 
  if (!hit.exist) 
    return Black; // Nothing was hit. 

  Material m = hit.material; 
  if (m.isLight())
    return Black; // different from simple path tracing

  // create shadow ray,
  // eval brdf and compute visibility
  //               
  float3 lpos  = LightSample(light);
  float shadow = Visibility(hit.pos, lpos); 
  
  float  R     = dist(hit.pos, lpos);
  float3 sdir  = normalize(lpos - hit.pos);
  
  // calc implicit probability 
  //
  float cos_theta1 = -dot(sdir, light.norm);
  float cos_theta2 = dot(sdir, hit.norm);
  floart PI = 3.1415926535f;
  float implicitP  = light.area*cos_theta1*cos_theta2/(PI*R*R);
 
  // get brdf value
  Color sbrdf  = m.reflectance;
  
  // summ up all together
  Color shade = shadow*sbrdf*light.emittance*implicitP;

  // Pick a random direction from here and keep going. 
  // This is cosine-weighted distribution 
  //
  Ray newRay; 
  newRay.origin    = hit.pos; // + hit.norm*epsilon
  newRay.direction = RandomCosineVectorOf(hit.norm);

  Color BRDF = m.reflectance; // cos_theta/cos_theta = 1. 
  return shade + BRDF * TracePath(newRay, depth + 1); 
}
\end{lstlisting}


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

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


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


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

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

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

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

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

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

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

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

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

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

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

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

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

\noindent $ \theta $ - угол между нормалью и отраженным лучом. $ p_{e} $ - плотность вероятности для неявной стратегии. Обратите внимание: в данной формуле считается, что сэмплирование ДФО было произведено при помощи косинусодидального распределения (в формуле стоит косинус в первой степени). Если это не так, то $ p_{e} $ нужно вычислять исходя из того, как выбирался случайный отраженный луч. $ p_{i} $ -  плотность вероятности для явной стратегии. $A$ - площадь поверхности источника света. $ w_{e} $ - вес для учета яркости полученной по неявной стратегии. $ w_{i} $ - вес для учета яркости полученной по явной стратегии. Обратите внимания, что $ p_{i} $ зависит от площади источника света, которая, в свою очередь, зависит от глобального масштаба сцены. Результат будет корректен если вклад теневого луча $ L_{i} $ учитывается как было описано ранее с ипользованием вероятности $ implicitP $ . В этом случае, изменение масштаба не повлияет на результат, поскольку площадь сократится в формуле \ref{eq:mis_pt_diffuse}.

\noindent Для того чтобы реализовать многократную выборку по значимости нужно отдельно хранить яркость пришедшую от источника света в направлении теневого луча ($L_{i}$) и в случае если отраженный луч попал в источник света (и дал яркость $ L_{e} $) результирующая яркость для данного к-ого сэмпла (выборки) должна быть вычислена по формуле \ref{eq:mis_pt_diffuse}:

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

\noindent Если в каждой точке пространства построить гистограмму падающего освещения (по углам $ \phi_{i} $ и $ \theta_{i} $) с учетом BRDF и косинуса  угла $ \theta $, то такая гистограмма может быть использована в качестве функции плотности вероятности. Однако для того чтобы не хранить реальную гистограмму в каждой точке, применяется специальный алгоритм оптимизации случайных параметров путя на основе мутаций (идея сродни генетическому алгоритму) после чего все выборки учитываются определенным образом. Для детального изучения алгоритма мы рекоммендуем обратиться к \cite{Cline05apractical}.
 
\begin{lstlisting}[label=MLT, caption=Metropolis Light Transport.]
Image MLT() 
{ 
  Path  X     = InitialPath();
  Image image = ArrayOfZerows();
  
  for(int i=0;i<N;i++)
  {
    Path  Y = mutate(X);
    float a = AcceptProb(Y | X);
    if (rnd(0,1) < a) 
    {
      X = Y;
      RecordSample(image, X);
    }
  }
  
  return image;
}
\end{lstlisting}

\noindent Функция mutate получает следующий путь Y посредством мутации пути X. Функция AcceptProb вычисляется по следующим формулам:

\begin{eqnarray*}
T(Y|X) = p(Y \rightarrow X) \\
a(Y|X) = min\{1,\dfrac{f(y)T(X|Y)}{f(x)T(Y|X)}\}
\end{eqnarray*}

\noindent Где $ T(Y|X) =  p(Y \rightarrow X) $ - вероятность пути X мутировать в путь Y.

\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{mlt1}           
  \caption{\label{fig:mlt1}
           Иллюстрация алгоритма Metropolis Light Transport.}   
\end{figure}

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

\begin{lstlisting}[label=photon0, caption=Простейшее представление фотона. Нормаль сжата в одно 30-ти битное целое число в целях экономии памяти. Сущетствующих способов сжатия нормали достаточно много. Для фотонных карт точности в 30 бит более чем достаточно. ]
struct Photon // 16 bytes per photon
{ 
  float3 pos;
  uint compressedNorm : 30;
  uint colorId : 2; // 4 colors max
} 
\end{lstlisting}

\noindent Авторы \cite{Jensen:2007:HQR:1281500.1281593} используют представление из листинга \ref{photon1} для хранения фотона. В переменных (phi, theta) вместо нормали к поверхности авторы \cite{Jensen:2007:HQR:1281500.1281593} хранят направление, с которого фотон пришел. Это необходимо для вычисления матовых отражений (glossy reflection). Однако, если фотонная карта используется для вычисления диффузной составляющей, вместо этого направления удобно хранить нормаль к поверхности в точке удара (в разделе про сбор освещенности мы покажем для чего).

\begin{lstlisting}[label=photon1, caption=Структура фотона из \cite{Jensen:2007:HQR:1281500.1281593}.]
struct Photon // 20 bytes per photon
{ 
  float3 pos;
  char   color[4]; 
  char   phi, theta;
  short  flag;
} 
\end{lstlisting}

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

\begin{lstlisting}[label=photon2, caption=Структура фотона для CPU.]
struct Photon // 40 bytes per photon
{ 
  float3 pos;
  float3 norm;
  float3 color;
  float  lambda;
} 
\end{lstlisting}

\noindent Для графических процесоров мы рекомендуем использовать паттерн хранения 'Structure Of Arrays' вместо 'Array Of Structure'. Поскольку GPU подерживают 16-битную плавающую точку на аппаратном уровне, в целях снижения нагрузки на шину памяти во время сбора освещенности мы рекомендуем использовать тип half (16-ти битный float).

\begin{lstlisting}[label=photon3, caption=Структура хранилища фотонов для GPU.]
struct PhotonArray // 32 bytes per photon
{ 
  float4 pos[N];   // pos.w is free 
  half4  norm[N];  // norm.w is free
  half4  color[N]; // color.w is free  
} 
\end{lstlisting}


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

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

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

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

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

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

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

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

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

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

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


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

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

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

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

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


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

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

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

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

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

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

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

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

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

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

\end{enumerate}

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

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

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

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

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

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

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

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

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

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

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

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


\noindent Алгоритм марширования по лучу:

\begin{lstlisting}[label=raymarch, caption=Марширование луча внутри объема.]
float3 RayMarching(Ray ray)
{ 
  float tmin, tmax; 
  if(!RayBoxIntersection(ray, boxMin, boxMax, &tmin, &tmax))
    return Black;
  
  float dt = 0.05f;
  float t  = tmin;
	
  float alpha  = 1.0f;
  float3 color = float3(0,0,0);
	
  while(t < tmax && alpha > epsilon)
  {
    float3 p = ray.pos + ray.dir*t;
    float  a = GetOpacity(p);
    color += a*alpha*GatherPhotons(p);
    alpha *= (1.0f-a);
    t += dt;
  }
	
  return color;
}

\end{lstlisting}


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


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

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

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

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

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

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

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

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

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

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

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

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

\begin{lstlisting}[label=icalg, caption=Применение кэша освещенности.]
float3 IrradianceCaching(float3 p, float3 n)
{
  if (InterpolationIsPossibe(p,n))
  {
    return IrradiacneCacheLookUp(p,n);
  }
  else
  {
    float3 I = EvalIntegral(p,n);
    float  R = CalcValidityRadius(p,n);
    InsertRecordInCache(I,R,p,n);
    return I;
  }
}
\end{lstlisting}


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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

\begin{lstlisting}[label=icintrerpolate, caption=Поиск в окто-дереве смножественными ссылками и интерполяция в кэше освещенности.]

float4 IrradianceCacheLookUp(float3 lookPos, float3 lookNorm)
{
  float3 nodePos  = g_octreeCenter;
  float  nodeSize = g_octreeBoxSize;
  OctreeNode node = GetOctreeNode(0); // 0 is root node
  float3 nextNodePos;

  // root to leaf strackless traversal
  //
  while(!node.Leaf())
  {
    float3 diff = lookPos - nodePos;
    int childOffset = GetOctIndex(diff);
    node = GetOctreeNode(node.GetOffsetToChild() + childOffset);

    nextNodePos.x = diff.x > 0.0f ? 1.0f : -1.0f;
    nextNodePos.y = diff.y > 0.0f ? 1.0f : -1.0f;
    nextNodePos.z = diff.z > 0.0f ? 1.0f : -1.0f;

    nodePos += 0.25f*nodeSize*nextNodePos; 
    nodeSize *= 0.5f;
  }

  // interpolation
  //
  float3 E = make_float3(0,0,0);
  float summWt = 0.0f;

  for(int i = 0; i < node.GetNumPoints(); i++)
  {
    float3 pointPos  = node.records(i).pos;
    float3 pointNorm = node.records(i).norm;
    float validityRadius = node.records(i).r;
  
    // perr = ||p1-p2||/(Ri)
    // nerr = (1-sqrt(n1*n2))/(sqrt(1-cos(10))
    //
    float perr = length(lookPos-pointPos)/validityRadius; 
    float nerr = 9.0124*sqrtf(1.f - dot(lookNorm, pointNorm)); 
    float err  = fmaxf(perr, nerr);

    if (err < 1.0f)
    {
      float wt = (1.0f - err);
      summWt += wt;
      E += wt*color;
    }
  }

  E *= (1.0f/summWt);
  float errorEstimate = 1.0f/(summWt + epsilon);

  return float4(E, errorEstimate);
}
\end{lstlisting}

\noindent После обхода всех точек в цикле, полученную освещенность E нужно делить на сумму весов summWt, чтобы получить правильное значение, учтя вклады от каждой точки соответственно их весам. Функция GetOctIndex вычисляет индекс узла на основании 'линейного' представления окто-дерева в памяти (листинг \ref{getoctreeindex}). Функция 'GetOffsetToChild' возвращает смещение в массиве узлов окто-дерева до первого потомка. Следующие 7 потомков лежат сразу за первым в соответствии с листингом \ref{getoctreeindex} и рисунком \ref{fig:octreenodesnumbersee2}.

\begin{lstlisting}[label=getoctreeindex, caption=Вычисления индекса узла окто-дерева.]

int GetOctIndex(float3 pos)
{
  int x1 = (int)(pos.x < 0);
  int y1 = (int)(pos.y < 0);
  int z1 = (int)(pos.z < 0);

  return x1 + (y1 << 1) + (z1 << 2);
}
\end{lstlisting}


\begin{figure}[htb]
  \centering
  \includegraphics[width=0.75\linewidth]{octreenodesnumbers}           
  \caption{\label{fig:octreenodesnumbersee2}
           Нумерация узлов окто-дерева в соответствие с функцией GetOctIndex из листинга \ref{getoctreeindex}. }   
\end{figure}

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

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

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

\begin{enumerate}
\item Кэш освещенности хорошо работает только на сценах с преимущественно гладкими поверхносями. На сценах с большим количеством геометрических деталей (например карты нормалей для имитации микрорельефа) кэш освещености не даст преимущества и может даже замедлить процесс редеринга.
\item Кэш освещенности достаточно сложно сделать параллельным. 
\item Кэш освещенности дает мерцание при анимации. Если стохастическая трассировка путей дает шум, который эффективно удаляется при помощи размытия в движении, то пятна вызванные кэшированием освещенности удалить гораздо труднее. 
\item Кэш освещенности все же снижает точность решения и является смещенным методом. Чем ниже точность, тем больше выигрыш в скорости, но тем более заметны артефакты.
\end{enumerate} 
 
\bibliographystyle{utf8gost705u}  %% стилевой файл для оформления по ГОСТу(украинскому??)
\bibliography{biblio}     %% имя библиографической базы (bib-файла) 



\end{document}