Нелинейная экстраполяция излучения Текст научной статьи по специальности «Математика»

Нелинейная экстраполяция излучения

К. А. Востряков Отдел компьютерной графики Институт Прикладной Математики М. В. Келдыша РАН vostryakov@gin.keldysh.ru

Аннотация

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

Уравнение визуализации

Расчет освещения на поверхности можно записать и решить методом Монте-Карло достаточно просто, но эффективность прямого решения низкая. Математически рассеивание света на поверхностях сцены можно записать уравнением визуализации [1]:

L(x, w) = Le (x, w) + J p(x, w, w')L(x, w') cos в dw'

W

Отраженное излучение L в точке x в направлении rn равно эмиссии излучения Le и интегралу по всем направлениям на полусфере Q по произведению двулучевой

функции отражения р (далее ДФО), падающего излучения Ь и косинуса угла между нормалью к поверхности и направлением а'. Существует и другая эквивалентная формулировка этого уравнения с интегралом не по телесному углу, а по площади поверхности всей сцены А:

с V (х, у)оо$,(пхху)оо$,(п уух) Ь(х а) = Ь(ха) +1 Р(У ® ха)Ь(У ® х)-г-^---

А II х - У ||2

где Ь(у ® х) — падающее излучение из точки у в точку х, р(у ® х,а) — ДФО в точке х, V(x, у) — функция видимости, которая принимает значение 1, если нет объектов между х и у, и 0 иначе, сов(ихху)— косинус угла между нормалью в точке х и направлением из х в у, соБ(иуух) — косинус угла между нормалью в точке у и

направлением из у в х. Ясно, что подынтегральное выражение не равно нулю только в точках на поверхностях источников света (первичных или вторичных), которые находятся в прямой видимости от точки х.

Когерентность

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

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

Алгоритмы расчета освещения от протяженных источников света

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

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

В 1997 году Keller представил новый алгоритм для расчета глобального освещения, который был назван алгоритмом быстрой излучательности [2]. Генерируется небольшое количество виртуальных точечных источников света, используя метод квази-случайного блуждания. Затем происходит визуализация сцены

с использованием затенения от этих виртуальных источников, как от обычных. Большинство вычислений может быть выполнено на графическом процессоре. Этот алгоритм использует дискретизацию только вторичных источников света, поэтому он не имеет артефактов связанных с дискретизацией геометрии, но создает свои артефакты из-за дискретизации виртуальных источников, что проявляется в виде множества жестких теней, где должна быть мягкая тень. Позднее алгоритм был усовершенствован [3] и назван алгоритмом быстрого глобального освещения, где затенение выполнялось выполнялось трассировкой лучей с перемежающейся выборкой виртуальных источников света. Быстрая излучательность и быстрое глобальное освещение плохо подходит для неламбертовских ДФО и сцен со множеством преград.

Классический алгоритм теневых объемов для теней с жесткой границей был обобщен для мягких теней [4]. Поскольку алгоритм основывается на поиске силуэта объекта, который отбрасывает тень, то скорость мягких теневых объемов зависит от количества ребер создающих тень. Поэтому данный метод плохо подходит для сложных сцен с большими протяженными источниками света.

Методы, которые основаны на теневых картах [5], лучше масштабируются под размер и сложность сцен в сравнении с теневыми объемами. Сцена растеризуется из источника света, а в буфер записывается глубина пикселей. Далее, при растеризации из камеры, с помощью z теста можно определить: лежит ли точка поверхности в тени от данного источника света или нет. Немаловажную роль играет аппаратная поддержка растеризации, которая есть в современных графических процессорах. Ранние работы по мягким теневым картам просто усредняли несколько жестких теней, что естественно сильно ограничивало возможный размер «правильно» выглядящей полутени. Поздее стали вычислять затенение как долю заграждения геометрии сцены, спроецированной обратно на источник света. Методы в большинстве случаев пренебрегают физической корректностью, поэтому их трудно применять для больших протяженных источников света.

Часто для расчета освещения от протяженных источников света используют различные схемы предвычисления освещения на поверхностях сцены.

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

которая согласуется с градиентом освещения. Использование неламбертовских ДФО или невозможно, или очень ограничено. Главным достоинством метода является возможность визуализации карты освещения на графическом процессоре с использованием вершинного освещения.

Фотонные карты [6] сохраняют освещение как облако точек столкновения частиц с поверхностями сцены, которые были выпущенны из источников. Фотонные карты гораздо более гибкий метод, чем карты освещения и получили очень широкое распространение на практике. Однако, оценка освещения, полученная с помощью фотонной карты, имеет низкочастотный шум и для его устранения используется финальный сбор. Это трассировка лучей по полусфере согласно уравнению визуализации, где подынтегральная оценка освещения уже берется с помощью фотонной карты. Финальный сбор занимает подавляющую часть времени расчета.

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

Кэш облучения впервые был предложен в работе [8], что позволило в несколько раз ускорить визуализацию, поскольку только в небольшой части точек, видимых через пиксели, облучение вычислялось трассировкой лучей по полусфере, а в остальных точках оно экстраполировалось. Для того, чтобы вычислить значение излучения в некоторой точке, ищем в кэше уже вычисленные подходящие значения для экстраполяции. Если экстраполяция возможна, то выполняем ее и переходим к следующей точке. Если нет, то вычисляем освещение, трассируя лучи по всей полусфере, и записываем вычисленное значение в кэш. Для оценки погрешности экстраполяции первоначально [8] было использовано среднее гармоническое расстояние (далее СГР) до ближайших объектов сцены, которое вычислялось при трассировке лучей по полусфере. Чем меньше это расстояние, тем больше ошибка и следовательно больше плотность точек. Плотность точек в кэше облучения адаптируется под градиент освещения, тем не менее, трудно добиться отсутствия артефактов, поскольку градиент освещения вычисляется достаточно приблизительно.

Krivanek [10] предложил использовать многопроходный адаптивный алгоритм, который, при добавлении новой точки в кэш, сравнивает СГР и относительное отклонение облучения в соседних точках, и уменьшает радиус влияния новой точки, а также радиусы влияния точек в кэше. Понятно, что после этого часть точек, которые были успешно экстраполированы раньше, оказываются экстраполированы некорректно. Поэтому необходимо несколько проходов, пока все точки не окажутся вычисленными «корректно». Многопроходность существенно уменьшает количество артефактов изображения, но приводит к множеству чрезмерных вычислений, что снижает его эффективность.

В работе [9] метод кэша облучения [8] был усовершенствован добавлением градиентов облучения. Вовлечение градиентов позволяет осуществить линейную экстраполяцию или кубическую интерполяцию освещения.

Кэш излучения

В работе [10] было предложено расширить метод Ward'а, чтобы использовать любые ДФО, а не только Ламбертовские. Это потребовало хранение освещения как функции от направления на полусфере. Для этого Krivanek использовал сферические гармоники. Кроме сжатия данных, этот базис позволил быстро осуществлять свертку. Сферические гармоники были выбраны, поскольку проекция на базис выполняется быстро, а свертка представляет собой скалярное произведение векторов коэффициентов. Это позволило применять не только ДФО Ламберта, но и умеренно направленно-диффузные ДФО, однако метод не применим к остронаправленным ДФО.

Ближнее и дальнее освещение

В работе [11] освещение было разделено на ближнее и дальнее. Для дальнего освещения использовался метод кэша излучения [10] с базисом из сферических гармоник. Поскольку дальнее освещение изменяется достаточно плавно, то кэш излучения гораздо более редкий, чем для полного освещения. Для ближнего освещения применялся метод излучательности без учета функции видимости для расчета форм-факторов, поэтому ближнее освещение вычисляется очень быстро — достаточно найти все треугольники в окрестности точки. Хотя игнорирование видимости дает физически некорректное решение, визуально в большинстве случаев оно кажется «правильным». Для сцен средней сложности с диффузными поверхностями алгоритм подходит хорошо, создавая изображения по визуальному качеству сравнимые с [9], но значительно быстрее.

Репроекция

В курсе [6] описывается использование репроекции вместе с кэшем излучения. Точки пересечения лучей и поверхностей сцены, которые были вычисленны ранее и записаны в кэш, проецируются на полусферу в новой точке поверхности, в которой необходимо вычислить освещение. Как сообщает Toskiaki Kato, в визуализаторе «Kilauea» применение репроекции часто не давало ускорения. Репроеция — не быстрая процедура, поэтому необходима нормировка вектора и вычисление обратных тригонометрических функций (если используется сферические координаты для параметризации полусферы). Кроме того, необходимо использовать z-буфер и что-то делать с пустотами, которые останутся после репроекции. В данной работе мы попытаемся устранить недостатки алгоритма репроеции применительно для кэша излучения.

Выбор алгоритма

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

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

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

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

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

Алгоритм верхнего уровня

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

Для того чтобы хранить излучение, падающее с полусферы, в линейной памяти нужно иметь отображения с полусферы на плоскость и обратно. Было выбрано отображение [12], поскольку оно непрерывно, сохраняет площадь и имеет небольшое искажение. Если отображение из квадрата на полусферу будет с плотностью cosd/ж, тогда интеграл по полусфере превращается в интеграл по квадрату с координатами u, v е [0,1] х [0,1].

Создание записи кэша

На вход подается точка на поверхности, где было вычислено освещение, и нормаль к поверхности. На выходе имеем текстуру падающего излучения и массив п-элементов. П-элемент — это сокращение от поверхностный элемент (аналог термина surfel в английском языке, сокращение от surface element). Под этим термином понимается то, что в точке на поверхности имеется некоторая плоская дифференциальная площадка, параллельная поверхности в данной точке.

Вначале трассируем n х n лучей по полусфере (типичное значение 32 х 32, 64 х 64), находим точки столкновения с поверхностями сцены (или бесконечно удаленным фоном). Вычисляем падающее излучение, которое идет вдоль лучей, сохраняем в текстуру. Вычисленные точки разделяем на слои по глубине. Каждая точка становится п-элементом. Сортировка выполняется по K слоям, внутри слоя сортировать п-элементы не обязательно, поэтому сортировка может быть выполнена за K ■ O(n х n) времени, что достаточно быстро в сравнении с трассировкой n х n лучей. Слои выбирают по логарифмической шкале, поскольку ближние п-элементы нужно

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

Репроекция

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

Репроекция осуществляется слой за слоем, начиная от ближнего слоя и заканчивая самым дальним. Это похоже на обратный алгоритм художника. При этом полноценного Z-буфера не нужно, можно обойтись лишь маской заполненности. Проецируем все п-элементы в новой точке. Для каждого из них имеем в пространстве четыре угловые точки. В общем случаи они образуют некоторый четырехугольник (возможно даже не выпуклый). При этом по четырем точкам строим прямоугольник, координаты угловых точек которого — это максимальные uv координаты четырех репроецированных точек. Этим мы вычисляем приближенную проекцию п-элемента. Для п-элемента создаются маски п-элемента и маски неизвестности. Маска п-элемента — это битовая маска, где 1 расположены в тех местах ыу пространства, куда п-элемент был спроецирован. Далее следует создание битовой маски следа движения проекции п-элемента. Эта маска вычисляется как объемлющий прямоугольник координат п-элемента в точке из кэша и новой проекции п-элемента. Маска следа позволяет определить область ^ пространства, где открылась неизвестная зона. Маска неизвестности получается как маска следа минус маска заполнения. Маска неизвестности запрещает заполнение более дальних слоев в области неизвестности. Объект, который может находится в такой зоне, может перекрыть фон, видимый из новой точки. Маской неизвестности мы препятствуем тому, что перекрытый фон будет учтен неправильно. Следует отметить, что маска

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

^ « 7 2

пикселей маски деленное на число пикселей к .

Реализация и заключение

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

Алгоритм может быть использован для расчета финального сбора в методах фотонных карт и излучательности, освещения от карты окружения, заданной с большим динамическим диапазоном, теней от больших площадных источников света, «размытых» отражений. Алгоритм физически корректен и может использоваться для визуализации очень сложных сцен с высокочастотных ДФО и освещения, чего аналогичные алгоритмы кэширования обеспечить не могут. Алгоритм дает ускорение 2-10 раз в сравнении с классической реализацией метода Монте-Карло.

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

1. James T. Kajiya. The rendering equation. In SIGGRAPH '86: Proceedings of the 13th annual conference on Computer graphics and interactive techniques, pages 143-150, New York, NY, USA, 1986. ACM.

2. Alexander Keller. Instant radiosity. In SIGGRAPH '97: Proceedings of the 24th annual conference on Computer graphics and interactive techniques, pages 49-56, New York, NY, USA, 1997. ACM Press/Addison-Wesley Publishing Co.

3. Ingo Wald, Carsten Benthin, and Philipp Slusallek. Interactive global illumination using fast ray tracing. In Rendering Techniques 2002 (Proceedings of the Thirteenth Eurographics Workshop on Rendering), June 2002.

4. Ulf Assarsson and Tomas Akenine-Möller. A geometry-based soft shadow volume algorithm using graphics hardware. ACM Trans. Graph., 22(3):511-520, 2003.

5. Lance Williams. Casting curved shadows on curved surfaces. Computer Graphics, 12(3):270-274, August 1978.

6. Per H. Christensen, Per H. Henrik Wann Jensen, Toskiaki Kato, and Frank Suykens. A practical guide to global illumination using photon mapping. In SIGGRAPH 2002 Course Notes. Association for Computing Machinery, ACM SIGGRAPH, August 2002. Course 43.

7. P.P. Sloan, J. Kautz, and J. Snyder. Precomputed radiance transfer for real-time rendering in dynamic, low-frequency lighting environments. ACM Transaction on Graphics, 21(3):527-536, 2002.

8. Gregory J. Ward, Francis M. Rubinstein, and Robert D. Clear. A ray tracing solution for diffuse interreflection. In ACM SIGGRAPH '88, pages 85-92, New York, NY, USA, 1988. ACM.

9. Gregory J. Ward and Paul Heckbert. Irradiance gradients. pages 85-98, May 1992.

10. Jaroslav Krivânek. Radiance Caching for Global Illumination Computation on Glossy Surfaces. Ph.d. thesis, Université de Rennes 1 and Czech Technical University in Prague, December 2005.

11. Okan Arikan, David Forsyth, and James F. O'Brien. Fast and detailed approximate global illumination by irradiance decomposition. In ACM SIGGRAPH 2005, August 2005.

12. Peter Shirley and Kenneth Chiu. A low distortion map between disk and square. Journal of graphics tools, 2(3):45-52, 1997.