Компьютерная Графика и Мультимедиа Сетевой журнал

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

Содержание

1. Введение
1. Когерентная трассировка лучей
2. ДФОС
3. Формулировка проблемы
2. Набор команд SSE
1. Экскурс в историю
2. Краткий обзор архитектуры SSE
3. Использование SSE-команд
4. Используемый псевдокод
5. Организация данных при работе с SSE
3. Реализация некоторых алгоритмов на SSE
1. Вычисление тригонометрических функций
2. Одновременный двоичный поиск 4-х элементов в массиве
3. Индексирование и интерполяция
4. Заключение
5. Литература

 

Введение

Когерентная трассировка лучей

Когерентная трассировка лучей [1] (ее еще называют SSE-трассировкой лучей) представляет собой подход к трассировке лучей с использованием возможностей большинства современных процессоров выполнять команды, работающие одновременно с несколькими элементами данных. Такие команды обычно называют ОКМД (одна команда - множество данных) - командами. Они позволяют трассировать несколько лучей параллельно одним потоком команд. <Когерентная> означает, что лучи одного набора, трассируемые одновременно, обычно выбираются так, чтобы при работе с ними в большинстве случаев требовались одни и те же (или отстоящие недалеко друг от друга) данные сцены. В таком случае говорят, что лучи <когерентны> (это слово взято в кавычки, так как оно не имеет ничего общего с физическим определением когерентности). Меру сходства между лучами одного пакета иногда называют <когерентностью>.

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

Относительно бурное развитие когерентная трассировка лучей получила с появлением SSE-команд в общедоступных процессорах. Первым таким процессором стал в 1999 году Pentium III Katmai, который обладал набором инструкций SSE (Streaming SIMD Extension, потоковое ОКМД расширение), и содержал команды для обработки наборов из 4-х 32-битных вещественных чисел. В это время AMD выпустила свое собственное расширение, 3D Now!, которое состояло из команд для работы с наборами из 2-х чисел с плавающей точкой. Однако первые работы по когерентной трассировке лучей выполнялись именно на процессорах Pentium и использовали SSE. Впоследствии SSE стало поддерживаться и в процессорах Athlon (все современные процессоры от AMD поддерживают SSE), что сделало SSE фактическим стандартом для набора ОКМД - команд современного процессора. В том числе и для реализации трассировки лучей. Лишним подтверждением этому может служить [2], где описан довольно любопытный трассировщик, который работает интерактивно и использует SSE. Поэтому сама когерентная трассировка стала называться SSE-трассировкой лучей. В дальнейшем в настоящей статье эти понятия считаются синонимами.

ДФОС

ДФОС (Двунаправленная функция отражающей способности) является наиболее общим способом представления отражающих свойств поверхности материала (без учета внутреннего рассеяния). ДФОС для каждой точки поверхности объекта определяет коэффициент переноса энергии между любой парой направлений (направление падения и направление отражения) в этой точке. В общем случае она зависит от свойств материала, длины волны (то есть, цвета) падающего света, его поляризации и т. п. Детальный обзор дан в [3]. Для простоты предположим, что волновые эффекты в сцене отсутствуют (для большинства сцен это довольно реалистичное предположение), и что цвет (луча или ДФОС) можно представить в формате RGB (с вещественными компонентами). Тогда уравнение переноса энергии (действие ДФОС на падающий свет) можно описать формулой:

         (1)

Смысл обозначений этой формулы и более подробные комментарии к ней приведены в [3].

Рисунок 1. Слева: пример экранизации объекта со сложной ДФОС. Справа - пример экранизации того же объекта (та же геометрия, те же источники света) без ДФОС.

ДФОС в общем случае позволяет довольно точно моделировать материалы со сложными оптическими свойствами. К таковым, например, относятся ткани, краски автомобилей и т.п. Если сравнивать с более привычными объектами, то можно сравнить ДФОС с текстурой. Только текстура определяет изменение цвета с изменением положения на поверхности материала, а ДФОС - изменение цвета с изменением направления освещения или наблюдения. Вообще говоря, ДФОС может зависеть и от пространственной координаты, но, как правило, это не так. Возникающие при этом многомерные структуры данных могли бы просто не вместиться в память современных вычислительных машин. Да и в большинстве реальных приложений это не требуется. Поэтому, всюду в дальнейшем мы будем считать, что ДФОС зависит только от направлений входящего и исходящего луча и не зависит от выбора точки на объекте.

Понимая ДФОС именно в смысле функции от двух направлений в одной точке, можно задавать ее различными способами. Простейший и наиболее давно известный способ - задавать ДФОС как простую с математической точки зрения функцию с некоторыми параметрами, которые и определяют свойства материала. Таковы, например, модели ДФОС Фонга или Блинна ([4], [5]). Зачастую говорят просто о <фонговском материале> (например, в OpenGL или в простых трассировщиках), не вводя даже понятия ДФОС, хотя используются те же самые формулы. При помощи простых параметризованных моделей, однако, невозможно задать сложные оптические свойства поверхностей. Иногда, когда требуется добиться реализма, используются более сложные функции. Такие функции иногда называют <ДФОС-шейдерами>, по аналогии с текстурными шейдерами ([6]). Они могут моделировать более сложные свойства, однако не годятся для случаев, когда требуется использовать реальные измеренные ДФОС.

Наиболее универсальным является сеточное (или табличное) представление ДФОС. Для этого вводится некоторая сетка (чаще всего сетка по значениям углов входящего и исходящего луча), в узлах которой задаются значения ДФОС. Между узлами значения ДФОС интерполируются (как показывает опыт, линейной интерполяции вполне достаточно). В отличие, например, от текстур, сетки в ДФОС являются неравномерными. Связано это с тем, что сечение ДФОС для заданного входящего направления имеет профиль с ярко выраженным максимумом. Вблизи этого максимума ДФОС требуется приближать с большой точностью (с мелким шагом), а вдали от него можно использовать достаточно крупный шаг (например, полградуса у максимума и около 30 градусов вдали от него).

Формулировка проблемы

Задача заключается в реализации выборки нужных значений таблично заданной ДФОС c помощью команд SSE для использования в когерентном трассировщике лучей.

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

Набор команд SSE

Экскурс в историю

Первым процессором с набором команд SSE был процессор Pentium III Katmai, вышедший в 1999 году. Создание набора команд SSE было логичным продолжением созданного ранее набора команд MMX (даже новые регистры назывались с точностью до наоборот - XMM). В то время как MMX позволял работать только с целыми числами, SSE добавил возможность работы с наборами из 4-х 32-битными вещественными числами. Примерно в это же время фирма Athlon выпускает свое расширение подобного рода, названное 3DNow! Оно позволяло работать только с наборами из 2-х 32-битных вещественных чисел, и не вводило новых регистров процессора. По ряду причин, SSE получило намного большее распространение, чем 3DNow! Процессоры AMD сейчас также поддерживают SSE. В 2001 году, с выходом процессоров Pentium 4, появился набор команд SSE 2. Новых регистров не появилось, однако появилось более 140 новых команд, которые теперь позволяли на тех же регистрах работать с вещественными числами с двойной точностью, а также с целыми числами (последние по названию не отличались от команд для работы с MMX, только количество элементов увеличилось вдвое). На сегодняшний день SSE 2 поддерживается как и большинством процессоров Intel (Pentium 4), так и новейшими процессорами AMD. Можно считать, что и SSE 2 скоро станет фактическим стандартом, однако в данной статье будет преимущественно использоваться SSE.

Краткий обзор архитектуры SSE

Рассмотрим архитектурную реализацию SSE. А именно - какие для этого вводятся регистры, и какие существуют команды.

По сравнению с предшественниками, в процессорах Pentium III появились 8 новых регистров, которые предназначены для SSE. То есть, в отличие от MMX или 3DNow!, регистры которых физически совпадают с регистрами сопроцессора, SSE-регистры совершенно самостоятельные и не зависят ни от MMX, ни от работы с сопроцессором, что позволяет ускорять реальную работу процессора, если смешивать SSE и MMX-команды, которые при этом будут выполняться почти параллельно. Каждый из этих регистров является 128-битным, и может содержать одновременно набор из 4-х 32-битных вещественных чисел с одинарной точностью (более известные как float).

Имеются команды для загрузки данных регистров из памяти. Причем можно загрузить либо только младшее слово (32 бит) регистра, либо одну из его половинок, либо весь регистр целиком. При загрузке из памяти выполняется традиционное для x86 <переворачивание> по принципу <младший байт в начале>. При загрузке адреса наборов из 4-х вещественных чисел должны быть выровнены по 16-байтной границе в памяти, а наборы из 2-х чисел - по 8-байтовой границе. Нарушение одного из этих правил приводит к Access Violation. Существуют специальные команды, которые позволяют обойти это ограничение, однако они работают медленнее и не рекомендуются к использованию, кроме крайних случаев.

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

addps xmm0, xmm1

Здесь addps - мнемокод команды (эта конкретная команда выполняет параллельное сложение 4-х вещественных чисел), xmm0 - название регистра первого операнда и регистра результата (как и в обычных командах x86, они совпадают), xmm1 - название регистра второго операнда. В качестве второго операнда может выступать операнд в памяти. В качестве первого операнда ячейка памяти может выступать только в командах загрузки операндов, но не в арифметических командах.

К сожалени, объем данной статьи не позволяет рассказать обо всех подробностях архитектуры SSE. Желающие узнать большее могут прочесть книги по программированию на ассемблере [7], [8], [9], где, в частности, рассказывается и про SSE-команды.

Немного о тех командах, которые нам понадобятся. Все SSE-команды можно разбить на следующие группы:

• Команды загрузки данных
• Команды арифметических операций (+, -, *, /, min, max, sqrt, rsqrt, rcp)
• Команды операций сравнения
• Команды операций сравнения
• Команды перестановки данных в регистрах
• Команды логических операций (&, |, ^, andnot)
• Прочие команды (например, получение старших битов из каждого из слов)

Заметим, что, в отличие от сопроцессора, XMM не имеет встроенных команд вычисления некоторых математических функций. Их приходится реализовывать вручную, что и сделано ниже.

Опишем особенности каждого из видов команд. С командами загрузки данных все понятно. Команды арифметических операций, помимо традиционных (выполняющихся одновременно над 4-я вещественными числами), содержат еще команды взятия минимума и максимума и вычисления квадратного корня. Кроме того, существуют две команды rcpps и rsqrtps. Они предназначены для быстрого вычисления аппроксимации обратной величины и обратного квадратного корня. Будучи примерно раз в 8 быстрее, они, однако, дают куда меньшую точность. Возможность применения их в том или ином контексте следует устанавливать экспериментально.

Команды сравнения выполняют все 6 традиционных видов сравнения (<, >, <=, >=, ==, !=) одновременно над 4-я вещественными числами. В отличие от традиционных команд, они не выставляют флаги, а сохраняют результат в регистре первого операнда. Результат состоит либо целиком из 0, либо целиком из 1 (поскольку значение из 32 единиц не является правильным вещественным числом, мы о подобном значении будем говорить как о маске). Подобное решение, с одной стороны, позволяет (через команду movemsk) получить целочисленную маску, а с другой - использовать полученные значения для маскирования в последующих вычислениях, избегая зачастую использования условного оператора. Существует много команд перестановки данных в регистре, однако нам будет нужна только одна из них, да и ей пользоваться непосредственно мы не будем.

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

• Собственно логические операции с масками
• Условная установка значений части из элементов данных. Например, элементы 0 и 2 получаются по одному правилу, а 1 и 3 - по другому. Тогда, при помощи логических операций эти 2 значения можно будет эффективно скомбинировать, даже не выходя за регистры.
• Выполнение некоторых простейших вещественных операций (выделение знака или модуля)

Команда получения старших битов получает старший (знаковый) бит каждого из 4-х слов в регистре общего назначения. У нас она также будет использоваться.

Использование SSE-команд

Существует 3 способа работать с SSE-командами:

• Писать на ассемблере
• Использовать специальные встроенные функции (intrinsics)
• Использование собственных структур для работы с SSE-типами
• Использование классов, которые скрывают использование SSE от пользователя

Поясним эти методы несколько подробнее.

Первый из них - использование ассемблера. При этом можно использовать как чистый ассемблер (например, MASM v7), или же использовать ассемблерные вставки. Поскольку использование ассемблерных вставок облегчает взаимодействие с кодом на языке высокого уровня (например, C++), то лучше использовать ассемблерные вставки, тем более что по эффективности это примерно одно и то же.

При использовании ассемблера требуется использовать стандартные мнемокоды операций и названия регистров (названия регистров могут меняться от компилятора к компилятору). Visual Studio 7 уже поддерживает все мнемокоды SSE, для VC 6 же требуется поставить последний Service Pack или же установить Intel C++ для поддержки компиляторов мнемокодов SSE. Программирование на ассемблере позволяет добиться наибольшей эффективности, однако это наиболее сложный путь.

Второй вариант - использовать специальные встроенные функции, или intrinsics. Условия их использования - те же, что и для ассемблерных мнемокодов. Однако, в отличие от ассемблера, внешне они выглядят как обычные функции на языке высокого уровня, которые принимают в качестве параметра один или несколько аргументов типа __m128, который хранит наборы из 4-х вещественных чисел. Компилятор сам заботится о выравнивании переменных, о размещении их по регистрам и об оптимизации. При этом производительность компилятора Intel C++ в среднем в полтора раза выше, чем у Microsoft C++. Сами по себе встроенные функции почти один в один (за редким исключением) соответствуют командам процессора.

Использование intrinsics примерно на 10 % менее эффективно, чем использование ассемблерного кода. Однако оно и более удобно. В итоге даже пионеры в области SSE-трассировки лучей используют intrinsics, не прибегая к ассемблеру.

Третий вариант - это использование собственных структур и собственных типов данных. На основании Intrinsics можно будет написать свои классы, которые по сути являются одним SSE-значением (содержат единственное поле типа __m128), однако будут содержать еще и операции, облегчающие работу с ними. Можно будет написать несколько классов: один из них будет позволять работать с SSE-значением, как с четверкой вещественных чисел, как с четверкой целых чисел и как с четверкой логических значений (маска). Это придаст программированию более понятный вид, кроме того, введет дополнительный контроль типов, что позволит избежать многих досадных ошибок. Именно в этом ключе все и будут даваться дальшие пояснения. Заметим, что Intel по умолчанию предоставляет некоторые подобные классы, однако по ряду причин они нас не устраивают (в частности, операции сравнения возвращают маску в целочисленном формате, а не в SSE-формате).

Четвертый вариант - использование уже готовых библиотек, которые внутри используют SSE. Можно, например, использовать Intel Performance Primitives (IPP) или уже готовый трассировщик лучей, который поддерживает SSE, и нисколько не заботиться о том, что делается внутри. Хотя это сгодится для многих случаев, тем не менее, для наиболее эффективной реализации именно вашего алгоритма потребуется использование самого SSE.

Используемый псевдокод

Далее для работы с SSE будут использоваться типы intq, floatq и boolq. Подробно они описываются в прилагаемых заголовочных файлах, а здесь лишь вкратце опишу смысл основных операций. Итак, intq - это тип <четверка целочисленных значений>, floatq - <четверка вещественных значений>, а boolq - четверка масок (четверка из 32-битных значений, каждое из которых может быть либо набором нулей, либо набором единиц). Над типами intq и floatq определены операции сложения, вычитания, умножения, максимума, минимума, операции сравнения (они возвращают тип boolq). Кроме того, для вещественного типа определена операция деления, а для целочисленного типа - операция сдвига (при этом второй аргумент может быть как константой, так и intq).

Для типа boolq определены логические операции и, или, исключающее или, а также операция отрицания. В отличие от обычного набора команд, в SSE операция отрицания аппаратно не реализована, так что она определена только для полноты картины. Кроме того, для логических типов определены специальные операции < и >. Поясним их смысл. Мы говорили о том, что в наборе команд есть специальная команда andnot, которая определена следующим образом:

andnot(a, b) = a & (~b).

Подобная операция выполняется над каждым битом. Если понимать логическое значение false как 0, а логическое значение true как 1, то можно заметить, что andnot(a,b) = a > b для всех логических значений a и b. Аналогично andnot(b, a) = a < b. Также можно заметить, что a > b = ~(a <= b), а a <= b - это импликация a -> b. Поэтому в этом смысле определяются еще и операции a <= b и a >= b для логических значений. Однако следует понимать, что они реализованы реально как последовательность команд.

Для типа floatq дополнительно определяются одноместные функции rsqrt и rcp, которые вычисляют соответственно аппроксимацию обратного квадратного корня и обратной величины вещественного числа.

Для всех типов определена операция индексации sse_value[int_i], где sse_value - переменная SSE-типа, а int_i - целочисленное значение (не обязательно постоянное) от нуля до трех. Данная операция возвращает значение простого типа в соответствующей позиции четверки. Для изменения части значений служит операция set(sse_value, boolq_value, new_sse_value), которая заменяет на значения из new_sse_value только те позиции из sse_value, для которых в boolq_value стоит 1. Операция mask(sse_value) соответствует процессорной команде movmskps и позволяет целое число, младшие 4 бита которого - это знаковые биты соответствующих SSE-регистров, а остальные биты равны нулю. Заметим, что между значениями типа boolq и целыми числами от 0 до 15 существует однозначное соответствие, так что можно будет создавать boolq из таких чисел. Кроме того, значения SSE-типов можно будет создавать из одного или 4-х значений обычных типов, при этом последний аргумент функции соответствует младшему разряду регистра. Помимо этого, будет операция shuffle(sse_value, const_int_mask). Она возвращает значение того же типа, что и sse_value образованное по следующему правилу. const_int_mask должно быть константным значением от 0 до 255. Каждая группа из 2-х бит соответствует позиции в SSE-регистре. Так, биты 0 и 1 соответствуют 0-ой (младшей) позиции, а биты 6 и 7 - 3-ей (старшей). Значение из 2-х бит для каждой позиции показывает, из какой позиции в sse_value будет помещено значение в позиции возвращаемого значения. Так, если const_int_mask будет 0, это значит, что значение из 0-го элемента распространяется во все 4 элемента нового значения. С помощью этой функции можно выполнять любые перестановки. Заметим, однако, что реально такой функции не существует, так как _mm_shuffle_ps - intrinsic, который соответствует этой операции, принимает только константные значения, так что оформлять ее в реальном коде придется в виде шаблона и писать что-то типа

shuffle<0>(x).

Однако для ясности не будем перегружать псевдокод подобными сложностями.

По мере надобности будут определяться другие операции, которые потом и будут использоваться.

Организация данных при работе с SSE

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

Наверняка уже написано много кода, которых работает с этими самыми векторами. И переписывать его не хочется. Попытаемся каким-либо образом переписать сами операции, чтобы выполнять их за счет SSE. Заметим, что как правило обрабатываются большие массивы векторов. При этом традиционно используется организация данных, называемая массивом структур. То есть, в памяти идут сначала 3 компоненты первого вектора, потом 3 компоненты второго, и т. д. Вот так, к примеру:

x1, y1, z1; x2, y2, z2; x3, y3, z3; :

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

x1, y1, z1, 0; x2, y2, z2, 0; x3, y3, z3, 0; :

Что ж, потребление памяти выросло на треть. Пусть мы даже этим решили пренебречь. Теперь операция сложения векторов выполняется быстрее, за одну SSE-операцию (хотя КПД только 75 %, ведь один из элементов не задействован). Но попробуем теперь написать скалярное произведение (на ассемблере). Выглядеть оно будет примерно так:

movps xmm0, v0

movps xmm1, v1

mulps xmm1, xmm0

movps xmm2, xmm1

shufps xmm2, xmm2, 0xB1

addps xmm2, xmm1

shufps xmm1, xmm1, 0x02

addss xmm2, xmm1

movss result, xmm2

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

Что мы делали не так? Мы выбрали неправильную структуру данных. Для получения большей эффективности необходимо ввести тип, содержащий в себе одновременно трехмерных вектора. Такая организация называется структурой массивов. Название объясняется тем, что каждое SSE-значение можно рассматривать как маленький массив из 4-х элементов. Что самое главное - подобные структуры очень удобно создавать, используя ранее определенные типы floatq, intq и boolq (маленькое замечание: тип boolq используется исключительно для вычислений или параметров функций, хранить лучше mask(boolq_value), занимающую в 4 раза меньше места). Например, можно достаточно просто написать SSE-класс для трехмерного вектора:

struct vector3dq

{

// : definitions of operators and functions



public:

// vector components

floatq x, y, z;

}; // end of struct vector3dq 

Операции можно определить так же, как и для обычных векторов, только вместо типа float используется тип floatq. Дополнительно потребуется определить операцию типа set (ее семантика определена выше).

Далее можно точно так же определять классы для луча, цвета и т.п. Заметим, что при реализации SSE-трассировки не требуется менять сами структуры хранения данных, так как изменения происходят в алгоритме трассировки - трассируется 4 разных луча, а значения из структур данных просто размножаются в 4 элемента при необходимости. Как следствие - не возникает никаких проблем с выравниванием данных в этих структурах. Однако для достижения максимальной производительности может потребоваться дополнительная реорганизация. Например, может потребоваться <сплющить> многомерный массив в одномерный, чтобы обеспечить эффективное вычисление индексов. Или реорганизовать данные с учетом особенности кэшей процессора, чтобы рядом находились только те элементы данных, которые используются одновременно.

Реализация некоторых алгоритмов на SSE

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

Вычисление тригонометрических функций

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

Однако сначала о том, как можно реализовать такие простые функции, как вычисление модуля и извлечение и изменение знака.

Вспомним, что наш тип floatq содержит внутри значение типа __m128 (назовем его val). Предположим, что объявлена константная переменная (поскольку SSE-команды не предусматривают размещения констант непосредственно в коде, их надо явно хранить в памяти) M128_ABS, которая содержит 4 одинаковых значения, которые имеют битовое представление 0x7fffffff (0 в старшем бите и 1 во всех остальных). Выполняя and этой константы и нашего числа, мы можем получить его модуль. Вот как это выглядит:

floatq abs(floatq x)

{

return floatq(_mm_and_ps(x.val, M128_ABS));

}  // end of abs

Следует сделать несколько замечаний. Во-первых, это касается передачи аргументов. Значения SSE-типов должны всегда передаваться по ссылке. Во-первых, потому что они большие, во-вторых, чтобы не было проблем с выравниванием. То есть, надо было писать не floatq x, а const floatq &x. Но для краткости опустим это. Второе замечание касается функции _mm_and_ps, которая является встроенной функцией и соответствует операции <и>. Если еще принять во внимание вид константы и формать 32-битного вещественного числа, то можно понять, что мы просто ставим 0 в знаковый бит, оставляя все остальные такими, как они есть.

Для извлечения и смены знака потребуется другая константа, M128_SIGN. Она содержит 4 одинаковых значения, имеющих битовое представление 0x80000000 (1 в старшем разряде и 0 во всех остальных). Используя операцию and, можно извлечь знак (при этом они будет представлять т. н. <нуль со знаком>, знаку <+> будет соответствовать битовое представление 0x00000000, а знаку <-> - 0x80000000), при помощи операции xor его можно изменить. Псевдокод здесь приводиться не будет, заметим только, что функции будут иметь вид соответственно floatq sign(floatq x) и floatq xorsign(floatq x, floaq sign).

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

Значения коэффициентов можно узнать в [10]. Вычисление полинома можно организовать по схеме Горнера. В случае выхода x за диапазон (0, 1) считаем, что значение функции не определено. Заметим, что (даже в отладчике) не имеет смысла проверять все элементы на попадание в заданный диапазон, так как в данный момент могут быть <активны> не все элементы. Поясним, что это значит. Например, трассируются 4 луча. В определенный момент только три из них попали в объект, а четвертый - нет. Значит, в маске в четвертой позиции устанавливается false, и действия, которые требуются для этого луча, либо не выполняются, либо выполняются и игнорируются, в зависимости оттого, что эффективнее. В случае с арксинусом эффективнее игнорировать результат (не выполнять что-то можно в том случае, когда есть условные операторы).

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

Ниже приведен псевдокод функции arcsin x:

floatq arcsin(floatq x)

{

floatq sgn = sign(x);

floatq y = abs(x);

floatq res = PI/2 - sqrt(1-y) * (PI/2 + x * (A1 + x * (A2 + x * A3)));

return xorsign(res, sgn);

}  // end of arcsin

Данная функция работает в 18 раз быстрее, чем реализация из стандартной библиотеки языка C). Поясним, что значит <в 18 раз быстрее>. Это значит, что для одинакового количества аргументов SSE-версия требует для вычисления всех значений в 18 раз меньше времени, чем стандартная версия. При этом надо помнить, что SSE-версия вычисляет 4 значения за раз, соответственно, количество вызовов будет в 4 раза меньше. Именно в таком смысле выражение "SSE-версия в x раз быстрее> будет пониматься в дальнейшем.

Немного слов о точности. Точность подхода, приведенного выше, составляет примерно 10^-4. Для нашей конкретной задачи этого вполне достаточно. Если же нужна большая точность, то в той же книге [10] можно найти еще несколько коэффициентов этой формулы.

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

Запишем разложение в ряд Тейлора для sin(x):

Аппроксимация первыми n членами будет иметь вид (член n+1 выделен отдельно):

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

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

Одновременный двоичный поиск 4-х элементов в массиве

Пусть дан массив вещественных элементов, упорядоченный по возрастанию, и дано некоторые вещественное число. Требуется найти такой i для массива arr, что для нашего вещественного числа f выполнено arr[i] <= f < arr[i + 1]. Эта классическая задача решается при помощи стандартного метода двоичного (бинарного) поиска, который давно уже стал классическим алгоритмом, реализован во многих библиотеках и который можно найти в большинстве пособий по компьютерным алгоритмам.

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

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

Для решения подобной проблемы вводится стек текущих масок. В любой момент времени у нас имеется маска, в которой единицы соответствуют активным в настоящий момент элементам искомого значения, а нули - неактивным. Начальная маска передается как параметр функции. Далее, если в какой-то момент времени оказывается, что разные элементы требуют разных сегментов, то одни помещаются в стек, а с другими продолжается поиск. После того, как для них поиск завершен, мы проверяем стек. Если он пуст, то мы закончили поиск. Если стек не пуст, то мы извлекаем маску из стека, и продолжаем с ней поиск. На самом деле в стек помещается не только маска, но и текущий сегмент (его начало и конец). Заметим, что стек для поиска имеет постоянный размер (3, поскольку элементов в SSE-значении 4), что позволяет реализовать его эффективно. Ниже приводится псевдокод процедуры для двоичного поиска SSE-значения:

intq sse_binarysearch(array
 arr, floatq f, bool q qmask) {

i = 0, j = arr.length();

intq res;

int imask = mask(qmask);

while(true) {

if(i == j) {

res.set(i, boolq(imask));

if(stack.empty())

return res;

st.pop(i, j, imask);

continue;

}

k = (i + j) >> 1;

int cmp = mask(val < floatq(arr[k]));

if(all_true(cmp, imask)) 

j = k;

else if(all_false(cmp, imask))

i = k;

else {

push(false_mask(cmp, imask), i, k);

j = k;

imask = true_mask(cmp, imask);

}

}

}

В псевдокоде используются еще не определенные функции all_true, all_false, false_mask и true_mask. Каждая из них имеет 2 параметра: cmp и imask. imask является текущей маской, а cmp - результат сравнения. Функция all_true возвращает true тогда и только тогда, когда для всех позиций, в которых в imask стоит 1, в cmp тоже стоит 1. Функция all_false возвращает true тогда и только тогда, когда для всех позиций, в которых в imask стоит 1, в cmp стоит 0. Функция true_mask возвращает новую маску, в которой 1 стоит только в тех местах, в которых стоит 1 стоит и в imask, и в cmp (то есть - операция <и>). Функция false_mask возвращает новую маску, в которой 1 стоит только в тех местах, в которых стоит 1 в imask и 0 в imask (по сути, это операция andnot, которая была описана ранее; теперь становится понятным, зачем такая операция введена в SSE). Заметим, что в данном примере эти функции определены для целых чисел, и операции производятся над отдельными битами. Их также можно будет определить для типа boolq, при этом в качестве отдельных элементов будут рассматриваться отдельные значения из четверки.

Данный подход позволяет ускорить двоичный поиск примерно в 2,25 раз. Впрочем, здесь следует сказать несколько слов. Реальное ускорение зависит от разброса значений в четверке. Если разброс небольшой, то стек придется использовать только в самом конце, и выигрыш будет больше, чем в случае большого разброса компонентов четверки. В случае сильного разброса можно даже получить потерю в производительности. Но ведь программирование для SSE предполагает когерентность (см. выше), которая как раз и означает небольшой разброс значений. А для отсечения неактивных в текущий момент значений используется маска.

Индексирование и интерполяция

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

Начнем с индексирования. В не-SSE случая оно достаточно тривиально: просто выполняется операция адресной арифметики. В случае SSE возникают некоторые проблемы, которые надо решить. Для начала введем некоторые обозначения. Пусть arr - это вещественный массив, i - переменная типа intq, каждая позиция которой представляет индекс, m - маска в текущий момент, а f - значение, которое необходимо получить.

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

Сначала мы сравниваем все элементы с нулевым элементом (если нулевой элемент в настоящий момент неактивен, то этот шаг просто пропускается). Далее мы выделяем те элементы, которые соответствуют полученной маске. По данной маске производится первая выборка, и этим значением инициализируется число f.

Полученная на первом шаге маска <вычитается> из начальной маски. Для того, что осталось, просто производится поштучная выборка. Псевдокод приведен ниже.

floatq index(array
 arr, intq i, boolq m)

{

floatq res;

boolq m1 = i == intq(i[0]);

int im1 = mask(true_mask(m, m1));

if (im1 != 0)

res = arr[i[0]];

else

res = 0;

int im2 = mask(false_mask(m, m1));

while(m2 != 0)

{

int m3 = least_mask(m2);

int i3 = least_bit(m3);

set(res, arr[i[i3]], boolq(m3));

m2 = false_mask(m2, m3);

}	

}  // end of operator [] 

Функция least_mask() для данной маски находит младший бит, в котором стоит 1, и возвращает маску, в которой 1 стоит только в этом бите. А функция least_bit возвращает номер этого наименьшего бита. Любая такая операция выполняется за один просмотр таблицы. Заметим, что если подобные операции будут выполняться часто, то сама таблица из 16 значений попадет в кэш, что еще более ускорит работу.

Интерполяция используется обычная билинейная (поскольку мы рассматриваем двумерную ДФОС). После того, как все 4 значения выбраны из таблицы, интерполяция заключается всего лишь в применении простейших арифметических операций, и полностью аналогична не-SSE случаю.

Однако все-таки имеются некоторые нюансы, касающиеся реализации интерполяции. Вот они:

1. Многомерный массив лучше хранить в виде одномерного, а индексы вычислять при помощи целочисленного умножения на SSE. Так все будет работать быстрее.
2. Не следует вызывать отдельные функции для индексации массива для получения всех значений, между которыми производится интерполяция. Это приведет к выполнению одной и той же работы по несколько раз. Вместо этого лучше явно переписать код, приведенный выше, и вместо заполнения одной переменной заполнять все сразу (в двумерном случае - 4, в трехмерном - 8).
3. Следует отдельно определять случай, когда все индексы совпадают. В этом случае удается избежать накладных расходов на операции с масками и на хранение переменных веса - можно выполнять операции сразу с существующими данными.

В общем случае тяжело дать универсальные советы относительно индексации и интерполяции. В нашей конкретной задаче подходы, приведенные выше, работали хорошо. В остальных задачах эти подходы могут работать хуже. Например, если вероятность того, что все индексы совпадут, очень мала, то можно не проверять этот случай отдельно. Если достаточно вероятно, что 4 индекса разбиваются на 2 группы по 2 элемента в каждой, то этот случай следует обрабатывать отдельно. Возможно, придется прибегнуть и к низкоуровневому программированию и использовать intrinsics.

Заключение

Полученные результаты приведены в двух таблицах ниже. Результаты даны отдельно для трехмерных и четырехмерных ДФОС.

Для трехмерных ДФОС:

Количество вычисленных значений	100 000	200 000	400 000
не-SSE (сек)	0.085	0.177	0.350
SSE (сек)	0.024	0.050	0.101
Ускорение	3.54	3.54	3.46

Для четырехмерных ДФОС:

Количество вычисленных значений	100 000	200 000	400 000
не-SSE (сек)	0.137	0.248	0.495
SSE (сек)	0.040	0.078	0.156
Ускорение	3.43	3.17	3.17

При этом визуальных отличий между получаемыми изображениями нет.

В данной статье рассказывалось, что такое SSE (вкратце), как с ним можно работать, как с его помощью можно запрограммировать некоторые привычные операции и алгоритмы. Эти алгоритмы потом комбинируются для написания SSE-интерфейса для работы с ДФОС, который позволяет работать с ДФОС SSE-трассировке лучей, и позволяет ускорить вычисление ДФОС более чем в 3 раза. Надо сразу сказать, что можно было бы получить и большее ускорение, если все запрограммировать на ассемблере. Но в общем это дало бы не очень большой выигрыш, да и код стал бы менее читабельным.

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

Литература

1. Ingo Wald, Carsten Benthin, Markus Wagner, Philipp Slusallek. Interactive Rendering with Coherent Ray Tracing. Computer Graphics Forum (Proceedings of EUROGRAPHICS 2001) Vol. 20, № 3, pp. 153-164.
2. Real-Time Ray Tracing Realization. SSE Optimization. http://www.digit-life.com/articles/rtraytracing/index.html. Или просто www.virtualray.ru.
3. Игнатенко Алексей. Использование двухлучевой функции отражательной способности (ДФОС) для моделирования освещения. http://cgm.graphicon.ru/issue3/brdf_intro/index.html.
4. B. Phong, "Illumination for Computer Generated Pictures", Communications of the ACM, vol. 18, no. 6, pp. 311-317, 1975.
5. J. F. Blinn. Models of light reflections for computer synthetized pictures. InComputer Graphics (SIGGRAPH'77 Proceedings), pages 192-198, 1977.
6. Carsten Benthin, Ingo Wald, Philipp Slusallek. A Scalable Approach to Interactive Global Illumination. EUROGRAPHICS, 2003.
7. В. И. Юров. Assembler. Учебник для вузов. Питер, 2003.
8. В. И. Юров. Assembler. Практикум. Питер, 2004.
9. Михаил Гук, Виктор Юров. Процессоры Pentium 4, Athlon и Duron. Питер, 2001.
10. Mathematical Handbook for Scientists and Engineers, Second Edition. McGraw-Hill Book Company, 1968.

Дополнительная информация

Ссылка: 

Андрей Адинец. Применение SSE для вычисления значений ДФОС в когерентном трассировщике лучей. Компьютерная графика и мультимедиа. Выпуск №3(2)/2005. http://cgm.computergraphics.ru/content/view/78

Выпуск: 

Выпуск №3(2)/2005