CC BY
72
Известия Тульского государственного университета Естественные науки. 2012. Вып. 3. С. 142-152
= ИНФОРМАТИКА :
УДК 519.711.3; 519.40
Применение GPU для вычисления коэффициентов затенения в задаче моделирования интерактивного освещения
С.С. Свистунов
Аннотация. Рассматривается задача моделирования интерактивного первичного освещения на основе уравнения рендеринга Kaijia. Для учета функции видимости используются варианты Ambient occlusion. Предлагается способ быстрого вычисления коэффициентов затенения при помощи возможностей современных GPU.
Ключевые слова: уравнение рендеринга, низкочастотное
освещение, функция видимости, ДФОС, сферическая гармоника, сферический дизайн, коэффициент затенения, графический процессор, буфер глубины, CUDA, OpenCL.
Введение
Эта работа продолжает исследования, начатые в [1, 2].
Изучается актуальная задача ЭБ-графики, заключающаяся в быстром вычислении коэффициентов затенения непрозрачной модели. Она возникает в задаче моделирования интерактивного первичного освещения на основе уравнения рендеринга (интеграла освещенности) Кау1а [3]:
Бр(у) = Цх)Ур(х)рр(х,у) йц,(х). (1)
Js2
Здесь Б2 — единичная сфера, й^(х) = йх/(4п) — нормированная мера на сфере, Вр(ь) — функция яркости отраженного освещения из точки модели р в направлении единичного вектора V € Б2, Ь(х) — функция яркости удаленного окружающего освещения (свет), Ур(х) — функция видимости, рр(х^) — двунаправленная функция отражающей способности (ДФОС или БИБР) с учтенным ламбертовским множителем (прх)+ = тах(0, прх), где пр — нормаль поверхности модели в точке р. Далее обозначение точки р часто будем опускать.
Работа состоит из теоретической и прикладной частей. В первой из них приводятся способы аппроксимации интеграла освещенности (1)
для рендеринга в реальном времени (real-time rendering) и приводятся два определения коэффициента затенения, включая общеизвестный Ambient occlusion [4, 5]. Вторая часть посвящена методам вычисления коэффициента затенения при помощи возможностей современных графических процессоров (GPU). Используются общеизвестные факты из 3Б-графики, гармонического анализа на сфере и программирования для GPU (см. обзор [3] и библиографический список в нем, а также [5-9]).
Свет L. Рассматривается функция L: S2 ^ R+, для которой
\\L — LW\ ^ £IIlII2> (2)
где L(x) = i=o^2m=-i LimVim(x) — частичная сумма сферического ряда Фурье функции L, Vim — ортонормированные действительные сферические гармоники, е — малое положительное число, ||L||2 — L2(S2)-норма функции L (ее квадрат ||L||2 имеет смысл энергии L).
Если порядок частичной суммы s = deg L — небольшое целое число, то освещение называется низкочастотным (основная энергия L передается низкими частотами). Оно очень часто встречается на практике.
В приложениях освещение, как правило, задается значениями на густой дискретной сетке (например, картой окружения). Для вычисления в промежуточных точках используется билинейная интерполяция.
BRDF р. Функция pp(x,v) характеризует свойства материала поверхности модели по рассеянию в направлении v, приходящего в точку p с направлений x окружающего освещения. В интерактивном рендеринге различают диффузный (diffusion), зеркальный (specular) и анизотропный (anisotropic) случаи.
В диффузном случае p(x,v) = p(nx). Поскольку здесь нет зависимости от v, то и освещение B не зависит от v. В этом случае свет в точке рассеивается одинаково по всем направлениям. Основной диффузной BRDF является BRDF, отвечающая закону Ламберта: p(x) = 4(nx)+, где n — нормаль, значение 4 получается из условия fs2 р dц = 1. BRDF здесь приводится без множителя. На практике она умножается на коэффициент поглощения Cd (альбедо поверхности).
В зеркальном случае свет, приходящий в точку, рассеивается в концентрировано в малом телесном угле, зависящем от направления просмотра v (глаз это видит в виде блика на поверхности): p(x, v) = p(rv), где r = 2(nx)n — x — зеркальный к x вектор. Основным примером зеркальной BRDF, активно используемой в компьютерной графике (несмотря на ее физическую необоснованность), является BRDF Фонга p(x) = (2а + 2)(rv)+, где а > 0 — показатель Фонга (чем он больше, тем уже блик).
Если поверхность обладает диффузными и зеркальными свойствами, то ее BRDF приближенно можно записать линейной комбинацией ламбертовской и фонговской BRDF. BRDF для анизотропного случая в [8].
Функция видимости V. Для непрозрачной модели функция видимости Vp(x): S2 ^ {0,1} определяется следующим образом: она равна
0, если луч, выпущенный из точки модели p по направлению x пересекает модель и 1 иначе. Считаем, что Vp(x) = 0 при npx ^ 0.
Если модель выпуклая, то в интеграле освещенности можно опустить функцию видимости и его аппроксимация значительно упрощается. В общем случае функция видимости может быть очень сложной и ее точный учет становится проблематичным.
Для сведения к выпуклому случаю воспользуемся полиномиальной аппроксимацией функции видимости
V(x) = V0 + Vix + V2x ■ x + ...,
где V0 — константа, V1 — 3-вектор, V2 — 3-матрица. Отметим, что здесь ряд не обязательно отвечает сферическому ряду Фурье функции V(x), поскольку могут использоваться разные метрики аппроксимации.
На практике, как правило, основную часть энергии функции видимости V несет коэффициент V0, называемый коэффициентом затенения. Положим V(x) ~ V0. В этом случае из уравнения (1) (как в теореме о среднем значении) получаем
/ L(x)V(x)'p(x,v) d^(x) V0 L(x)p(x,v) d^(x). (3)
J S2 Jnx>0
1. Рассмотрим диффузный случай и в качестве метрики аппроксимации потребуем точности формулы (3) на постоянном освещении L = const. Тогда получим известную формулу для вычисления Ambient occlusion [4]
V0 = ^ V (x)(nx)+ d^(x).
S2
В выпуклом случае Vo = 1.
2. Потребуем точности формулы (3) при Lf) = const. Тогда получаем следующую известную формулу для коэффициента затенения
A = V0 = % V(x) d^(x).
S2
В выпуклом случае, по-прежнему, Vo = 1. Коэффициент затенения A зависит только от внутренней геометрии модели.
Сферические кубатурные формулы. При моделировании освещения на основе уравнения рендеринга возникает задача приближенного вычисления интеграла If = fs2 f (x)u(x) d^(x) с неотрицательным весом u(x), в котором f может быть как гладкой функцией (например,
сферическим полиномом), так и функцией, задаваемой дискретно. Применим для этого метод кубатурных формул, когда
N
If ~ Qf = Y^ Xkf (Xk), k=l
где Xk > 0 и Xk E S2 — соответственно веса и узлы кубатурной формулы. Потребуем точности для f = const, тогда £N=l Xk = и0 = fs2 и d^. Если и = = 1, то и0 = 1. В случае равных весов Xk = u0/N.
Приведем два способа выбора весов.
1. Метод Монте-Карло для равных весов и и = 1. В этом случае точки Xk
выбираются равномерно распределенными на сфере S2. В [10], например, генерируются два независимых равномерно распределенных на (—1,1) числа Cl, $2, для которых s = ^2 + < 1, и полагается
Xk = (2$i(1 — s)l/2, 2$2(1 — s)l/2,1 — 2s).
Если f E C(S2), то справедлива оценка погрешности
\If — Qf I = O(N-l/2).
Из нее следует, что для получения точности порядка 10-2 нужно брать порядка 104 точек Xk.
2. Если функция f является достаточно гладкой, то более выигрышным по сравнению с методом Монте-Карло является метод сферических дизайнов (или интерполяционных кубатурных формул). Множество узлов и весов называется сферическим s-дизайном, если формула If = Qf точна для всех сферических полиномов f порядка s [11, 12]. Дизайн называется минимальным, если он при заданном порядке s содержит наименьшее число узлов N. Нахождение минимальных дизайнов является сложной математической задачей. В случае и = 1 и равных весов дизайны, близкие к минимальным, можно найти в [11], а для неравных весов они построены в серии работ В.И. Лебедева и его последователей (см., например, [13, 14]).
Для гладких функций использование лебедевских дизайнов при одинаковой точности позволяет на порядок сократить число узлов по сравнению с методом Монте-Карло.
При использовании сферических кубатурных формул для коэффициентов затенения Ap, зависящих от точки модели р, приближенно имеем
N
Ap & 2 XkVp(Xk). (4)
k=1
Методы PRT и DPRT. В заключении первой части работы приведем два способа аппроксимации интеграла освещенности (1). Первым из них является известный метод PRT (precomputed radiance transfer) [3, 6, 7]. Он использует разложения функций по сферическим гармоникам и особенно
эффективен в диффузном случае. Второй способ предложен в работе автора [1] и называется DPRT (design precomputed radiance transfer). Он использует взвешенные дизайны небольшого порядка и эффективен для зеркального и анизотропного случаев.
В обоих случаях рассматривается низкочастотный свет L. Если функция видимости учтена в виде коэффициента затенения, то
B(v) и AB(v), B(v) = L(x)p>(x,v) dp.(x).
JS2
В диффузном случае B(v) = В не зависит от v, но зависит от нормали n в точке р:
/» Л S l
B(n)= L(x)p>(nx) dp,(x) и / L(x)p(nx) dp,(x) = EE Ll mTlm (n),
JS2 ,'S2 l=0 m=-l
где Tlm(n) = (2l + 1)-l/2ploVlm(n) — коэффициенты функции транспорта. Эти коэффициенты один раз рассчитываются в каждой вершине модели р и сохраняются. Также один раз рассчитываются сферические коэффициенты Фурье света Lim. Если впоследствии модель вращается (без изменения формы), то коэффициенты Llm легко перерасчитываются исходя из матрицы повората, а Tim не меняются. В этом состоит идея метода PRT.
Метод PRT испытывает сложности в зеркальном случае, когда
,, s l
B(v) = L(x)p(rv) dp,(x) и^ ^2 LimTlm(n,v).
JS2 l=0 m=-l
В этом случае коэффициенты транспорта дополнительно зависят от заранее не известного вектора v. Метод же DPRT позволяет выполнять расчеты так же эффективно, как и в случае диффузного PRT. Пусть {xk, Xk}N=l — s-дизайн, построенный по весу и(х) = р(хз) (для малых s значение N также мало, подробности см. в [1]). Тогда
~ f ~ i ~ N'
B(v) и / L(x)p(rv) dp.(x) = L(gnvх)р(х3) dp.(x) = V XkL(gnvXk),
Js2 Js2 k=l
где gnv — вращение, переводящее вектор v в n. Значения L(gnvXk) для каждой вершины модели вычисляются эффективно [1].
Вычисление коэффициентов затенения вершин модели
при помощи ОЕИ
Рассмотрим задачу практического вычисления коэффициентов затенения точек трехмерной модели, задаваемой полигональным мешем, по формуле (4). Используется дизайн нужного порядка, например, один
из лебедевских [13], самый большой из которых имеет порядок 131 и состоит из 5810 точек. Тогда для вычисления в каждой вершине модели p коэффициента затенения Ap необходимо найти значения функции видимости Vp(xk) для k = 1,..., N. Это требует определения факта пересечения луча, выпущенного из вершины p по направлению Xk, с мешем модели.
При наивном подходе данная задача сводится к полному перебору полигонов модели. Его выполнение на центральном процессоре (CPU) для моделей с сотнями тысяч вершин и полигонов занимает продолжительное время. Поэтому для ускорения расчетов на CPU используют ускоряющие структуры, основанные, например, на BVH-деревьях, пространственном хешировании и т.п. Это позволяет на порядок снизить временные затраты, но пока не позволяет приблизиться к real-time из-за ограничений CPU.
Мы преследуем цель вычисления коэффициентов затенения Ap (и близких задач) в реальном времени. Для этого естественно воспользоваться возможностями современных GPU, включая буферы глубины (Z-buffer), CUDA и OpenCL. Использование Z-buffer позволяет быстро выполнить операцию определения пересечения луча и меша. Применение технологий CUDA и OpenCL позволяет выполнять вычисления коэффициентов затенения в параллельном режиме очень быстро, приближаясь к real-time.
Вначале кратко опишем применяемые нами технологии.
Используемые технологии. Во всех приложениях для работы с 3D-графикой использовался OpenGL — открытая графическая библиотека (графическое API).
Для расчета удаленного окружающего освещения в интерактивной сцене, в том числе с учетом функции видимости, было реализовано приложение на языке программирования Python с использованием Panda3D — фреймворка для трехмерной визуализации и разработки игр. В этом приложении для вычисления функции видимости использовалась кубическая карта окружения из буферов глубины (см. далее). Это привело к значительному ускорению расчетов. Приложение представляет из себя демонстрационную реализацию метода DPRT. При реализации шейдеров (программа для одной из ступеней графического конвейера) использовался язык шейдеров GLSL.
Для более полного использования возможностей GPU был реализован плагин к Blender — свободному пакету для создания трехмерной компьютерной графики, включающему в себя средства моделирования, анимации, рендеринга, постобработки видео, а также создания интерактивных игр. В качестве языка программирования использовался Python. Для вычислений на GPU использовалась технология NVIDIA CUDA (программно-аппаратная архитектура параллельных вычислений), а именно PyCUDA — библиотека для использования CUDA на языке Python. Плагин позволяет вычислять и отображать значения функции видимости для моделей используемых в Blender.
Использование Z-buffer. Для ускорения расчета факта пересечения луча и меша воспользуемся возможностью GPU сохранять буфер глубины. Представим, что в каждой вершине модели расположено шесть камер, ориентированных по нормали и направленных по соответствующим координатным осям. При помощи каждой камеры можно снять (отрендерить) буфер глубины. Таким образом, получаем кубическую карту окружения со значениями глубины в данной вершине модели. Важно, что эти операции на GPU за счет параллелизма выполняются очень быстро.
Как известно, карта окружения — это изображение 3D-сцены в виде текстуры. В случае кубической карты окружение такая текстура представляет собой развертку куба, на каждой грани которого запечатлена окружающая 3D сцена.
Теперь для поиска пересечений луча с моделью достаточно найти необходимую сторону карты окружения и с помощью простых матричных преобразований определить пиксель, отвечающий лучу. Далее по цвету этого пикселя в текстуре буфера глубины несложно определить, было пересечение или нет (если цвет не черный, то пересечение есть). Это процесс похож на известный метод построения теней от точечных источников света — теневыми картами (Shadow mapping).
Данный подход позволяет избежать вложенного цикла по полигонам при поиске пересечений луча с мешем и заменить его созданием карты глубины с помощью стандартных средств OpenGL/DirectX.
Использование CUDA и OpenCL. С появлением технологий NVIDIA CUDA и OpenCL стало гораздо проще использовать параллельные возможности GPU. Для выполнения расчетов с использованием CUDA использовались различные способы распараллеливания вычислений: по вершинам модели (ядро с одномерной топологией), по вершинам и лучам (ядро с двумерной топологией), по полигонам и лучам (ядро с двумерной топологией). Наиболее удобным и результативным подходом оказалось распараллеливание по вершинам и лучам.
При реализации возникла проблема ограничения максимального времени выполнения одного запуска ядра (около одной секунды). Для обхода этого ограничения вычисления разбивались на несколько этапов: обработка по несколько лучей, обработка частей модели. Стоит уделять особое внимание работе с различными типами памяти (особенно на бюджетных видеокартах). Самый простой, но наименее эффективный способ реализации — это трансляция алгоритма однопроцессорной версии и использование глобальной памяти видеокарты. Такой подход дает ускорение в расчетах, но простая замена глобальной памяти на текстурную позволяет в разы улучшить производительность.
Предпринимались попытки использования различных типов памяти видеокарты — от глобальной и текстурной до реализации управляемого кэша в разделяемой памяти и побитового хранения результатов вычислений.
При использовании разделяемой памяти в ней сохранялась информация о полигонах, что позволило (как и положено) снизить количество чтений из глобальной памяти и ускорить расчет. Однако наибольшей производительности удалось достичь при использовании разделяемой памяти и битового формата хранения значений функции видимости.
Реализация. Рассмотрим версию трансляции алгоритма однопроцессорной версии. Результат представляет собой float Vis[N_v], где N_v количество вершин; используется одномерная blockDim(128), gridDim(N_v/128) :
// нить -- это вершина
int tid = threadIdx.x + blockIdx.x*blockDim.x sum = 0
// цикл по всем полигонам модели for i=0..Np:
// отсеиваем явно противоположное направление луча if (dot(Normal[tid], Direction) > 0.0):
sum += weight * intersect(Polygons[i], Vertex[tid], Direction) Vis[tid] = sum
Даже такая реализация дает значительное ускорение при вычислении функции видимости.
Для ускорения расчетов можно воспользоваться разделяемой памятью и закешировать полигоны модели:
— результат: bool Vis[N_v*N_d], где N_d — количество лучей;
— топология: одномерная, blockDim(128), gridDim(N_p/128), где N_p — количество полигонов;
— количество запусков равно количеству блоков в сетке;
— в разделяемую память загружаем 128 полигонов на блок;
— при следующем запуске делаем циклический сдвиг для подгружаемых полигонов;
— отдельное ядро для получения конечного результата: среднее значение по всем лучам из каждой вершины.
// нить -- полигон
int tid = threadIdx.x + blockIdx.x*blockDim.x
SharedPolygons[threadIdx.x] = Polygons[pi%Nv]
__syncthreads()
// на каждый полигон по три вершины for j=0..3:
// цикл по направлениям for di=0..Nd:
// если уже было пересечение то продолжаем
if (!Vis[i(tid, di)]) continue // цикл по подгруженным полигонам for pi=0..blockDim.x:
Vis[i(tid, di)] &&= intersect(SharedPolygons[pi], Vertex[vi(tid, j)], Direction[di]);
Для достижения большего ускорения можно кешировать в разделяемой памяти дополнительную информации, а результат упаковывать в битовом представлении, используя свойство функция видимости принимать только два значения 0 или 1:
— результат: uint16 Vis\_bits[N_v*N_d/16], запись ведется побитово;
— топология: многомерная; blockDim(16,16), gridDim(N_d/16, P*16), где, например, P = 10;
— топология: по x считываются направления, по y полигоны и вершины;
— запуск ядра осуществляется в двойном цикле для перебора всех вершин и полигонов:
verts_steps = Nv/(blockDim.y*gridDim.y) poly_steps = Np/(blockDim.y*P) for v_index = 0..verts_steps: for p_index = 0..poly_steps: run kernel
— каждая нить по y грузит в разделяемую память P полигонов (всего 160 полигонов); каждый блок грузит одни и те же полигоны; перебор осуществляется при очередном запуске ядра.
// Кеширование в разделяемой памяти
int tx = threadldx.x + blockDim.x*blockIdx.x; int ty = threadIdx.y + blockDim.y*blockIdx.y;
// индекс вершины
int vi = ty + blockDim.y*gridDim.y*step_vert;
// индекс полигона
int pi = 3*P*(threadIdx.y + step_poly*blockDim.y);
if (threadldx.x == 0): for ii = 0..P:
SharedPolygons[i(threadIdx.y, ii)] = Polygons[i(pi, ii)]; SharedIntersect[threadIdx.y] = Vis_bits[i(vi, blockIdx.x)]; SharedVertex[threadIdx.y] = Vertex[vi];
SharedNormals[threadIdx.y] = Normals[vi];
__syncthreads();
// расчет // текущий бит
unsigned short int current_bit = (1 << threadIdx.x)
// проверка бита и направления
if (current_bit & SharedIntersect[threadIdx.y]) == 0 && dot(Direction, Normal) > 0.0f):
// цикл по полигонам for i = 0..P*blockDim.y:
if (!intersect(SharedPolygons[i], SharedVertex[threadIdx.y], Direction[i(tx)])):
SharedIntersect[threadIdx.y] |= current_bit; break;
__syncthreads();
if (threadIdx.x == 0 && SharedIntersect[threadIdx.y] != 0): Vis_bits[i(vi,blockIdx.x)] |= SharedIntersect[threadIdx.y];
Результаты. Была подготовлена тестовая сцена, состоящая приблизительно из 5к вершин и полигонов. Использовался бюджетный компьютер с видеокартой NVIDIA. Результаты приведены в следующей таблице: _______________________________________
Метод Время расчета АО
CPU более 1 часа
Буфер глубины около 11 сек
NVIDIA CUDA около 0.8 сек
Они подтверждают наши выводы.
Список литературы
1. Свистунов С.С. Интерактивный рендеринг при помощи сферических дизайнов (SDPRT) для низкочастотного окружающего освещения и BRDF Фонга // Изв. ТулГУ. Естественные науки. 2011. Вып. 1. С. 188-199.
2. Горбачев Д.В., Иванов В.И., Странковский С.А. Моделирование освещения в интерактивной графике при помощи сферических дизайнов // Изв. ТулГУ. Естественные науки. 2007. Вып. 1. С. 17—36.
3. Ramamoorthi R. Precomputation-Based Rendering // Foundations and Trends in Computer Graphics and Vision. 2009. V. 3, № 4. P. 281-369. http://www.cs.berkeley. edu/~ravir/
4. Kontkahen J, Samuli L. Ambient occlusion fields // http://www.tml.tkk.fi/~janne/ aofields/aofields.pdf
5. Shanmugam P., Arikan O. Hardware accelerated ambient occlusion techniques on GPUs // http://perumaal.googlepages.com/ao.pdf
6. Green R. Spherical Harmonic Lighting: The Gritty Details // http://www.research. scea.com/gdc2003/spherical-harmonic-lighting.html
7. Sloan P., Kautz J., Snyder J. Precomputed Radiance Transfer for Real-Time Rendering in Dynamic, Low-Frequency Lighting Environments // http://web4.cs.ucl.ac. uk/staff/j.kautz/publications/prtSIG02.pdf
8. Ramamoorthi R. A signal-processing framework for forward and inverse rendering: Ph. D. dissertation. Stanford University. 2002.
9. Akenine-Moller T, Haines E, Hoffman N. Real-Time Rendering. Massachusetts: Wellesley, 2008.
10. Marsaglia G. Choosing a Point from the Surface of a Sphere // Ann. Math. Statist. 1972. V. 43, № 2. P. 645-646.
11. Hardin R.H., Sloane N.J.A. McLaren’s Improved Snub Cube and Other New Spherical Designs in Three Dimensions // Discrete and Computational Geometry. 1996. V. 15. P. 429-441. http://www.research.att.com/~njas
12. Мысовских И.П. Интерполяционные кубатурные формулы. М.: Наука, 1981.
13. Лебедев В.И., Лайков Д.Н. Квадратурная формула для сферы 131-го алгебраического порядка точности / / Докл. РАН. 1999. Т. 366, № 6. С. 741—745; http://www.mathworks.com/matlabcentral/fileexchange/27097-getlebedevsphere
14. Попов А.С. Кубатурные формулы на сфере, инвариантные относительно группы вращений икосаэдра // Сиб. журн. вычисл. матем. 2008. Т. 11, № 4. С. 433—440.
Свистунов Сергей Сергеевич ([email protected]), аспирант, кафедра прикладной математики и информатики, Тульский государственный университет.
Application for GPU computing shading coefficient in the problem modeling of interactive lighting
S.S. Svistunov
Abstract. The problem of modeling of interactive primary lighting by rendering Kaijia equation is considered. Ambient occlusion variants are used to account visibility function. The method of quick calculation of shading coefficients with the help of modern GPU’s capabilities is suggested.
Keywords: rendering equation, low-frequency light, visibility function, BRDF, spherical harmonics, spherical design, Ambient occlusion, GPU, Z-buffer, CUDA, OpenCL.
Svistunov Sergey ([email protected]), postgraduate student, department of applied mathematics and computer science, Tula State University.
Поступила 14-10.2012
