Интерактивный рендеринг при помощи сферических дизайнов (sdprt) для низкочастотного окружающего освещения и BRDF Фонга Текст научной статьи по специальности «Математика»

Известия Тульского государственного университета Естественные науки. 2011. Вып. 1. С. 188-199

Информатика

УДК 004.925

Интерактивный рендеринг при помощи сферических дизайнов (8ВРКГ) для низкочастотного окружающего освещения и БИБЕ Фонга

С.С. Свистунов

Аннотация. Предлагается новый метод интерактивного рендеринга SDPRT для низкочастотного окружающего освещения.

Он дает результаты, аналогичные получаемым по методу PRT.

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

Ключевые слова: интерактивный рендеринг, сферический

дизайн, освещение, низкочастотное окружающее освещение, precomputed radiance transfer (PRT).

Введение. В настоящее время актуальной является задача физически обоснованного интерактивного (real-time) рендеринга сцены, состоящей из модели и окружающего освещения, задаваемое Cubemap [1]. Методы, применяемые в неинтерактивном (offline) режиме, например, трассировка лучей по методу Монте-Карло, пока не пригодны для применения в современных GPU в силу большой вычислительной сложности [2].

Изучается случай первичного освещения. Для решения поставленной задачи был разработан метод PRT (Precomputed Radiance Transfer) [3]. Он основывается на аппроксимации уравнения освещенности Kaijia [1, с. 15-25]:

Здесь Вр(ь) — отраженное освещение из точки модели р по направлению V е Б2; Б2 — единичная сфера, йц,(х) = ^х/(4п); Ь(х) — окружающее

1. Метод PRT

(1)

освещение; Vp(x) — функция видимости в точке p; P>p(x, v) — BRDF в точке p с учтенным ламбертовским множителем (nx)+ = max(0,nx), где n = np — нормаль в точке p. Для непрозрачных поверхностей Vp(x) = 0 при nx ^ 0.

Рассматривается низкочастотное (е-низкочастотное) освещения L, для которого

\\L — L||2 ^ e||L||2, (2)

где

8 l

L(x) = Y.Y. Llmylm(x) l=0 m=-l

— частичная сумма ряда Фурье функции L, е — малое положительное число (в приложениях часто полагают е < 0.1). Если использовать в (1) вместо L ее частичную сумму L и воспользоваться неравенством Коши-Буняковского, то имеем

s l

Bp(v) Bp(v) = ЕЕ LlmTplm(v) (3)

l=0 m=-l

с абсолютной погрешностью

\bp(v) — BP(v)l < e\\Lh\\pp{ • ,v)\\2 (4)

(она оценивает качество итогового изображения). Здесь

Tplm(v) = ylm(x)Tp(x,v) d^(x), Tp(x,v) = Vp(x)pp(x,v). (5)

JS2

Функция Tp(x,v) называется функцией транспорта (transfer), Tplm(v) — ее коэффициенты Фурье по зональным сферическим гармоникам.

В интерактивном рендеринге различают случаи диффузного (diffusion), зеркального (specular) и анизотропного (anisotropic) отраженного освещения.

В диффузном случае BRDF для каждой точки p имеет вид: p)(x, v) = = Cd(nx)+ (зависимость от p понятна и опускается). Она не зависит от направления v. Следовательно, коэффициенты Tplm также не зависят от v. Их можно вычислить для каждой вершины p в offline (на CPU) и сохранить. Интегралы в offline вычисляются методоми Монте-Карло или дизайнов [4]. В низкочастотном случае используются значения s = 2, 3, 4. В этом случае в (3) достаточно просуммировать 9,16,25 слагаемых соответственно. Эта операция реализуется в шейдере. Если модель или Cubemap вращаются, то коэффициенты Llm для каждого кадра «подкручиваются» на CPU [5]. Для сжатия текстур применяется метод LPCA, основанный на SVD-разложении. В этом состоит суть метода PRT.

Однако, если BRDF зависит от направления v (зеркальный или анизотропный случай), то метод PRT значительно усложняется.

Далее в работе рассматривается (зеркальный) случай BRDF Фонга, отвечающей за блики на модели:

p(x,v) = CP (2а + 2)(rv)+ (6)

(зависимость от p понятна и опускается). Здесь r — зеркальный к x вектор, а > 0 — показатель Фонга, cp Е (0,1] — константа Фонга. Для этого случая применение метода PRT требует массивных вычислений, связанных с «подкручиванием» коэффициентов Фурье Llm для каждой вершины в шейдере (GPU) [5]. Становится сложным учет функции видимости. Поэтому ее либо не учитывают [1], либо аппроксимируют константой Ap = 2Vpo (метод Ambient Occlusion [6, 7]), Vp0 — нулевой коэффициент Фурье функции видимости. Для преодоления этих сложностей используется модификации метода PRT, например, основанные на дискретизации v [3]. Однако это приводит к необходимости использования больших объемов дополнительной памяти для хранения предрасчитанной информации.

2. Метод SDPRT

Автором совместно в Д.В. Горбачевым разработана модификация метода PRT — метод SDPRT (Spherical Design PRT). Этот метод позволяет обойтись без дополнительной памяти. В SDPRT реализуется аналог «подкрутки» коэффициентов с помощью простых вычислений. Он основывается на использовании взвешенных сферических дизайнов.

Взвешенным сферическим дизайном порядка s (s-дизайном) называются узлы xv и веса Xv кубатурные формулы

N

If = \ f (x)u(x) dp.(x) =Y^ Xv f (xv),

JS2 V=

точной для всех многочленов f (x) степени s [8, 9]; u(x) — неотрицательная весовая функция. Минимальным называют s-дизайн содержащий при заданном s наименьшее число узлов. Выделяют класс чебышевских дизайнов с равными весами Х\ = ... = Xn.

Дизайны из большого числа точек практически равномерно распределены на сфере и могут быть использованы в качестве узлов для метода Монте-Карло [10]. Например, в работе [4] показано, что результаты из [5], полученные методом Монте-Карло с 10000 точками, сгенерированными

по формуле ______

Pv = 2n£v, 0v = 2 arccos(^/T—ni), (8)

(£v,nv — равномерно распределенные на [0,1] точки, pv,dv — сферические углы точек xv), с той же точностью достигаются при использовании чебышевского дизайна из (всего) 48 точек.

Дизайны изучаются, в основном, для случая единичной весовой функции: и = 1 [9]. Однако нам принципиально потребуются дизайны небольшого порядка с произвольной весовой функции.

2.1. Построение дизайнов для зональной весовой функции.

Пусть и(х) = ио(ах) = ио(совв) — зональный вес, а € £2 — произвольная точка. Рассмотрим в этом случае известный метод построения дизайнов на основе одномерных квадратурных формул типа Гаусса-Маркова.

Запишем интеграл в (7) в сферических координатах:

1 г2п гп

I/ йр / (р,в)и0(еО8 в) вШ в(1в, (9)

Уо Уо

где / — многочлен степени в. Кубатурную формулу (7) достаточно проверить на мономах

/(х) = х^1 х22х33 = сов"1 р ■ вт^2 р ■ ‘$т1'1+1'2 в ■ сов"3 в, щ + и2 + и3 ^ в. (10) Имеем

1 г2п 1 гп

I/ = — сов1"1 р ■ вт^2 рйр ■ - в\п}'1+}'2 в ■ слов1'3 в и0 (сов в) вт в йв =

2п Уо 2 Уо

= II ■ 12 (11)

Используя формулу прямоугольников для тригонометрических многочленов степени VI + ^2, получим

1 г2п 1 5

II = ----- сов1"1 р ■ вШ^2 рйр = ---------- сов^1 рг ■ вШ^2 рг, (12)

2п ^ в+1

где

2пг

рг - ------

8 + 1

При этом І1 = 0, если VI + У2 — нечетное число. Отсюда следует, что интеграл І2 достаточно рассматривать, если + V2 — четное число.

В этом случае случае

+и2 в ■ совиз в = (вш2 ф^1+У2>/2 С08из в =

= (1 — сов2 в)(У1+У2)/2 соб”3 в = д(сов в), (13)

где д — тригонометрический косинус-многочлен порядка Р\ + V2 + из ^ ,в. Выполняя замену переменного ^ = соб в, получим:

1 Г 1 Г1

І2 = ^ д(сов в)ш0(сов в) біп вйв =- д(і)ш0(і) М, (14)

2 Уо 2 .1-1

где д(і) — алгебраический многочлен степени 8.

Применим к (14) квадратурные формулы Гаусса-Маркова [8] с фиксированными узлами:

г+1

І2 — ^2 7и д{і1ц), Ь,г+і — I,

3 = 1

г+2

І2 — ^ І2Ід{і2]), £2,г+1 — 1, £2,г+2 — — 1-

3 = 1

(15)

(16)

Здесь 2г ^ в и 2г + 1 ^ в соответственно, І1ц — нули ортогонального многочлена Q^1’0\і) с весом {1 — і)и0{і), і2] — нули ортогонального многочлена Qrl’1)(і) с весом {1 — і2)и0{і). При этом для неотрицательной весовой функции ио все веса будут положительными.

Для построения ортогональных многочленов на отрезке [—1,1] с весовой функцией р{і), можно воспользоваться методом ортогонализации системы степеней [11].

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

Л, —

в+1

7и;

для I — 1

для I — 2

хV — {рі, вц), І — 0, 1,...,в, І — 1,...,Г,

хм — {0,п), N — {в + 1)г + 1;

хV — {рі, в 2] ), І — 0, 1,...,в, І — 1,...,Г,

хм-1 — {0,п), хм — {0,п), N — {в + 1)г + 2

(17)

(18)

В результате получим дизайн порядка в из N точек с положительными весами. Если в = 2г — четное, то используем узлы (17) N = + 1), если

в = 2г + 1 — нечетное, то используем узлы (18) N = + 2).

2.2. Построение дизайнов небольшого порядка для произвольной весовой функции. Рассмотрим случай произвольной весовой функции и. Выписывая формулу (7) для сферических гармоник, получим систему нелинейных уравнений:

N

ит =/ ут (х)и(х) й^(х) = V ут (XV), т = = 0,...,в.

^2 5=1

1

Может оказаться, что для выбранного N система (19) не имеет решения. Поэтому будем минимизировать квадрат ошибки:

Е = ^ ^ АиУ\т(ху) - ит ^ тт . (20)

\ I >о,

1т 4 и ' XV еЯ2

Если оптимальное значение Е = 0, то соответствующее ему множество {XV, Хи}М=1 будет дизайном порядка в.

Преобразуем выражение для Е:

Е=

'У ' Аиу1т(хи) I 2 ^ ' Аиу1т(хи)и1т + и1т

= Б1 - 2Б2 + Б3.

(21)

Пусть

в I

и (х) = Е Е и1ту1т(х) ■ (22)

1=0 т=—I

Для многочленов степени в справедливо следующее представление через сферическое ядро Дирихле Кв:

и(х) = и(х')Кв(х,х') йц,(Х) = и(х')Кв(х,х') й^(х'), (23)

Зя2 Зя2

х \х\

Кв(х,х') = ^^(21 + 1)\х\1 Р\(х х'~\ , х Е М3, х Е Б2, 1=0 ^\х\ '

ЩЬ) — многочлены Лежандра [11].

По переменной х ядро Кв(х,Х) является сферическим многочленом степени в. Если х Е Б2, то [4]

Кв(х, х) = кв(ххГ) = ^2 ^2 У1т(х)У1т(х').

1=0 т=—1

Возьмем дизайн {хти, Ати}М=т1 порядка ^ 2в. Тогда

и(х) — ^ ' Атии(хти')Кв(х,хти). (24)

и=1

Для приближенного вычисления значений

и(хти)= / и(х')кв(хтиX) й»(х) (25)

я2

2

воспользуемся дизайном из большого числа точек {Аьи,хьи}М= 1:

• и=

N

и(хти) ~ ^2 АЬии(хЬи)кв(хтихЬи). (26)

V=1

Эти значения вычисляются один раз.

Нам потребуется градиент функции и(х). Он равен:

и (х) — ^ ' Атии(хти')Квх(х,хти)> (27)

и=1

где Квх(х, хг) — градиент ядра Кв(х, хг) по переменной х:

К3х(х,хг) =

= (21 + 1) 1\х\1-2хрА х х^ + \х\1-г(х'\х\2 — х(хх'))Р! ( х х'^

1=0 L \\х\ / \\х\ /

Если х Б2, то

Квх(х, х') = ^2,(21 + 1) [1хР1(хх') + (х' — х(хх'))Р' (хх')\ 1=0

в

= х£ (21 + 1)1Р1 (хх') + (х' — х(хх'))кв (хх').

— х

1=0

Имеем:

Бз = ^ и2т = \\и\\2 = I и2(х) й^(х) = ^ Атк™2(хтк) (28)

1,т ^я2 к=1

(это значение вычисляется один раз),

Б2 "У ' Аиу1т(хи)и1т У ' АуУ ' и1ту1т(хи) ^ ' Аии(хи^ (29)

1,т,и и 1,т V

Б1 — ^ ' I ^ ' Аиу1т(хи) I — ^ ^ ' АиА^у1т(хи')у1т(х^) —

1,т \ V ) 1,т и,^

У ' Аи А^У ' у1т(хи ')у1т(х^) = У ' Аи Ацкв(хихц). (30)

и,^ 1,т и,^

Таким образом, получаем, что

N N

Е =^2 Аи А^кв(х„ хц) — 2^2 Аи и(хи)+ Бз. (31)

и,^=1 и=1

Для поиска минимума функции Е воспользуемся методом линеаризации [12]. В качестве начальных значений весов берем числа иоо/^, где

, N

и00 — / u(x) d^(x) Xbvи(Xbv)

JS2 v=l

Начальные узлы выбираются случайный образом, например, по формуле (8) На каждой итерации полагаем

Av —— Av + aAAv, xv

xv + f3Axv \xv + l3Axv |

(32)

Линеаризуем Е, считая ДА, Ах малыми. Для этого подставим Аи + ААР хи + Ахи в (30) и отбросим квадратичные члены. Получаем

■ N / N \

Е ^ Е + 2 £АЧ£А ^к.3(хихц) - ш(хи)\ +

.и=1 \^=1 /

N / N '

+ У \ Аи Ахи {у ' А^Кзх(хи, х^) ^ (хи)

и=1 \р=1 /

(33)

Чтобы уменьшить значение E положим:

AXv — —

Axv — -

A^kg(xuxj) - u(xu) j ,

У ' AjKsx(,xv, xj) u (xv )^ •

(34)

(35)

Вычисления завершаются, когда текущее значение E будет удовлетворять заданной точности.

Для применения описанного метода необходимо для заданного s выбрать параметры N (число точек) и а, в (шаги по весам и узлам). Д.В. Горбачевым предложено брать эти параметры теми же самыми, что и для тестового случая и = 1. В частности, при нахождении взвешенного 2-дизайна по произвольной весовой функции можно взять N — 4 — значение минимального чебышевского дизайна для и = 1.

2.3. Расчет освещения при использовании Ambient Occlusion.

Рассмотрим случай е-низкочастотного окружающего освещения L (2) и BRDF Фонга (6).

Аппроксимируем функции видимости Vp в каждой точке модели p величиной Ap (Ambient Occlusion): Vp(x) & Ap. В наиболее важных случаях Ap ~ 0 (почти полное затенение) или 1 (почти полная открытость)

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

Вр(у) = Ь(х)Ур(х)р>р(х, у) й^(х) к> Бр(у) = Ар Ь(х)р>(х,у) й^(х), (36)

Зя2 Зя2

где р(х, у) = вР(2а + 2)(гу)+.

Аппроксимируем Ь частичной суммой Фурье Ь. Тогда аналогично (3)

имеем _

Вр(у) Вр(у) = Арвр Ь(х)ш0(ху) й^(х), (37)

зя2

где ш0^) = (2а + 2)(1)+ При этом для любого у в силу (4)

\Вр(у) - Вр(у)\ ^ еЦЬуАрвр||шо||2 < 5.

Здесь 0 ^ Ар ^ 1, 0 < вР ^ 1,

2 = (2а + 2)21 -()+йі = а+/2 ^ а + 2, а> 0-

Например, если 5 = 10~2, ||Ь||2 ^ 1, є = 10-3, то Vа + 2 ^ 10 и а ^ 102. Таким образом, метод работает при а < 102.

Запишем (37) в локальной системе координат:

Бр(у) = Арвр Ь(д„,ух)шо(хз) йц,(х). (38)

Зя2

Здесь дпо — вращение, переводящее у в п.

Построим дизайн порядка в {Ха^,хаис весовой функцией ш0(х3) по алгоритму, описанному в п. 2.1. Тогда

^ Г ~ Ма ~

Бр(у) = Арвр Цд^ох)шо(хз) йр,(х) = Арвр Е Ь(дпо хаи ). (39)

^Я'2 и=1

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

Формулу (39) можно применять при расчетах на ОРИ, если предварительно дискритизировать многочлен Ь, сохранив результат в текстуру достаточного объема, после чего в шейдере брать значения Ь(дпохи) из текстуры.

Предложим способ, как обойтись без дополнительной текстуры. Из равенства

1з(х) = [ /(x,)ks(xx,) й^(х') я2

имеем:

Ь(дха„)= / ks(gxаvх)Ь(х) йц,(х). (40)

я2

Здесь интеграл можно вычислить, воспользовавшись 8-дизайном [\ь^,хь^}^\ с весовой функцией Ь(х), построенным по методу, описанному в п. 2.2:

Мь

Ь(дх) = ^2 Хь^кз(дххь^). (41)

1л=1

Окончательно получаем следующую простую расчетную формулу для построения итогового изображения:

_ Ма Мь

В р(у) = Ар вр ЕЕ к3 (дхаи • хЬ^)- (42)

и=11-1=1

При в = 2 достаточно взять Ыа = = 4.

2.4. Вычисление Ар. Требуется рассчитать значения Ар для каждой вершины модели. Это сводится к вычислению интеграла

Ар = 2т I Ур(х) йх.

2п -)пх>0

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

Для ускорения расчетов предлагается использовать следующий подход [13]. В текущей вершине модели поместим 6 камер ориентированных по нормали п. Сделав снимок буфера глубины (или черно-белое изображение — черная модель на белом фоне) каждой камерой, получим СиЪеМар из текущей вершины модели. После этого поиск пересечения луча и модели сводится к считыванию информации из соответствующего пикселя СиЪеМар. Но для поиска среднего значения Ар требуется цикл по всем пикселям текстуры, что может быть затруднительно. Например, цикл по 6 • 128 • 128 = 98304 пикселей в каждой из 500 000 вершин модели. Поэтому на данном этапе предлагается считать значение Ар с помощью дизайнов. Эту операцию также можно выполнить быстро, распараллелив вычисления на ОРИ, т.е. воспользоваться такими технологиями как ОрепСЬ, КУГО1Л Quda, ЛМБ ИгеБ^еаш. Описанный подход позволяет в несколько раз ускорить вычисления.

3. Результаты

Как видно из табл. 1 для вычислений требуется меньше видеопамяти, меньше время для предрасчетов, а также вычисления, производимые в шейдере, значительно проще.

Результаты работы программы можно также увидеть на видеороликах: goo.gl/vEksM, goo.gl/ndEzK

Таблица 1

Сравнение результатов

SDPRT PRT [3]

Модель 130 тысяч 50 тысяч

Память глобальная: 64 float; в вершинах: 1 faloat в вершинах 25x25 float

Предрасчеты 35 минут 2.5 часа

Фраймрейт 3.6/16/125 300-400

Список литературы

1. Ramamoorthi R. A signal-processing framework for forward and inverse rendering // Ph. D. dissertation. Stanford University. 2002.

2. Akenine-Moller T, Haines E., Hoffman N. Real-Time Rendering. Third Edition, 2008.

3. 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.

4. Горбачев Д.В., Иванов В.И., Странковский С.А. Моделирование освещения в интерактивной графике при помощи сферических дизайнов // Изв. ТулГУ. Сер. Естественные науки. 2007. Т.1. Вып.1. С.17-36.

5. Green R. Spherical Harmonic Lighting: The Gritty Details // http://www.research. scea.com/gdc2003/spherical-harmonic-lighting.html.

6. Kontkahen J., Samuli L. Ambient occlusion fields // http://www.tml.tkk.fi/~janne/ aofields/aofields.pdf.

7. Shanmugam P., Arikan O. Hardware accelerated ambient occlusion techniques on GPUs // http://perumaal.googlepages.com/ao.pdf.

8. Мысовских И.П. Интерполяционные кубатурные формулы. М.: Наука, 1981.

9. 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.

10. Андреев Н.Н., Юдин В.А. Наименее уклоняющиеся от нуля многочлены и кубатурные формулы чебышевского типа // Труды МИРАН. 2001. Т.232. С.45-57.

11. Бейтман Г., Эрдей А. Высшие трансцендентные функции. Т.2. М.: Наука, 1974.

12. Амосов А.А., Дубинский Ю.А., Копченова Н.В. Вычислительные методы для инженеров. М.: Высшая школа, 1994.

13. http://http.developer.nvidia.com/GPUGems/gpugems_ch17.html.

Свистунов Сергей Сергеевич ([email protected]), аспирант, кафедра прикладной математики и информатики, Тульский государственный университет.

Interactive rendering using spherical designs (SDPRT) for the low-frequency lighting and Phong BRDF

S.S. Svistunov

Abstract. A new SDPRT method for real-time rendering for the low-frequrency lighting. It gives results similar to those obtained by the method of PRT. The method is based on the use of specially constructed spherical designs. This is achieved by simplifying the formulas for the PRT, which allows you to implement them effectively in the shader. The proposed method is less demanding of computational resources and requires fewer predraschitannoy information.

Keywords: real-time rendering, spherical designs, radiance, low-frequency light, precomputed radiance transfer (PRT).

Svistunov Sergey ([email protected]), postgraduate student, department of applied mathematics and computer science, Tula State University.

Поступила 05.02.2011