КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ

И МОДЕЛИРОВАНИЕ 2012 Т. 4 № 1 С. 7584

ЧИСЛЕННЫЕ МЕТОДЫ И ОСНОВЫ ИХ РЕАЛИЗАЦИИ

УДК: 519.2, 519.6

Локальные оценки метода Монте-Карло в решении

уравнения глобального освещения с учетом

спектрального представления объектов

В. П. Будакa, В. С. Желтовb, Т. К. Калакуцкийc

Национальный исследовательский университет «МЭИ»,

Россия, 111250, Москва, Красноказарменная улица, дом 14 a [email protected], b [email protected], c [email protected] E-mail:

Получено 8 ноября 2011 г., после доработки 20 ноября 2011 г.

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

Ключевые слова: глобальное освещение, Монте-Карло, локальная оценка, двойная локальная оценка Local estimations of Monte Carlo method with the object spectral representation in the solution of global illumination V. P. Budak, V. S. Zheltov, T. K. Kalakutsky National research university “MPEI”, Russia, 111250, Moscow, Krasnokazarmennay st., Abstract. – The article deals with the local and double local estimation of the Monte Carlo method for solving the equation of global illumination. The local estimation allows calculating the illumination at any point at the approximation of diffuse reflection, whereas the double local estimation allows calculating directly the luminance at a given point in a given direction. The article presents the mathematical basis of local estimations and the basic stages of the software implementation. The representation of three-dimensional objects in the basis of spherical functions and the possibility of using them in the local estimations are also considered.

Keywords: global illumination, Monte Carlo, local estimation, double local estimation Citation: Computer Research and Modeling, 2012, vol. 4, no. 1, pp. 75–84 (Russian).

© 2011 Владимир Павлович Будак, Виктор Сергеевич Желтов, Тимофей Кириллович Калакуцкий 76 В. П. Будак, В. С. Желтов, Т. К. Калакуцкий Введение Визуализация трехмерных сцен производится на основе решения уравнения глобального освещения, которое представляет собой интегральное уравнение второго рода [Будак В. П., 2000]:

L(r, ) = L0 (r, ) + L(r, ) (r;, ) (N, ) d, ll l l l l l (1) где L(r, ) – яркость в точке r в направлении ; (r;, ) – двунаправленная функция рассеяния l l ll (отражения или преломления); L0 – яркость прямого излучения непосредственно от источников света; N – нормаль в точке r к поверхности.

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

Трассировка лучей получила очень широкое распространение и включена в такие известные программы моделирования, как 3D Studio Max и Maya. При этом различают прямую и обратную трассировки, где, соответственно, трассируются лучи либо от источников света, либо от приемников (камеры).

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

К решению уравнения глобального освещения также применялся метод конечных элементов, получивший свое название в теории глобального освещения - метод излучательности (radiosity). Метод основан на допущении, что все элементы сцены имеют диффузное отражение, тогда уравнение (1) может быть записано в виде M (r ) F (r, r )(r, r )d 2r, M (r ) = M 0 (r ) + (2) где M(r) – светимость точки поверхности r; M0(r) – светимость точки r, полученная непосредственно от прямых источников света; (r, r) – функция видимости элемента d 2r из точки r ;

( N(r ),(r r )) (N(r ),(r r )) F= – элементарный форм-фактор, где N(r ) – нормаль в точке r (r r ) к поверхности.

В нашей работе мы предлагаем использовать локальные оценки метода Монте-Карло, хорошо известные при решении уравнения переноса излучения в атмосферной оптике [Марчук Г. И., 1980; Ермаков С. М., Михайлов Г. А., 1976]. Также в работе рассматривается возможность визуализации трехмерных объектов, заданных в спектре сферических гармоник.

Локальная оценка Уравнение (1) не является удобным для статистического моделирования, так как искомая функция под интегралом стоит в точке r, а определяется в точке r. При этом переменные r и l не являются независимыми, а связаны соотношением = r r.

l (3) r r КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Локальные оценки метода Монте-Карло в решении уравнения глобального освещения Соответственно мы можем записать уравнение (1) в виде Ядро уравнения глобального освещения (4) содержит в себе -функцию, что определяет особенность в его угловом распределении и делает фактически невозможным моделирование яркости методами Монте-Карло. Устранить особенность можно проинтегрировав по пространству, что при замене отражения на диффузное будет эквивалентно уравнению (2). В результате оценка I для L(r, ) примет вид где k(r,r) – ядро уравнения (4); Qn – вес марковской цепи; M – математическое ожидание различных случайных траекторий лучей.

Марковская цепь моделирует последовательный ряд случайных лучей, блуждающих по сцене. Вес первого луча Q1 определяется начальной яркостью источника L0 (r, ). При этом k(r,r) определяет вероятность перехода из точки марковской цепи в заданную точку r. [Марчук Г. И., 1980; Ермаков С. М., Михайлов Г. А., 1976]. Со всех узлов траектории, где произошло ее пересечение с поверхностью сцены, рассчитывается вклад в освещение сцены в точке r на основе ядра уравнения k(r,r) в выражении (5). При этом все распределения вероятностей должны отвечать условиям нормировки.

Величины L0 (r, ) и k(r,r) будут определять статистические веса марковской цепи Qn.

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

Выражение (5) получило название локальная оценка метода Монте-Карло [Марчук Г. И., 1980; Ермаков С. М., Михайлов Г. А., 1976]. Она позволяет проводить оценку освещенности в заданной точке сцены. Таким образом, для расчета освещенности в некоторой заданной точке r необходимо построить марковскую цепь и на каждом акте отражения вычислять ядро k(r,r) для всех заданных точек. Математическое ожидание полученной величины и будет равно освещенности. На рисунке 1 представлена общая схема построения марковской цепи и расчетов освещенности методом локальной оценки в заданной точке.

Марковская цепь – это последовательность случайных событий с конечным или счётным числом исходов, характеризующаяся тем свойством, что, говоря нестрого, при фиксированном настоящем будущее не зависимо от прошлого. В данном случае цепь строится от источника света, при этом начальный вес кладется равным силе света источника света с учетом нормировки на поток. После чего находится точка пересечения луча с элементом сцены и вес умножается на коэффициент отражения. Далее вычисляется ядро уравнения (4) для каждой из исследуемых точек, умноженное на вес луча, и суммируется с предыдущими значениями. После чего производится розыгрыш луча в соответствии с диффузным законом отражения и ищется следующее пересечение. Процесс повторяется до тех пор, пока луч либо не покинет пределы сцены либо его вес не станет ниже заданного порога. Разыграв большое количество лучей и усреднив, мы получим значения освещенности в исследуемых точках. Особо отметим, что локальная оценка позволяет по одному лучу вычислять значение освещенности сразу в нескольких точках. Это принципиальное отличие от прямого моделирования и трассировки лучей.

Точность локальной оценки Идеальным вариантом для сравнения любого численного метода является наличие точных аналитических решений частных случаев. Для уравнения (2) существуют два аналитичеТ. 4, № 1, С. 75–84 _ ских решения. Первое из них – фотометрическая сфера. Однако для сравнения точности она слабо подходит из-за полной симметрии. Вторым частным случаем являются две бесконечные параллельные плоскости и точечный источник между ними, получившим название задача Соболева [Соболев В.В., 1944]. Полученное в статье Соболева решение не является удобным для непосредственных расчетов. Рассмотрим решение задачи Соболева, позволяющее получить результат в более приемлемой аналитической форме.

Рис.1. Схема построения марковской цепи и расчета локальной оценки При рассмотрении задачи Соболева, уравнение (2) преобразуется в систему из двух интегральных уравнений. Каждое из уравнений описывает распределение освещенности по одной из плоскостей. Соответствующее уравнение для первой плоскости примет вид (для второй будет аналогичное выражение с другими индексами):

где Ei(r) – освещенность i-ой плоскости (i=1,2); i – коэффициент i-ой плоскости; hi – расстояние от i-ой плоскости до источника. Мы положили силу света источника равную 1 и h1+h2=1.

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

где K1(k) – модифицированная функция Бесселя чисто мнимого аргумента или функция МакДональда первого порядка; J0 – функции Бесселя нулевого порядка.

Выражение (7) является удобным для расчетов на компьютере и позволяет проводить сравнение математических методов моделирования. На рисунке 2 представлено сравнение освещенности полученной нами локальной оценкой метода Монте-Карло и точным решением задачи Соболева (7).

В расчете локальной оценкой мы использовали 2000 лучей, при этом время расчета составило менее 1 секунды на компьютере с процессором AMD Athlon 64 X2 5200.

КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Рис. 2. Распределение освещенности в задаче Соболева: h1 = h2 = 0.5 и коэффициенты отражения 1 = 2 = 0. Также в рамках нашей работы было реализовано моделирование уравнения (2) с помощью метода излучательности. В результате при одинаковой заданной точности локальная оценка сходится быстрее в среднем в 80–90 раз на примере задачи Соболева.

Двойная локальная оценка Уравнение глобального освещения может быть представлено в операторной форме Решение этого уравнения представимо в виде ряда Неймана, что позволяет провести следующие преобразования в нормальной форме это примет вид Локальная оценка, соответствующая выражению (10), может быть названа двойной локальной оценкой [Марчук Г. И., 1980; Ермаков С. М., Михайлов Г. А., 1976] и будет иметь вид где В выражении (12) угловая особенность пропадает в результате интегрирования, а независимые переменные r,, r,, r, соответствуют геометрии распространения луча [Ермаков С. М., Михайлов Г. А., 1976].

Таким образом, двойная локальная оценка позволяет моделировать непосредственно уравнение глобального освещения (4) и вычислять яркость в заданной точке в заданном наТ. 4, № 1, С. 75–84 _ правлении от кратностей выше первой. Первая кратность содержится в члене KL0 уравнения (9) и может быть вычислена непосредственно.

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

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

Одним из перспективных направлений в представлении трехмерных объектов стало использование разложения объектов в базисе сферических функций [Mousa M., Chaine R., Akkouche S., 2006]. Подобное описание объектов предоставляет возможность регулирования качества воспроизведения объектов. Так, объекты, находящиеся на большом удалении от камеры при визуализации могут воспроизводиться с низким качеством. С точки зрения фотометрии, существенным преимуществом данного подхода является непрерывное воспроизведение нормалей без аппроксимации.

Рассмотрим математический аппарат, лежащий в основе сферических гармоник.

Сферическая гармоника {Ykm (, ) : m k } – это специальная функция, определенная на единичной сфере где – зенитный угол [0 ]; – азимутальный угол [0 2]; Akm, Bkm – некоторые коэффициенты;

Pkm (cos ) – присоединенные полиномы Лежандра. Основное свойство сферических гармоник – их ортогональность:

где – дельта функция Кронекера.

Система сферических функций является полной, и любая дважды непрерывно дифференцируемая функция, заданная на сфере, может быть разложена в ряд в спектре сферических гармоник [Тихонов А. Н., Самарский А. А., 1976]:

где Akm, Bkm – коэффициенты Фурье, определяемые как КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ и Ykm – норма, определяемая выражением Особую сложность при разложении составляет вычисление факториалов и, как следствие, существенное отличие коэффициентов по порядку величины, приводящее к дополнительной погрешности вычислений. При разложениях по сферическим функциям более удобно использовать полиномы Шмидта, определяемые как Нетрудно видеть, что при разложении по ним отпадает необходимость вычислять факториалы в норме, что существенно повышает производительность при вычислениях. Таким образом, любой объект, однозначно определенный по радиусу относительно некоторой точки, может быть разложен по сферическим гармоникам.

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

Нормаль к функции (r0, 0,0 ), заданной в сферических координатах, равна градиенту в этой точке:

Для поверхности, заданной в сферической системе координат, можем записать Найдем частные производные в соответствии с уравнением (19). Производная по r находится элементарно Найти производную по также не представляет сложности:

В свою очередь, производная потребует взятия производной по полиномам Шмидта:

Известно рекуррентное соотношение для полиномов Лежандра [Виленкин Н. Я., 1965] _ 2012 Т. 4, № 1, С. 75–84 _ с учетом соотношения между полиномами Шмидта и Лежандра [Виленкин Н. Я., 1965] проведя ряд преобразований, окончательно мы можем получить выражение для производной полиномов Шмидта Отдельно необходимо рассмотреть случай m=0. В этом случае с учетом ряда известных соотношений [Mishchenko M.I., Travis L.D., Lacis A.A., 2002] можно получить выражение Для дальнейшей работы необходимо записать нормаль в декартовой системе координат (e r, e, e ), не связанной с положением точки. С учетом ряда известных соотношений можно получить окончательное выражение для нормали к объекту заданному сферическими:

Поиск пересечения с объектом, заданным сферическими гармониками, также не тривиальная задача. Рассмотрим объект, заданный сферическими гармониками и расположенный в центре декартовой системы координат, и луч, заданный в этой же системе в векторной форме Косинус угла точки пересечения поверхности и луча относительно центра объекта может быть выражен как Вектор соответственно Косинус и синус угла примут вид В результате мы получили зависимость углов и от одной переменной. В точке пересечения имеет место равенство КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ Локальные оценки метода Монте-Карло в решении уравнения глобального освещения Уравнение (33) содержит множество решений, в данном случае нас интересует только первая точка пересечения. Она может быть найдена с помощью различных методов, начиная от низкоэффективных, таких как последовательное приближение, так и более сложными алгоритмами, например, генетическими. Применение конкретного алгоритма зависит от требований производительности и точности. В нашей реализации мы последовательно с определенным шагом вычисляли значения и локализовывали положение первой точки пересечения, после чего методом вилки находили его с заданной точностью.

Реализация локальной оценки с учетом спектральных объектов Мы реализовали алгоритмы локальной и двойной локальной оценки с возможностью представления трехмерных объектов в базисе сферических функций в среде Matlab. На рисунке 3 представлена визуализация простейшего помещения с одним изотропным источником света и головой человека, представленной сферическими гармониками. Изначально голова была представлена сеткой из 32 654 граней. В сцене голова была преобразована в 32 сферические гармоники. Мы провели сравнение времени расчета в случае представления объекта сферическими гармониками и стандартной сеткой. Время расчета при использовании сферических гармоник выросло на 40%, однако при этом объем хранимой информации снизился более чем в 1000 раз.

Рис. 3. Трехмерная сцена, визуализированная с помощью локальной оценки Заключение Локальные оценки метода Монте-Карло позволяют проводить физически адекватное моделирование уравнения глобального освещения. При этом, в отличие от методов прямого моделирования, они позволяют проводить оценку по одному лучу сразу во всех исследуемых точках. Это позволяет значительно повысить эффективность расчетов.

Локальная оценка позволяет вычислять освещенность в заданных точках трехмерной сцены. Благодаря своей эффективности, существенно превосходящей используемый в настоящее время в светотехнических расчетах метод конечных элементов (программные продукты DIALux, Relux), этот метод становится новым вектором в развитии светотехнических расчетов.

Двойная локальная оценка впервые в светотехнической практике позволяет вычислять непосредственно яркость в заданной точке по заданному направлению. Это открывает принципиально новые горизонты светотехнического проектирования и возможность пересмотра существующей нормативной базы, основанной на существующих методах моделирования. Двойная _ 2012 Т. 4, № 1, С. 75–84 _ локальная оценка может стать новым численным методом, используемым при реалистичной визуализации трехмерных сцен, с учетом того, что добиться реалистичной визуализации можно только при физически адекватном моделировании уравнения глобального освещения.

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

Список литературы Будак В. П. Визуализация распределения яркости в трехмерных сценах наблюдения. – М.:

Виленкин Н. Я. Специальные функции и теория представлений групп. – М.: Наука, 1965.

Ермаков С. М., Михайлов Г. А. Статистическое моделирование. – М.: ФИЗМАТЛИТ, 1982.

Метод Монте-Карло в атмосферной оптике. / Под ред. Г. И. Марчука. – Новосибирск: Наука Соболев В. В. Точечный источник между параллельными плоскостями // ДАН СССР, 1944.

Т. 42, № 4. С.176–177.

Тихонов А. Н., Самарский А. А. Уравнения математической физики: Учебник. – М.: МГУ, 2004.

Mishchenko M. I., Travis L. D., Lacis A. A. Scattering, Absorption, and Emission of Light by Small Particles. – Cambridge: Cambridge University Press, 2002.

Mousa M., Chaine R., Akkouche S. Frequency-Based Representation of 3D Models using Spherical Harmonics // The 14-th International Conference in Central Europe on Computer Graphics, Visualization and Computer Vision "WSCG'2006", Plzen, Czech Republic, January 31 – February 2, 2006. Full Papers Proceedings, p. 193–200.

КОМПЬЮТЕРНЫЕ ИССЛЕДОВАНИЯ И МОДЕЛИРОВАНИЕ