Интерактивная трассировка лучей c использованием SIMD инструкций

Создать новую статью

Дата последнего изменения :   14.08.2009 05:31
Рейтинг
 


Узнайте, как новейшие разработки, такие как Hyper-Threading Technology и двухъядерные процессоры, позволяют оптимизировать расчеты трассировки лучей и добиться трех- и даже пятикратного прироста производительности.

Вступление

С момента появления алгоритмов расчета отражения и преломления света (ray tracing, трассировка лучей) их вычисление обычно рассматривалось как статический процесс и представляло собой весьма трудоемкий способ получения красивых изображений. В последнее время появились возможности расчета лучей в интерактивном режиме, даже при использовании довольно доступного аппаратного обеспечения. С учетом повышения эффективности таких расчетов и благодаря преимуществам современных разработок в области микропроцессоров, в частности, таких как технология Hyper-Threading и многоядерных процессоров, можно ожидать, что в ближайшем будущем эти технологии будут применяться как в интерактивных приложениях, так и в приложениях реального времени. В данной статье рассматриваются возможности низкоуровневой оптимизации вычислений для алгоритмов трассировки лучей света, позволяющие добиться трех- и четырехкратного прироста производительности по сравнению с существующими реализациями.

Как правило, алгоритм трассировки лучей включает следующие этапы:

  • Расчет лучей
  • Трассировка хода лучей
  • Пересечение с примитивами
  • Затенение

В статье приводится обзор более ранних работ по данному вопросу, а также рассматривается реализация указанных элементов в пакетном алгоритме трассировки лучей, оптимизированном под процессорную архитектуру SIMD.

В первой части описывается процесс расчета пакетов лучей. Рассматриваются моменты, когда не следует использовать пакетные операции и как лучше поступать в таких случаях. Во второй части описывается структура kd-tree (k-мерное дерево) и ее применение для быстрого нахождения геометрических примитивов, с которыми должен пересечься луч. В третьей части обсуждаются методы эффективного расчета пересечений пакетов лучей и треугольников. В четвертой части описывается модель затенения, использованная в нашем алгоритме трассировки лучей.

В заключение приводятся результаты замеров производительности и предполагаемые направления дальнейших исследований.

Трассировка лучей

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

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

Основной алгоритм для расчета трассировки лучей выглядит следующим образом:

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

  • Определить ближайшую точку пересечения луча с некой поверхностью
  • Определить цвет поверхности в точке пересечения
  • Скорректировать этот цвет с учетом источников света в сцене.

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

Предыдущие исследования

Для повышения скорости расчетов трассировки лучей по классическому алгоритму Уиттеда (Whitted) предлагались различные решения.

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

Помимо уменьшения числа запросов на построение лучей, скорость работы алгоритма трассировки лучей можно увеличить за счет снижения затрат на отдельные запросы. Это достигается за счет использования методов пространственного деления. Одной из наиболее эффективных структур данных для ускорения расчетов построения лучей является kd-tree5.

В некоторых работах в качестве одного из подходов к решению рассматриваемой проблемы предлагается низкоуровневая оптимизация. По наблюдениям Уальда (Wald6), векторизация вычислений в алгоритме трассировки лучей позволяет приблизительно удвоить производительность типовых программных реализаций6. В исследовании рассматривается пример алгоритма трассировки лучей, который в единицу времени обрабатывает четыре луча, используя набор инструкций SIMD (Single Instruction, Multiple Data) для параллельного выполнения требуемых расчетов. Уальд показывает, что применение для расчетов множества процессоров в сетевой среде позволяет рендерить очень сложные сцены в режиме реального времени; алгоритм трассировки лучей масштабируется практически линейно с увеличением количества процессоров.

Рассматриваемый в статье пример реализации ориентирован на одиночные компьютеры; он заимствует предложенные Уальдом идеи, чтобы продемонстрировать алгоритм трассировки лучей с высокой степенью оптимизации . Данный алгоритм обеспечивает производительность на уровне более 10 млн. лучей в секунду при рендеринге сложных сцен на вполне доступном оборудовании.

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

1. Andrew S. Glassner, An introduction to ray tracing, Academic Press Ltd., London, UK, 1989

2. Bruce Walter, George Drettakis, Steven Parker, Interactive Rendering Using the Render Cache, in Rendering Techniques '99, G. Larson, D. Lichinski (eds), Springer-Verlag, (Proc. 10th Eurographics workshop on Rendering, Granada, Spain.)

3. Bruce Walter, George Drettakis, Donald P. Greenberg, Enhancing and Optimizing the Render Cache, in Rendering Techniques 2002, (13th Eurographics Workshop on Rendering. Pisa, Italy), Springer-Verlag, pp. 37-42, 2002.

4. E. Haines and D. Greenberg, "The Light Buffer: A Shadow-Testing Accelerator," IEEE CG&A, Vol. 6, No. 9, Sept. 1986, pp. 6-16.

5. L.Szirmay-Kalos, V.Havran, B.Balazs, L.Szecsi: "On the Efficiency of Ray-shooting Acceleration Schemes", Proceedings of SCCG'02 conference, Budmerice, Slovakia, Published by ACM SIGGRAPH, PDF, pp. 89--98, April 2002.

6. Ingo Wald, Realtime Ray Tracing and Interactive Global Illumination - PhD thesis, Computer Graphics Group, Saarland University, 2004.

Пакетная трассировка лучей на основе SIMD

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

В случае применения Hyper-Threading Technology на двухъядерной платформе, запросы на построение лучей могут быть равномерно распределены между доступными процессорами. Для облегчения этой задачи можно инициировать два потока трассировки лучей, каждый из которых выполняет половину расчетов всей задачи. Такая схема может без труда применяться на более сложных конфигурациях; производительность операций трассировки лучей масштабируется почти в прямой зависимости от количества доступных процессоров.

Сложнее обстоят дела с использованием технологии SIMD. Алгоритм трассировки лучей становится зависимым от пропускной способности памяти, поскольку каждый луч рассчитывается как проходящий через некоторую пространственную структуру и проверяется на предмет пересечения с несколькими примитивами для определения ближайшей точки пересечения.

Рассмотрим вкратце пример расчетов на основе инструкций SIMD. Технология SIMD позволяет оперировать одновременно четырьмя значениями с плавающей запятой одинарной точности. Эти значения сохраняются в регистрах размером 128 бит. Если, например, мы воспользуемся входящей в набор SSE инструкцией ADDPS (r0, r1), то в результате получим 128-битное значение из четырех слагаемых: первое значение в r0 суммируется с первым значением в r1 и так далее.

Рисунок 1.

Следует заметить, что операция не делает горизонтальный сдвиг; таким образом, производится параллельное выполнение четырех операций, не затрагивающих друг друга.

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

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

Для расчета одного луча потребуются следующие данные:

  • Начальная координата луча
  • Направление луча

В коде эта структура данных следующим образом:

struct ray
{
float ox, oy, oz;
float dx, dy, dz;
};

Преобразуем структуру в векторную форму:

struct RayPacket
{
union
{
struct
{
union { real ox[4]; __m128 ox4; };
union { real oy[4]; __m128 oy4; };
union { real oz[4]; __m128 oz4; };
union { real dx[4]; __m128 dx4; };
union { real dy[4]; __m128 dy4; };
union { real dz[4]; __m128 dz4; };
};
};
};

Поскольку каждое 128-битное значение фактически включает четыре числа с плавающей запятой, мы можем получить доступ к каждому из них путем объединения (union) регистра с массивом чисел с плавающей запятой.

Структура данных алгоритма трассировки лучей теперь описывает четыре луча. На практике, эти лучи представляют собой квадраты 2х2 пикселя, расположенные в непосредственной близости друг к другу. Таким образом, четыре луча в большинстве случае оперируют одними и теми же данными. Кроме того, обычно лучи пересекаются с одним и тем же примитивом и в большинстве случаев проходят через одни и те же вершины дерева kd-tree. Это означает, что векторизация поможет сократить нагрузку на шину памяти. Вместо того, чтобы вызывать вершины структуры kd-tree и примитивы для каждого луча, мы вызываем эти данные для четырех лучей сразу, что теоретически позволяет сократить обмен данными на 75%. В этом случае алгоритм трассировки лучей становится менее зависимым от пропускной способности памяти.

После задания начальных координат и направлений лучей в пакете, мы воспользуемся новой структурой данных для быстрой нормализации лучей:

__m128 v1 = _mm_mul_ps( RP->dx4, RP->dx4 );
__m128 v2 = _mm_mul_ps( RP->dy4, RP->dy4 );
__m128 v3 = _mm_mul_ps( RP->dz4, RP->dz4 );
__m128 sum = _mm_add_ps( v1, _mm_add_ps( v2, v3 ) );
__m128 sqrroot = _mm_sqrt_ps( sum );
__m128 reciprocal = _mm_div_ps( one, sqrroot );
RP->dx4 = _mm_mul_ps( RP->dx4, reciprocal );
RP->dy4 = _mm_mul_ps( RP->dy4, reciprocal );
RP->dz4 = _mm_mul_ps( RP->dz4, reciprocal );

Обратите внимание, что в приведенном примере кода подразумевается существование переменной __m128 под названием 'one', которая содержит четыре значения '1.0f'.

Данный код является SIMD-адаптацией исходной версии кода на языке C:

float sum = RP->dx * Rp->dx + RP->dy * RP->dy + RP->dz * RP->dz;
float reciprocal = 1.0f / sqrt( sum );
RP->dx *= reciprocal;
RP->dy *= reciprocal;
RP->dz *= reciprocal;

SIMD-версия кода уже работает немного быстрее, однако и ее можно улучшить. В набор инструкций SSE-2 процессоров Intel добавлена быстрая функция 1/sqrt(x). Единственным ее недостатком является ограниченная точность. Нам необходима высокая степень точности, поэтому для уточнения результата мы используем метод Ньютона-Рапсона (Newton-Rhapson):

__m128 nr = _mm_rsqrt_ps( x );
__m128 muls = _mm_mul_ps( _mm_mul_ps( x, nr ), nr );
result = _mm_mul_ps( _mm_mul_ps( half, nr ), _mm_sub_ps( three, muls ) );

Данный пример кода предполагает существование переменной __m128 под названием 'half' (четыре раза по 0.5f) и переменной 'three' (четыре раза по 3.0f).

Структура kd-tree и пакеты лучей

Структура kd-tree (k-мерное дерево) представляет собой схему пространственного деления, аналогичную BSP-дереву (двоичное разбиение пространства): пространство разделяется плоскостью пополам, затем каждая половина обрабатывается рекурсивно по той же схеме. Вместо произвольно ориентированных плоскостей разделения, в структуре kd-tree используются выровненные по осям плоскости. С одной стороны, это уменьшает степень свободы при построении дерева, с другой стороны, полученное в результате дерево представляет собой весьма эффективную структуру для обработки запросов на построение лучей. В данной статье акцент будет сделан на прохождение пакетов лучей по kd-tree структурам, нежели чем на их построение.

Базовая операция трассировки с использованием kd-tree структур производится рекурсивно.

Рисунок 2.

На Рисунке 2 показана структура kd-tree с вертикальным разбиением по корневой вершине. Затем каждая половина пространства разбивается по горизонтали. Возможно три варианта прохождение луча через такую структуру:

Рисунок 3.

  • Первый случай иллюстрирует прохождение луча только через «дальний» потомок вершины («дальний» по отношению к направлению луча).
  • Второй случай иллюстрирует прохождение луча через оба потомка вершины.
  • В третьем случае луч проходит только через «ближний» потомок.

Каждая точка на пути луча определяется как p = Oray + t * Dray, где Oray это начальная координата луча, Dray – направление луча и t – расстояние от начала луча. Такое представление луча фактически позволяет хранить точки входа и выхода луча в/из текущей вершины в виде двух значений t: tnear и tfar. Луч входит в вершину в точке tnear и выходит в tfar. Три представленных выше случая теперь могут быть определены путем сравнения расстояния от точки пересечения луча с плоскостью разделения с tnear и tfar. В первом примере точка пересечения находится ближе к началу луча, чем точка вхождения. В третьем примере точка пересечения находится дальше точки выхода луча.

  • В корневой вершине (рисунок 2) луч пересекает плоскость разделения внутри вершины: точка пересечения находится ближе, чем tfar и дальше, чем tnear. Соответственно, мы обрабатываем оба потомка: сначала «ближнюю» вершину (в данном случае, правый потомок) и затем «дальнюю» вершину.
  • В правой вершине-потомке луч снова пересекается с плоскостью разделения внутри вершины. Мы снова обрабатываем оба потомка вершины.
  • Наконец, в левой вершине происходит пересечение только с одним потомком. Следовательно, другой потомок исключается из дальнейшего расчета.

В псевдокоде прохождение одного луча через kd-tree можно представить следующим образом:

while (!node->IsLeaf())
{
d = (node->splitpos – ray.origin[node->axis]) / ray.dir[node->axis]
if (d <= tnear) node = node->far
else if (d >= tfar) node = node->near
else
{
Push( node->far )
node = node->near
}
}

После определения листа, проверяется пересечение луча с примитивами в листе. Более подробное описание данного процесса приводится в следующем разделе.

Если мы попытаемся распространить такой код прохождения луча на пакеты лучей, мы столкнемся с определенной проблемой: лучи в пакете обычно имеют несколько различные направления и поэтому необязательно проходят строго через те же самые вершины дерева k-пространства.

Данная проблема решается путем расчета прохождения через вершину всех лучей пакета, если хотя бы один из лучей фактически через нее проходит. Результаты для “лишних” лучей исключаются по маске. “Лишние” лучи определяются путем сравнения tnear с tfar (лучи, которые не должны проходить через определенную вершину, не отвечают условию tfar >= tnea)r.

Ниже приводится пример псевдокода для расчета прохождения четырех лучей через kd-tree с использованием векторных инструкций.

mask = all_valid;
while (!node->IsLeaf())
{
splitpos[0…3] = node->splitpos )
dist[0…3] = (splitpos[0...3] – rp.origin[0...3][axis]) / (rp.dir[0...3][axis])
comparisonresult = (dist[0…3] <= tnear[0...3])
comparisonresult &= mask
if (all of comparisonresult[0...3] = 'true')
{
node = node->far
continue
}
comparisonresult = (dist[0...3] >= tfar[0…3])
comparisonresult &= mask
if (all of comparisonresult[0...3] = 'true')
{
node = node->near
continue
}
if (any of comparisonresult[0...3] > 0)
{
Push( node->far, dist[0...3], tfar[0...3] )
node = node->near
mask = (tnear[0...3] < tfar[0...3])
tfar[0...3] = dist[0...3] where mask = 'true'
}
}

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

Ниже показана оптимизированная реализация приводимого символического кода на языке C.

const int offs[4] = { (RP->dcell[0] >= 0)?1:0,  (RP->dcell[4] >= 0)?1:0,
(RP->dcell[8] >= 0)?1:0, 0 };

while (!node->IsLeaf())
{
const __m128 spos = _mm_load_ps( node->m_Split );
const int aidx = node->GetAxis();
const __m128 d4 = _mm_mul_ps( _mm_sub_ps( spos, RP->oc4[aidx] ), RP->rdc4[aidx] );
const KdTreeNode* restrict ln = node->GetLeft() + offs[aidx];
if (!_mm_movemask_ps( _mm_and_ps( _mm_cmpgt_ps( d4, tn4 ), mask ) )) { node = ln; continue; }
node = node->GetLeft() + (offs[aidx]^1);
if (_mm_movemask_ps( _mm_and_ps( _mm_cmplt_ps( d4, tf4 ), mask ) ))
{
const __m128 mask2 = _mm_cmpgt_ps( d4, tn4 );
const __m128 mask3 = _mm_cmplt_ps( d4, tf4 );
m_Stack[stackptr].tf4 = tf4;
tf4 = _mm_or_ps( _mm_and_ps( mask3, d4 ), _mm_andnot_ps( mask3, tf4 ) );
m_Stack[stackptr].node = (KdTreeNode*)ln;
m_Stack[stackptr].mask = mask;
m_Stack[stackptr++].tn4 = _mm_or_ps( _mm_and_ps( mask2, d4 ), _mm_andnot_ps( mask2, tn4 ) );
mask = _mm_cmplt_ps( tn4, tf4 );
}
}

В рассматриваемом выше примере статический массив 'offs' используется для определения «ближних» и «дальних» вершин с учетом направления лучей. Наш алгоритм трассировки лучей сохраняет значения для левого и правого потомков в последовательных ячейках памяти; «правый» потомок, таким образом, имеет смещение 1 от «левого» потомка. Следует заметить, что точка пересечения рассчитывается на основании заранее полученных обратных величин направления луча.

Поскольку деление потомков на «ближние» и «дальние» зависит от знака направления луча, все лучи в пакете лучей должны иметь одинаковые знаки. Обычно это и имеет место; в противном случае мы трассируем лучи независимо, с использованием альтернативного кода для «одиночного луча».

Напоминаем, что пакет лучей определяется следующим образом:

struct
{
union { real ox[4]; __m128 ox4; };
union { real oy[4]; __m128 oy4; };
union { real oz[4]; __m128 oz4; };
union { real dx[4]; __m128 dx4; };
union { real dy[4]; __m128 dy4; };
union { real dz[4]; __m128 dz4; };
}

Алгоритм для одиночных лучей применяется в пакете в тех случаях, когда знаки в точках dx4, dy4 или dz4 различаются. Это эффективно определяется с помощью встроенной функции _mm_movemask, которая обычно используется для конвертирования результата сравнения векторов в 4-битное целое число. Поскольку такая конвертация проверяет только знаковый бит каждого из четырех значений с плавающей запятой, ее использование для проверки знака весьма эффективно.

int m1 = _mm_movemask_ps( rp->dx4 );
int m2 = _mm_movemask_ps( rp->dy4 );
int m3 = _mm_movemask_ps( rp->dz4 );
if (((m1 == 0) || (m1 == 15)) && ((m2 == 0) || (m2 == 15)) && ((m3 == 0) || (m3 == 15)))
valid = true;

Пересечение пакетов с треугольниками

Последней составляющей запроса на построение луча является пересечение с примитивами. Данные о геометрических примитивах хранятся в листьях kd-tree. При нахождении листа кодом трассировки, пакет лучей пересекается с хранящимися в таком листе геометрическими примитивами.

Быстрым методом определения пересечений лучей с треугольниками является использование барицентрической (barycentric) системы координат:

for ( each primitive )
{
prim = node->prim[…]
k = prim->k
d = 1 / ray.D[k] + prim->nu * ray.D[ku] + prim->nv * ray.D[kv]
t = (prim->nd – ray.O[k] - prim->nu * ray.O[ku] - prim->nv * ray.O[kv]) * d
if (!(dist > t && t > 0)) continue
hu = ray.O[ku] + t * ray.D[ku] - prim->au
hv = ray.O[kv] + t * ray.D[kv] - prim->av
u = hv * prim->bnu + hu * prim->bnv
if (u < 0) continue
v = hu * prim->cnu + hv * prim->cnv
if ((v < 0) || ((u + v) > 1)) continue
// intersection found
}

Данный пример кода предполагает наличие следующих заранее полученных значений для каждого треугольника (векторы A, B и C описывают треугольник):

vector3 c = B - A;
vector3 b = C - A;
vector3 N = b.Cross( c );
// determine dominant axis
int k = 2, u, v;
if (_fabs( N.x ) > _fabs( N.y)) { if (_fabs( N.x ) > _fabs( N.z )) k = 0; }
else { if (_fabs( N.y ) > _fabs( N.z )) k = 1; }
u = (k + 1) % 3, v = (k + 2) % 3;
nu = N.cell[u] / N.cell[rayidk];
nv = N.cell[v] / N.cell[rayidk];
nd = N.Dot( A ) / N.cell[rayidk];
bnu = b[u] / (b[u] * c[v] – b[v] * c[u]);
bnv = -b[v] / (b[u] * c[v] – b[v] * c[u]);
cnu = c[v] / (b[u] * c[v] – b[v] * c[u]);
cnv = -c[u] / (b[u] * c[v] – b[v] * c[u]);
au = A[u], av = A[v];

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

В приводимом ниже примере кода показан способ реализации распараллеленной версии рассматриваемого выше псевдокода:

const static unsigned int modulo[] = { 0, 1, 2, 0, 1 };
const int k = prim->k;
const unsigned int _ku = modulo[k + 1] << 2;
const unsigned int _kv = modulo[k + 2] << 2;
const unsigned int _ka = k << 2;
const __m128 &da = *(__m128*)&RP->dcell[_ka];
const __m128 &du = *(__m128*)&RP->dcell[_ku];
const __m128 &dv = *(__m128*)&RP->dcell[_kv];
const __m128 &oa = *(__m128*)&RP->ocell[_ka];
const __m128 &ou = *(__m128*)&RP->ocell[_ku];
const __m128 &ov = *(__m128*)&RP->ocell[_kv];
const __m128 inu = _mm_load_ps ( &ta->nu );
const __m128 inv = _mm_load_ps ( &ta->nv );
const __m128 ind = _mm_load_ps ( &ta->nd );
const __m128 v2 = _mm_add_ps( _mm_add_ps( _mm_mul_ps( inu, du ), _mm_mul_ps( inv,dv ) ), da );
const __m128 v1 = _mm_sub_ps( _mm_sub_ps( _mm_sub_ps( ind, oa ), _mm_mul_ps( ou, inu ) ),
_mm_mul_ps( ov, inv ) );
const __m128 t4 = fastxovery( v1, v2 );
const __m128 comp_a = _mm_cmpgt_ps( ID->dist4, t4 );
const __m128 comp_b = _mm_cmpgt_ps( t4, zero );
const __m128 mask = _mm_and_ps(comp_a, comp_b);
if (_mm_movemask_ps( mask ) > 0)
{
const __m128 au4 = _mm_load_ps( ta->au );
const __m128 av4 = _mm_load_ps( ta->av );
const __m128 v1 = _mm_sub_ps( _mm_add_ps( ou, _mm_mul_ps( t4, du ) ), au4 );
const __m128 v2 = _mm_sub_ps( _mm_add_ps( ov, _mm_mul_ps( t4, dv ) ), av4 );
const __m128 ibnu = _mm_load_ps( ta->bnu );
const __m128 ibnv = _mm_load_ps( ta->bnv );
const __m128 icnu = _mm_load_ps( ta->cnu );
const __m128 icnv = _mm_load_ps( ta->cnv );
const __m128 vu4 = _mm_add_ps( _mm_mul_ps( v2, ibnu ), _mm_mul_ps( v1, ibnv ) );
const __m128 mask2 = _mm_and_ps( mask, _mm_cmpge_ps( vu4, zero ) );
const __m128 vv4 = _mm_add_ps( _mm_mul_ps( v1, icnu ), _mm_mul_ps( v2, icnv ) );
const __m128 mask3 = _mm_and_ps( mask2, _mm_cmpge_ps( vv4, zero ) );
union { __m128 mask4; __m128i imask; };
mask4 = _mm_and_ps( mask3, _mm_cmple_ps( _mm_add_ps( vu4, vv4 ), one ) );
const __m128i tacc4 = _mm_set1_epi32( (unsigned int)ta );
ID->tacc4 = _mm_or_si128( _mm_andnot_si128( imask, ID->tacc4 ), _mm_and_si128( imask, tacc4 ) );
ID->dist4 = _mm_or_ps( _mm_andnot_ps( mask4, ID->dist4 ), _mm_and_ps( mask4, t4 ) );
ID->u4 = _mm_or_ps( _mm_andnot_ps( mask4, ID->u4 ), _mm_and_ps( mask4, vu4 ) );
ID->v4 = _mm_or_ps( _mm_andnot_ps( mask4, ID->v4 ), _mm_and_ps( mask4, vv4 ) );
}

Теперь мы имеем все составляющие для параллельной трассировки четырех лучей в структурах kd-tree, при условии соответствия лучей определенным правилам. Лучи должны располагаться как можно ближе друг к другу и знаки их направлений должны совпадать. В этом случае становится возможным применение векторных инструкций для ускорения обработки типичных запросов на построение лучей, ценой некоторого усложнения кода. Для пакетов, не соответствующих указанным критериям, необходимо применять неоптимизированный альтернативный код построения одиночных лучей.

Производительность

Для определения производительности мы протестировали пакетный алгоритм трассировки лучей в нескольких сценах. Производительность пакетного алгоритма сравнивалась с алгоритмом расчета отдельных лучей в равные условиях.

Первая тестовая сцена – «игрушка» – довольно проста и состоит из 10 тыс. треугольников. Сцена включает пустое экранное пространство, что создает идеальные условия для пакетного алгоритма, так это обычно подразумевает сравнительно большие вершины kd-tree, что снижает вероятность отклонения лучей в пакете.

Вторая сцена – «галерея» – состоит из 81 тысячи треугольников. Данная сцена сравнительно сложна для пакетного алгоритма, так как характеризуется большими затратами на расчет затенений и имеет равномерно распределенную геометрию.

Третья сцена – «кухня» – состоит из 181 тысячи треугольников и имитирует реальную среду. Она содержит несколько сложных участков и несколько участков с низкой детализацией.

Производительность в сценах измерялась на ноутбуке с процессором Mobile Intel® Pentium® 4 с тактовой частотой 1700 МГц. Программа трассировки лучей выполняла рендеринг изображений с разрешением 800 на 600 пикселей.

Сцена

Пакетные расчеты

Одиночные расчеты

Игрушка

3277k лучей/с, 4.6 fps

1064k лучей/с, 1.9 fps

Галерея

1554k лучей/с, 1.4 fps

283k лучей/с, 0.2 fps

Кухня

2135k лучей/с, 2.5 fps

452k лучей/с, 0.6 fps


Как следует из приведенных выше цифр, различия между пакетным алгоритмом и алгоритмом расчета одиночных лучей весьма значительны, особенно в сцене «галерея». Крайне большой разброс результатов для данной сцены объясняется тем, что последовательный пакет также задействует код SSE для расчета затенения. Сцена «галерея» отличается мощной текстуризацией с применением билинейного фильтра. Для последовательных пакетов фильтрация выполняется по четырем пикселям одновременно.

При оценке прироста производительности запросов на построение лучей без учета затенения типичный прирост скорости составляет 300%.

Мы показали, что производительность алгоритма трассировки лучей может быть значительно увеличена за счет векторизации вычислений для параллельной трассировки четырех лучей. Распараллеливание привносит в алгоритм некоторые сложности, однако, способствует значительному улучшению производительности.