CC BY
110
Информационные технологии Вестник Нижегородского университета им. Н.И. .Лобачевского, 2011, № 3 (2), с. 180-186
УДК 004.925.3
МОДЕЛИРОВАНИЕ КАУСТИК В РЕАЛЬНОМ ВРЕМЕНИ НА ОСНОВЕ КОМБИНИРОВАННЫХ ВОЗМОЖНОСТЕЙ OpenCL И ШЕЙДЕРОВ
© 2011 г. Д.К. Боголепов, В.Е. Турлапов
Нижегородский госуниверситет им. Н.И. Лобачевского
denisbogol@gmail.com
Поступила в редакцию 27.01.2011
Рассмотрен метод фотонных карт для моделирования каустик в реальном времени. Предлагается быстрый алгоритм, адаптированный для ГПУ и реализованный на основе OpenGL/OpenCL. Для высокопроизводительной визуализации и построения фотонной карты используются шейдеры OpenGL. Эффективное использование фотонной карты обеспечивает воксельная ускоряющая структура, которая конструируется на ГПУ средствами OpenCL. Приводятся оценки производительности.
Ключевые слова: интерактивная визуализация, трассировка лучей, метод фотонных карт, моделирование каустик, GPGPU, OpenGL, GLSL, OpenCL.
Введение
Эффекты глобального освещения играют важную роль в задачах синтеза реалистичных изображений. К базовым эффектам можно отнести тени, многократные отражения и преломления. Для их моделирования в большинстве случаев применяется алгоритм трассировки лучей [1]. Расширенные эффекты, связанные с диффузным отражением1 света и каустиками2, требуют более сложных алгоритмов. К числу наиболее распространенных относится метод фотонных карт, который обеспечивает практическое и вычислительно эффективное решение для задач глобального освещения [2]. Алгоритмы глобального освещения и метод фотонных карт, в частности, традиционно полагаются на сложные программные реализации и статическую визуализацию. Тем не менее, появление производительной программируемой графической аппаратуры позволило сделать большой
1 Диффузное отражение света - его рассеивание неровной поверхностью среды по всем возможным направлениям. Пространственное распределение отражённого потока излучения и его интенсивность различны в разных конкретных случаях и определяются соотношением между длиной волны X и размерами неровностей, распределением неровностей по поверхности, условиями освещения, свойствами отражающей среды.
2 Каустика - геометрическое место всех фокусов
негомоцентрических пучков света (не сходящихся в одной точке). На практике каустики легко заметить: в них резко возрастает интенсивность светового поля.
шаг в направлении интерактивного применения данных методов.
Одна из первых успешных реализаций метода фотонных карт на графическом процессоре описана в работе [3]. Авторы отказались от использования традиционного kd-дерева для ускорения доступа к фотонной карте. Вместо этого, в работе используется регулярная сетка, и предлагаются два способа распределения фотонов по вокселям. Кроме того, авторами предложен метод поиска к ближайших фотонов сетки для заданной точки соударения. Данный метод оказался весьма ресурсоемким и занимал порядка 90% от времени визуализации. Безусловно, используемая на тот момент аппаратура класса NVIDIA GeForce FX 5900 Ultra не могла обеспечить интерактивной визуализации даже для самых простых сцен. Например, сцена с одним прозрачным объектом в окне 512 х 512 точек для 65 х 103 фотонов обрабатывалась более одной минуты. Хотя производительность уступала существующим реализациям для центрального процессора, работа была крайне интересной и послужила основой для дальнейших улучшений. К их числу относятся работы [4, 5], в которых авторы исследовали различные ускоряющие структуры, включая kd-деревья и иерархию ограничивающих объемов (BVH -bounding volume hierarchy). Однако применение современного оборудования не позволило достичь реального времени даже на простых тестовых сценах. Например, в работе [5] на простых сценах в окне 512 х 512 точек для
128 х 103 фотонов получена производительность порядка 4-6 кадров в секунду на ГПУ NVIDIA GeForce GTX 280.
Целью настоящей работы является реализация упрощенного метода фотонных карт для исполнения на графическом процессоре в реальном масштабе времени. В работе показано, каким образом можно дополнить трассировку лучей для корректного моделирования отражающих и прозрачных объектов, которые вызывают каустики - области с резко возрастающей интенсивностью светового поля. Основой программной реализации служат инструменты OpenGL/GLSL [6] и OpenCL [7], которые являются стандартами в области графики и гетерогенных вычислений и обеспечивают поддержку оборудования различных производителей. Благодаря широким возможностям, включающим запись в текстуру по произвольному адресу и работу с массивами данных на основе указателей, OpenCL является оптимальной платформой для задач сортировки и формирования ускоряющих структур. С другой стороны, максимальная производительность в задаче визуализации достигается с использованием шейдерного языка GLSL (OpenGL Shading Language), специально спроектированного для нужд компьютерной графики.
1. Моделирование каустик 1.1. Простейший подход
Для моделирования каустик может использоваться простой двухпроходный алгоритм. На первом проходе в соответствии с расположением и ориентацией источников света генерируются фотоны3, число которых определяется заранее. В простейшем случае можно рассматривать точечные (испускают фотоны равномерно по всем направлениям) или прямоугольные (испускают фотоны в произвольном направлении полупространства) источники света.
В данной работе использовался точеный источник. Однако для выделения значимых областей сцены выбирались только те направления, которые проходят через заданную прямоугольную площадку. Вместо использования генератора случайных чисел площадка покрывалась равномерной сеткой, для каждого узла которой генерировался фотон (рис. 1). При таком подхо-
3 В данном случае понятие фотона используется в контексте компьютерной графики: это виртуальная частица, переносящая некоторую порцию световой энергии и распространяющаяся по законам геометрической оптики.
де размерность сетки определяет общее число фотонов Ы, а фотонная карта может быть легко представлена в виде прямоугольной текстуры [8], каждый тексель которой соответствует точке регулярной сетки.
Рис. 1. Генерация фотонов для регулярных точек на прямоугольной площадке
Для каждого фотона прослеживается траектория его движения и взаимодействия с объектами сцены. При соударении с диффузным объектом в фотонную карту записываются координаты найденной точки соударения. При соударении с объектами, обладающими отражающими или прозрачными свойствами, отслеживаются дальнейшие взаимодействия в направлениях отражения или преломления соответственно. Для фотонов, покинувших сцену без взаимодействия с объектами, записываются недопустимые координаты (любые координаты вне ограничивающего параллелепипеда сцены). Интенсивность всех фотонов считается одинаковой и равной параметру 1рио0п- Подход может быть обобщен для учета потерь интенсивности фотона при соударении с неидеальными отражающими или прозрачными объектами (необходимо сохранять не только координаты точки соударения, но и интенсивность на момент данного соударения) или учета хроматических аберраций4 (необходимо генерировать фотоны отдельно для каждой длины волны).
При втором проходе используется традиционный алгоритм трассировки лучей. При этом к вычисленной в точке соударения освещенности прибавляется интенсивность фотонов из некоторой окрестности, доступных посредством обращения к сформированной фотонной карте. В
4Хроматические аберрации заключаются в паразитной дисперсии света, проходящего через оптическую систему (фотографический объектив, бинокль, микроскоп, телескоп и т.д.). При этом белый свет разлагается на составляющие его цветные лучи, в результате чего изображения предмета в разных цветах не совпадают в пространстве изображений.
данной работе принимается, что на освещенность точки влияют все фотоны из ее окрестности вместо к ближайших (как принято в традиционном методе фотонных карт), что в значительной степени упрощает процедуру сбора. Кроме того, принимается тривиальное соображение, что влияние фотона уменьшается с увеличением расстояния до точки соударения. В простейшем случае можно воспользоваться линейной весовой функцией, что приводит к следующей формуле:
||^ - А-II
1 1 ray tracing+~1 photon
^ maxi 0.1 --
+
Іршоп Z max{o. 1 - smoothstep(o, е
> P - Pi
в лексикографическом порядке. Данная сортировка позволяет применять алгоритм бинарного поиска со сложностью 0(1о§ N) для определения индексов фотонов с минимальной и максимальной координатой в кубической Е-окрест-ности точки соударения. В результате вычисления на втором проходе сокращаются до перебора фотонов из непрерывного сегмента карты, заключенного между данными индексами:
(1)
1 =1ray tracing*1 photon Z тЯХ\ 0 1 -
i=Nі [
||P - Pi\\
(3)
є
Здесь р - координаты обрабатываемой точки соударения, р1 - координаты г-го элемента фотонной карты. Параметр е - радиус окрестности сбора фотонов, который подбирается эмпирически для каждой конкретной сцены.
Согласно работе [2] сбор всех фотонов из фиксированной окрестности вместо к ближайших может привести к худшей оценке освещенности в областях с малой плотностью фотонов или размытому результату в областях с высокой плотностью. Тем не менее, за счет применения излагаемых далее быстрых алгоритмов выигрыш производительности оказывается весьма существенным и оправдывает данный подход. Кроме того, в качестве альтернативной весовой функции можно использовать интерполяцию Эрмита вида 3? -2?, ? - параметр в диапазоне [0, 1], отдавая больший приоритет близким фотонам и меньший - более отдаленным (в сравнении с линейной весовой функцией). Данный вид интерполяции реализован в ГПУ аппаратно и доступен через встроенную функцию 8шоо1Ь81ер языка шейдеров ОЬ8Ь. С использованием данной функции формула (1) принимает вид:
(2)
Применение модифицированной формулы (2) позволяет придать каустикам дополнительную резкость.
1.2. Сортировка фотонов
Описанный подход неэффективен, поскольку предполагает перебор всех элементов фотонной карты. Для повышения производительности можно использовать предварительную сортировку фононной карты по координатам фотонов
Здесь N\ и N2 - индексы фотонов с минимальными и максимальными координатами соответственно. Сортировка может быть выполнена средствами ЦПУ или ГПУ. В первом случае целесообразно использовать быструю сортировку со сложностью O(N log N), во втором - поразрядную или битоническую сортировку с трудоемкостью O(N) и o(n log2 N) соответственно.
Данный подход значительно повышает производительность визуализации, однако также не вполне эффективен: сбор фотонов осуществляется не из Е-окрестности точки соударения, а из Е-полосы вида | px — pix | <s . Далее будет рассмотрен алгоритм, позволяющий сократить перебор до некоторой окрестности точки соударения. Данный алгоритм развивает подход, предложенный в [3], обеспечивая более эффективное построение и использование ускоряющей структуры.
1.3. Воксельная ускоряющая структура
Разделим ограничивающий параллелепипед сцены на воксели со стороной е (радиус окрестности для сбора фотонов). Тогда каждому фотону с координатами (xp, yp, zp) можно сопоставить координаты вокселя, в котором он находится:
(4)
Х'p Xmin vy її Уp У min
є є
z — z •
p min
Здесь через (хтщ, Ушін, ^тіп) обозначена вершина ограничивающего параллелепипеда сцены с минимальными координатами.
После формирования фотонной карты описанным выше способом все фотоны сортируются в лексикографическом порядке по координатам соответствующих вокселей. При таком выборе оператора сравнения все фотоны, принад-
N
є
i=0
i=0
Xv =
zv
є
лежащие одному вокселю, являются эквивалентными и занимают непрерывный сегмент в отсортированной фотонной карте.
Для отсортированной фотонной карты возможно создание воксельной ускоряющей структуры (рис. 2). Для этого необходимо сформировать массив двумерных векторов, каждый элемент которого соответствует одному вокселю. В первую компоненту вектора записывается индекс первого фотона в соответствующем воксе-ле, во вторую - индекс последнего фотона. Для представления данного массива на ГПУ удобно использовать трехмерную текстуру.
Текстура вокселей
о
2
3
4
i-ч и
Ез Ь Ї5 Ї5
5
5
Текстура фотонной карты
P J- 0X PlX P J- 2X P 3X P 4X P5X P J- 6X P7X P J- sx
P J- 0Y PlY P J- 2Y P3Y P J- 4Y P5Y P6Y P7Y P J- SY
P J- 0Z PlZ P J- 2Z P3Z P J- 4Z P5Z P«Z Pyz P J- SZ
Рис. 2. Воксельная ускоряющая структура
Использование ускоряющей структуры на этапе визуализации основано на тривиальном соображении: все фотоны в є-окрестности точки соударения принадлежат вокселям, координаты которых отличаются не более чем на 1 от координат вокселя для данной точки (в противном случае разность координат составляет не менее длины вокселя, которая равна є). В итоге для вычисления освещенности в конкретной точке необходимо произвести суммирование по всем фотонам из 33 = 27 вокселей. Окончательная формула будет выглядеть следующим образом:
1 1 гау 1гаст^+
х, +1 у, +1 г, +1 С II р р |П
+ /„„„.І І І I ™х|о,1
(5)
=x„ - l j=У, -1 k = z, - l І = Nf‘
Здесь (хи у,, 2У) - координаты вокселя, которому принадлежит рассматриваемая точка, а и
М'{к - индексы соответственно первого и последнего фотона в данном вокселе. В большинстве случаев значимыми окажутся лишь несколько вокселей из окрестности, в то время как пустые воксели будут быстро отброшены.
1.4. Формирование массива вокселей
Для типичной сцены значение параметра е лежит в диапазоне 0.5-1.5% от размера сцены (в качестве которого примем наибольшую сторону ограничивающего параллелепипеда). Такой выбор параметра е порождает множество вокселей размером не более 2563 элементов, что делает возможным создание и хранение данной ускоряющей структуры на ГПУ. Кроме того, в ряде случаев можно ограничиться двумерным массивом вокселей для моделирования каустик на плоской поверхности (например, на полу и стенах комнаты), что сокращает размер массива на два порядка. Для вычисления индекса первого и последнего фотона в каждом вокселе предлагается следующий параллельный алгоритм.
Для каждого элемента отсортированной фотонной карты выполняется проверка, является ли данный элемент первым в соответствующем вокселе (не лежит ли предыдущий элемент в другом вокселе). В случае положительного результата номер элемента записывается в соответствующий воксель в качестве индекса первого элемента. Аналогичная процедура выполняется для определения индекса последнего элемента вокселя. Поскольку в каждом вокселе существует единственный первый и последний элемент, не возникнет никаких коллизий при параллельном выполнении указанных операций. Данный алгоритм эффективно выполняется на ГПУ и имеет линейную сложность O(N ), где N - число генерируемых фотонов.
1.5. Схема результирующего алгоритма
В результате получаем следующий четы-рехпроходный алгоритм.
• Генерация фотонной карты с помощью прямой трассировки лучей. На ГПУ фотонная карта представлена двумерной текстурой размера PX х PY, каждый тексель которой соответствует ровно одному выпущенному фотону. Общее число фотонов N = PX х PY. Для реализации данного этапа использовался шейдерный язык GLSL.
• «Воксельная» сортировка фотонной карты на ГПУ. Эффективными алгоритмами являются битоническая сортировка со сложностью o(n log2 N) и поразрядная сортировка со сложностью o(N ). В данной работе реализована битоническая сортировка на базе OpenCL.
• Формирование массива вокселей для отсортированной фотонной карты на ГПУ. Пред-
ложенный алгоритм решает данную задачу за линейное время о(Ы). Для реализации этого этапа использовался ОреиСЬ.
• Визуализация сцены с помощью обратной трассировки лучей. Для обработки каустик выполняется сбор фотонов из окрестности точки соударения предложенным выше способом. Наряду с первым проходом для реализации выбран шейдерный язык ОЬ8Ь. Оба этапа используют общие шейдеры с минимальными изменениями.
Предлагаемый алгоритм обеспечивает эффективную визуализацию каустик, выполняя подготовку необходимой ускоряющей структуры за линейное время при использовании поразрядной сортировки.
Необходимо сделать несколько замечаний относительно программной реализации. Текущая версия стандарта ОреиСЬ не позволяет ядру выполнять чтение и запись в одну и ту же текстуру, что необходимо на этапе воксельной сортировки. Для обхода данного ограничения текстура с фотонной картой копируется в буфер (простой одномерный массив) перед этапом сортировки, содержимое которого вновь копируется в текстуру после окончания сортировки. На этапе формирования массива вокселей запись первого и последнего фотона осуществляется в различные компоненты текселя. Тем не менее, стандарт не позволяет обновлять отдельные компоненты текселя: записи подлежит сразу весь вектор цвета. Для обхода ограничения данный этап разбивается на два прохода: при первом проходе определяются индексы первых фотонов, а при втором - последних.
1.6. Комбинирование возможностей OpenCL и GLSL
Как следует из рассмотренной схемы, OpenCL используется только для программирования этапов, связанных с формированием ускоряющей структуры, в то время как этапы визуализации реализованы на основе OpenGL/GLSL. Формирование ускоряющей структуры требует записи в текстуру по произвольному адресу, что недоступно в шейдерном языке GLSL. С другой стороны, визуализацию удобнее и эффективнее выполнять на специализированном языке GLSL.
Для проверки данного тезиса было разработано несколько простых примеров трассировки лучей на базе GLSL и OpenCL. В процессе программирования OpenCL-версий учитывались все методики, изложенные в руководстве по оптимизации [9]. Для карты AMD Radeon 5870 1 Гб (R800) и ATI Stream SDK 2.1 реализация на основе GLSL оказалась эффективнее на 250% и 30% при визуализации неявно заданной поверхности и фрактала соответственно (рис. 3). В задаче построения фрактала преобладают более сложные вычисления, что существенно меняет соотношение производительности, хотя приоритет остается за шейдерной реализацией.
Таким образом, для извлечения максимальной выгоды целесообразно комбинировать оба инструмента. Данный подход приветствуется и самим стандартом OpenCL: эффективное взаимодействие двух интерфейсов возможно в рамках общего контекста OpenCL/OpenGL.
Рис. 3. Трассировка неявно заданной поверхности и 4D-фрaктaлa Жюлиа
2. Оценка производительности
Оценка производительности рассмотренных алгоритмов проводилась на примере визуализации неявно заданной поверхности - «Metaballs» или «Blobs» [10]. Поверхность визуализировалась непосредственным образом без предварительной триангуляции. При расчете учитывался эффект полного внутреннего отражения, глубина трассировки - 5 (рис. 4).
Для проведения эксперимента использовались графические ускорители AMD Radeon
5870 (R800) и NVIDIA Quadro FX 5600 (G80). Карта AMD работала под управлением драйвера Catalyst 10.4 и ATI Stream SDK 2.1, карта NVIDIA - под управлением драйвера Force Ware 197.13 и CUD A Toolkit 3.0. Для тестовой сцены ограничивающий параллелепипед имеет размеры 10 х 10 х 10, а значение параметра є = 0.04 = 0.4% от размера сцены. Визуализация выполнялась в разрешении 1920 х 1080 (Full HD). В следующих таблицах представлены результаты замера производительности.
Таблица І
Лексикографическая сортировка фотонов
Число фотонов Лексикографическая сортировка (мс) Производительность визуализации (к/с)
R5870 Q5600 R5870 Q5600
1282 = 16K 10.2 12.5 20 7.5
2562 = 65K 14.5 24.2 8 3.2
5122 = 262K 28.5 60.1 3 1.0
Таблица 2
Воксельная ускоряющая структура
Число фотонов Воксельная сортировка и заполнение массива вокселей (мс) Производительность визуализации (к/с)
R5870 Q5600 R5870 Q5600
1282 = 16K 10.5 16.8 49.5 19.5
2562 = 65K 17.5 30.3 35.8 14.5
5122 = 262K 39.4 98.5 18.2 6.5
10242 = 1M 140 381 5.2 1.8
Результаты экспериментов позволяют сделать вывод о высокой эффективности предложенного метода. Воксельная ускоряющая структура более чем в четыре раза быстрее метода на основе лексикографической сортировки и обеспечивает качественную визуализацию каустик в реальном времени с разрешением Full HD.
Заключение
Предложен вариант метода фотонных карт для выполнения на ГПУ в реальном масштабе времени. Он позволяет корректно моделировать каустики, возникающие при взаимодействии света с прозрачными или отражающими объектами. Эффективность визуализации достигается за счет воксельной ускоряющей структуры, которая формируется непосредственно на ГПУ средствами OpenCL. Формирование структуры включает в себя «воксельную» сортировку фотонов и определение для каждого вокселя индекса первого и последнего фотона. Данные задачи могут быть решены с помощью алгоритмов со сложностью o(n ). Предложенный подход обеспечивает высокую скорость работы при использовании фотонных карт вплоть до 1 миллиона элементов с разрешением Full HD.
Список литературы
1. Rushmeier H., Banks D., Shirley P. A Basic Guide to Global Illumination. SIGGRAPH 98 Course 5.
2. Jensen H.W. Realistic Image Synthesis Using Photon Mapping. AK Peters. Ltd., Massachusetts. 2001.
3. Purcell T.J., Donner C., Cammarano M., Jensen H.W. and Hanrahan P. Photon mapping on programmable graphics hardware // Proc. ACMSIGGRAPH/ EUROGRAPHICS Conference. 2003. Р. 41-50.
4. Fleisz M. Photon Mapping on the GPU // Master's thesis, School of Informatics, University of Edinburgh, 2009.
5. Fabianowski B. and Dingliana J. Interactive Global Photon Mapping // Computer Graphics Forum. 2009. V. 28, Iss. 4. P. 1151-1159.
6. Kessenich J., Baldwin D., Rost R. The OpenGL Shading Language (Revision 8, 7 September 2006). http://www.opengl.org/registry/doc/GLSLangSpec.Full.! .20.8.pdf
7. Khronos OpenCL Working Group. The OpenCL Specification (Version 1.0, Document Revision 29, Last Revision Date 12/8/08).
http://www.khronos.org/registry/cl/specs/opencl-L0.29.pdf
8. Спецификации расширения ARB_texture_rec-tangle (Revision 1.21, 4 October 2005). http://www. opengl.org/registry/specs/ARB/texture_rectangle.txt
9. ATI Stream Computing OpenCL Programming Guide (Revision 1.05, August 2010). http://developer. amd.com/gpu/ATIStreamSDK/assets/ATI_Stream_SDK _OpenCL_Programming_Guide.pdf
10. Ryan Geiss. Metaballs (Blobs) (3/10/2000). http://www.geisswerks.com/ryan/BLOBS/blobs.html
REAL-TIME CAUSTICS SIMULATION BASED ON COMBINED CAPABILITIES OF OpenCL AND
SHADERS
D.K. Bogolepov, V.E. Turlapov
The objective of this paper is to adapt photon mapping to real-time caustics simulation. A high-performance algorithm adapted to the GPU and implemented on the basis of cross-platform OpenGL and OpenCL APIs is proposed. OpenGL shaders are used for high-performance caustics rendering and building the photon map. Voxel acceleration structure constructed directly on the GPU by means of OpenCL provides the effective use of the photon map. Performance estimation for different equipment is given.
Keywords: interactive rendering, ray tracing, photon mapping method, caustics simulation, GPGPU, OpenGL, GLSL, OpenCL.
