Локальная модель взаимодействия света с изотропными и одноосными прозрачными средами Текст научной статьи по специальности «Математика»

УДК 004.925.3

В. А. Дебелов , Д. С. Козлов 2

1 Институт вычислительной математики и математической геофизики СО РАН пр. Акад. Лаврентьева, 6, Новосибирск, 630090, Россия

2 Новосибирский государственный университет ул. Пирогова, 2, Новосибирск, 630090, Россия

E-mail: debelov@oapmg.sscc.ru; kozlov@oapmg.sscc.ru

ЛОКАЛЬНАЯ МОДЕЛЬ ВЗАИМОДЕЙСТВИЯ СВЕТА С ИЗОТРОПНЫМИ И ОДНООСНЫМИ ПРОЗРАЧНЫМИ СРЕДАМИ *

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

Ключевые слова: изотропный кристалл, одноосный кристалл, оптическая ось, поляризация света, оптическая дисперсия, двойное лучепреломление, фотореалистическая визуализация, лучевая трассировка.

Введение

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

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

Одноосные среды являются оптически анизотропными, т. е. их оптические характеристики зависят от направления измерений. Законы преломления и отражения света в таких средах не всегда подчиняются классическому закону преломления Снеллиуса для изотропных сред. Так, например, в одноосных кристаллах наблюдаются эффекты двойного лучепреломления

* Работа выполнена при поддержке РФФИ по гранту № 09-07-00237-а.

Авторы благодарят сотрудника кафедры минералогии и петрографии ГГФ НГУ С. З. Смирнова за постоянное внимание к нашим разработкам, а также Ю. А. Данькина за квалифицированную огранку кристалла кальцита.

ISSN 1818-7900. Вестник НГУ. Серия: Информационные технологии. 2012. Том 10, выпуск 1 © В. А. Дебелов, Д. С. Козлов, 2012

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

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

Все кристаллы могут быть разделены на изотропные, одноосные и двуосные по числу оптических осей: 0, 1 или 2 соответственно. В природных условиях кристаллы могут быть подвержены действию различных сил (давление, электромагнитные поля и т. д.). В таких условиях одноосный кристалл может вести себя как оптически двуосный. Более того, даже вакуум в надлежащих условиях (эффект Керра) ведет себя как оптически анизотропная среда, в частности, наблюдается двойное лучепреломление. На практике такие влияния сторонних сил несущественны и могут быть исключены из рассмотрения. В данной работе мы рассматриваем только идеальные оптически неактивные немагнитные кристаллы в условиях отсутствия влияния существенных внешних сил.

В работах [2; 3] рассматривается взаимодействие света с границей между изотропной и одноосной средой, а случай взаимодействия луча с границей между двумя одноосными средами не рассматривается, и потому визуализация одноосных кристаллических агрегатов невозможна. В природе же наиболее часто встречаются именно кристаллические агрегаты, а не монокристаллы. Авторы книг по оптике кристаллов обычно не заботятся о приложениях компьютерной графики, поэтому они не содержат формул, которые могут быть легко интегрированы в алгоритмы визуализации, такие формулы должны быть специально разработаны, чему и посвящается данная работа.

В [3] используется тригонометрический подход при расчете лучей, в других работах, например, в [4], для вычислений используется матричный подход. Вычислительные формулы, изложенные в данной работе, основываются на ковариантном методе [5], применяя который, мы сконструировали формулы и алгоритм расчета взаимодействия луча с границей между изотропными и одноосными средами (во всех возможных комбинациях).

В данной работе мы рассматриваем только плоские электромагнитные волны [5; 6], являющиеся классическим решением уравнений Максвелла. Обращение только к такому типу волн является приемлемым и наиболее часто рассматриваемым случаем в приложениях компьютерной графики. В [3] для описания состояния поляризации луча используются векторы Джонса и матрицы Мюллера, мы же используем подход, основанный на матрицах когерентности и их модификаторах, предложенный в [7].

Краткое введение в тензорное исчисление

Трехмерным вектором будем называть вектор-столбец:

= {иа),

I и3 у где а е {1,2,3}.

Тензор второго ранга может быть представлен в виде матрицы

( и >

Ч1 а12 а > 13

а = а 21 а22 а23 = а аЬ ,

ча31 а32 а33 у

где а, Ь е {1,2,3}.

Если какие-либо величины (тензоры или векторы) стоят рядом без каких-либо значков между ними, то подразумевается их свертка, т. е. суммирование по соседним индексам. Например, ыу представляет скалярное произведение двух векторов: ыу = иаУа

и1У1 + и2 У2 + и3 V

■2У2

3 3 •

Свертка матрицы а c вектором и определяется как произведение матрицы на вектор у = аи.

Отметим, что величина иау = иаааЪУь является скаляром.

Для любых двух векторов определяется классическое векторное произведение:

[ыу ] = и

х у

и3 V - иъ3

V и1У2 - и2 V У

0 —и3

V и2

0

и1

2

-и1

и у,

а тензор и называется дуальным к вектору и.

Диадой двух векторов называется специальный тензор:

ии ЩУ2 ии

и •у=(и/ XV ) = и2 V и2 V и2 V

V и3^1 и3 V ЩУ3 у

Скаляр к эквивалентен следующей матрице:

к =

к00 0 к 0

V

0 0 к

У

Взаимной матрицей а к матрице а является матрица

( Л Л12 Л13

а = Л21 Л22 Л23

V Л31 Л32 Л33 33

где Ли- - алгебраическое дополнение элемента а]1 матрицы а.

.—. 1

сип ч^^^симт^^гчА^ч^ д^и^лп^ппс ^лсш^ша ]

Если ёе! а = |а| ф 0 , то обратная матрица равна а-1 = -

-а.

ёе! а

Следом матрицы а называется сумма элементов, стоящих на главной диагонали:

^ =а11 +а 22 +а33.

0

Обозначения

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

Рассмотрим декартову систему координат с центром в точке Q. Введем следующие обозначения:

ц - единичная нормаль к границе сред, направленная во вторую среду; Е - касательная плоскость в точке пересечения луча с границей; П - плоскость падения, образованная векторами ц и т;

[™я] П

и = |-г - единичная нормаль к плоскости П;

1М

V = [яи].

Векторы ц, и и V формируют базис системы координат (см. рис. 1). Направление единичного вектора оптической оси с выбирается таким образом, чтобы цс > 0.

Особенности распространения луча в анизотропной среде

Электромагнитное поле в любой среде определяется векторами электрического E и магнитного H полей и векторами электрической D и магнитной B индукции. Эти векторы описываются в пространстве и времени уравнениями Максвелла, которые в случае отсутствия токов и свободных зарядов могут быть представлены в следующем виде [5; 6]:

1 ^D 1 ЛВ

rot H =--; rot E =---; div D = 0; div В = 0. (1)

c dt c dt

Связь между векторами E, D и H, В определяется материальными уравнениями, которые описывают электромагнитные свойства среды:

D = sE, В = |H, (2)

где е - тензор диэлектрической проницаемости среды, а ц - тензор магнитной проницаемости среды. Поскольку мы рассматриваем лишь немагнитные среды, то ц = 1, и вектор магнитной индукции В совпадает с вектором напряженности магнитного поля Н.

Рассмотрим решения уравнений Максвелла в виде гармонических плоских волн. Векторы Е и Н в зависимости от времени и точки пространства представляются следующим образом:

Е = Е/Ф, Н = Иов", (3)

где ф = ю^t -— тг^; Е0 и И0 - амплитуды колебаний электрического и магнитного поля;

с - скорость света в вакууме; ю - круговая частота; т = пт - вектор направления распространения фазы; п - показатель преломления (отношение скорости распространения фазы к скорости света в вакууме); п - единичный вектор нормали к фронту распространения фазы; г - радиус-вектор в некоторой декартовой системе координат. При рассмотрении взаимодействия гармонических плоских волн в немагнитных средах (ц = 1) системы уравнений (1), (2) можно привести к виду [5]

Б = -[тН ], Н = [тЕ], Б = вЕ, В = Н. (4)

Из закона сохранения энергии следует, что тензор диэлектрической проницаемости должен быть симметричным [6], т. е. из девяти коэффициентов независимыми являются лишь шесть и е* = е*:

fS Ь11 S12 S ^ &13

s = S21 S22 S23

^S31 S32 S33 J

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

а) в =

0 ^

; б) 8 =

0 0 ^

0

0

; в) 8 =

0

0 ^ 0

8з у

(5)

Значения элементов матрицы, стоящих на главной диагонали, являются значениями главных диэлектрических проницаемостей. Для изотропной среды значения всех главных диэлектрических проницаемостей равны (5а), в результате векторы Б и Е оказываются колли-неарными. Для одноосных сред (5б) одинаковыми являются только два из трех коэффициентов. Для определенности в дальнейшем будем считать, что е1 = е2 =е0 и е3 = ев. Здесь о - обозначение параметра обыкновенного луча (волны) и в - необыкновенного луча. Для оптически положительной одноосной среды ев > ео, в противном случае среда оптически отрицательна. Для двуосных сред (5в) все три значения различаются.

Для изотропной и одноосной среды тензор может быть выражен через значения главных диэлектрических проницаемостей и направления оптической оси. Введем следующие обозначения: с - единичный вектор оптической оси в одноосной среде; п. - коэффициент преломления изотропной среды; по и пв - главные коэффициенты преломления одноосной среды. Известно, что

8. = П2, 8 = П2, 8 = П2. (6)

II* О О ' в в V ^

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

Тогда тензор диэлектрической проницаемости изотропной среды - это диагональная матрица

8 = 8.. (7)

А для одноосной среды верно следующее выражение:

8 = 8о-(8в -8о)(сс). (8)

Можно заметить, что для анизотропных сред, т. е. сред, в которых хотя бы одна из трех главных диэлектрических проницаемостей отличается от двух других, векторы Б и Е оказываются неколлинеарными, это влечет существенные изменения в структуре электромагнитной волны в целом (рис. 2, 3). Вектор направления распространения фазы т пропорционален [БН]. Вектор направления распространения луча (энергии) р пропорционален [ЕН] и в данном случае не совпадает с направлением распространения фронта волны. Электромагнитная волна в анизотропной среде не является строго поперечной, так как существует ненулевая проекция вектора Е на направление распространения фазы т.

Б

Е

т

Рис. 2. Структура электромагнитной волны в анизотропной I ^^ среде. Расположение векторов полей

Рис. 3. Различие направлений распространения энергии р и фазы т в анизотропной среде

Рассмотрим систему (4) более подробно. Исключив из нее вектор Е, получаем

Е = в-Б = -в-тхН.

Подставляя во второе уравнение системы (4), получаем уравнение относительно компонент вектора Н:

(1 + тв-1т) Н = 0. (9)

После некоторых преобразований получаем [5]:

т2твт -т(в -в")т + |в| = 0. (10)

Подставляя представление т = пп, получим биквадратное относительно п уравнение, известное как уравнение нормалей Френеля для изотропной, одноосной и двуосной сред:

п4пвп -п2п(в - в")п + |в| = 0. Для одноосных сред, учитывая (8), получаем уравнение на компоненты вектора распространения фазы:

(т2-в0)(твт-в0 ве ) =

Легко видеть, что уравнение разбивается на два уравнения.

Для обыкновенного луча для вектора распространения фазы то имеем

т2 - во = 0. (11)

Для необыкновенного луча для вектора распространения фазы те получаем

тевте -в0ве =0. 12)

Для изотропной среды уравнение (10) принимает следующий вид:

т2 -в,. = 0.

(13)

Определим направления колебаний векторов Е и Н в одноосной среде. Для этого подставим в уравнение (9) представление для тензора (8), тогда получаем [5]:

в0 - т2 + (ве-в°) ([тс] •[тс])) Н = 0.

V ве у

Тогда вектор Н для одноосных сред можно представить в виде

где Ь0 =[т0 [тос],

Но = Л0Ь0,

Не = ЛА,

(14)

(15)

где Ь е = [т ес].

Здесь Ло и Ле - вещественные константы, а то и те - векторы распространения фазы обыкновенного и необыкновенного луча соответственно.

Из системы уравнений (4) следует, что Е = -в 1 [шИ ], откуда получаем значение вектора Е для электромагнитной волны, распространяющейся в одноосной среде:

Ер = А еР, (16)

где ер = -в-1 [шрЬр] и Ре {о,е).

Из (14)-(16) следует, что все распространяющиеся в прозрачных одноосных средах лучи линейно поляризованы. Заметим, что в поглощающих средах лучи поляризованы слабо эллиптически [6].

Плотность потока энергии электромагнитного поля определяется как сумма энергий электрического и магнитного 2,т полей [5]:

Е = Е+Е £ (ЯеЕ)(ЯеР) Е =(КеИ)(КеБ)

8л 8л

Учитывая уравнения (4) и тот факт, что вектор ш является вещественным, получаем

1 2

^ = —(е И) . (17)

4л

Вектор плотности потока энергии (вектор Пойнтинга) определяется выражением

Р = ^-[е Е,Яе И ]. (18)

Вектором распространения энергии будем называть вектор

Р я

Р =4 = -, (19)

СС, 5

где я - единичный вектор направления распространения энергии; 5 - лучевой коэффициент преломления в среде; с - скорость света в вакууме.

Умножив уравнение (18) скалярно на ш и используя формулы (4) и (17), получаем

£=- шР.

с

Тогда, учитывая (19), получаем

шр = 1. (20)

Вектор распространения фазы ш будем называть Ф-нормализованным, если он удовлетворяет уравнению нормалей Френеля для той среды, в которой распространяется волна (10)-(13). В этом случае его длина равна коэффициенту преломления п (отношение скорости распространения фронта волны в данном направлении к скорости света в вакууме), т. е. |ш| = п. Вектор лучевой рефракции будем называть Ф-нормализованным, если его длина обратна лучевому коэффициенту преломления, т. е. |р| = 1/ 5.

Для Ф-нормализованных векторов ш и р можно выписать уравнения, позволяющие вычислять в явном виде ш по вектору р, и наоборот. Для этого необходимо подставить представления векторов Е и И (14)-(16) в определение вектора Пойнтинга (18) и учесть (20). Опустим вывод соответствующих соотношений, который достаточно громоздок (см. [5]). Для изотропных сред:

р, = вчш,, (21)

Для одноосных сред:

m, =sp,. (22)

Р. = ^, Р. = ^, (23)

s s s

. . .

m = s p , m = s_1p s s . (24)

Отметим, что в формулах (21)-(24) предполагается, что тензор s - это тензор (7) или (8), тензор среды, в которой распространяется луч.

Другой важной характеристикой луча является его интенсивность - среднее по времени количество энергии, приносимой им на единичную площадку, перпендикулярную направлению его распространения. Она вычисляется по формуле [5; 6]:

I = \mean ( Рц )|, (25)

где Р - вектор Пойнтинга; я - единичный вектор нормали к площадке, а теап() - усреднение по времени.

Упростим выражение для каждой среды. Поскольку векторы Н и Е вещественны в одноосной среде, то выражение для вектора Пойнтинга (18) приводится к виду

Р =-^-[ЕН ]. (26)

4—

Подставляя уравнение (26) в (25) и учитывая (14)-(16), получаем выражения для интенсивности луча, распространяющегося в одноосной среде:

1р= 8-лч, (27)

[ерЬр] я| и ре-^е}.

где ир =

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

Е, = Дт В, [п,т , (28)

Н, =[т,.Е,.] = Л, [т,тх]-пгВгт1, (29)

где п, - единичный вектор направления распространения луча; п - коэффициент преломления среды; Л, и В, некоторые комплексные константы. Тогда уравнение (25) сводится к следующему виду:

I' = ¿(1Л2 +1 В\2 )и, (30)

где и, = Кя|.

Модель взаимодействия света с границей двух сред

Геометрические соотношения для отражения и преломления света могут быть получены из условия равенства тангенциальных компонент вектора электрического поля Е и равенства векторов магнитного поля Н [5]:

[(Е1 - Е2) я ] = 0, (31)

Н1 - Н 2 = 0. (32)

Цифрами обозначены номер среды (первая или вторая). При падении волны на границу двух сред в общем случае возникает V отраженных и т преломленных волн. В этом случае уравнение (32) приводится к виду

Н 0 +(Н Г1 +... + Н „)-( Нд +... + Н,) = 0.

Здесь индекс 0 обозначает падающую волну, а г и ^ в индексах - отраженные и преломленные волны соответственно. Учитывая уравнение (3), получаем

И0е1ф° +((г1е+... + НуеЩг")-(иае,ф" +... + Ине,ф'х) = 0. (33)

Пусть уе{0, г1,..., гу, И,..., ¿х}. Поскольку фу=юу^ I - -1 т уг | и уравнение (33) должно выполняться во всех точках границы и во все моменты времени, получаем, что юу не изменяется во время взаимодействия и (т 0 - т у) г = 0.

Принимая во внимание, что уравнения выполняются лишь в точках границы между средами, т. е. гц = 0, тогда из произвольности радиус-вектора г вытекает, что (ш0 -шу)|. Это значит, что концы всех векторов шу на рис. 4 заканчиваются на одной вертикальной линии. Из коллинеарности векторов вытекает, что [ш 0ц] = |ш у ц ^.

Введем векторы а = [ш у ц ], g = [ца], и пу = ш у ц, тогда вектор распространения фазы ш у можно представить в виде

ш т= g + Яу Ч. (34)

Зная Ф-нормализованный вектор распространения фазы ш0 падающей волны, вычислим вектор g:

g = ш 0-(ш0я ) я. (35)

Легко видеть, что вектор g является проекцией вектора распространения фазы ш0 на плоскость границы раздела сред Е, причем важно то, что он одинаков для всех векторов распространения фазы падающей и всех порожденных волн. Этот факт является ключевым для получения геометрических соотношений взаимодействия луча света с границей двух сред. Также важно отметить, что из (34) следует, что все векторы ш у лежат в одной плоскости в

отличие от соответствующих им векторов р у (21)-(24).

Для того чтобы вычислить Ф-нормализованные векторы распространения фазы для всех отраженных и преломленных лучей, подставим (34) в уравнение нормалей для первой и второй среды соответственно. Если среда является изотропной, то ее уравнением нормалей будет (13), тогда

ш,2 = (+я ц )2 = я2 + 2Я ёц+ё2 = в2. (36)

Для одноосной среды получаем уравнение

ш 2 =(+я0 ц )2 = я2 + ^вц + ё2 =в2. (37)

шевше = ( + Яец)В ( + ^) = + 2Яе () + §В§ = ВоВе . (38)

Решая два уравнения (36)-(38) (для каждой из сред), в общем случае получим до 8 комплексных значений Яр (для границы двух одноосных сред: по 4 корня для каждой среды).

Для каждого значения Яр, где Р е {,,о,е), по формуле (34) вычисляем вектор шр . Используя формулы (21) и (23) в зависимости от типа среды, вычисляем вектор рр для каждого вектора шр .

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

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

Следующей задачей является задача получения амплитудных соотношений между падающим и порожденными лучами. Для этого необходимо разрешить систему уравнений (31)-(32), используя вычисленные вектора распространения фазы ш у и представление векторов Е и И для каждой из сред: формулы (14)-(16) - для одноосных, (28) и (29) - для изотропных. Так как в формулах (28), (29) вектор ш1 - произвольный единичный вектор, перпендику-

Изотропная

среда

mi mi •

g П1

\ 1 \ \ \ \ // ✓ / ✓ / * / ч/ /ч ^ , m t > П2 / / / / ✓ / / / ^ << I Изотропная ,4 среда

Изотропная среда m ii' g

\ I \ \ 4 \ \4 me) /f s// ' /I V . mio > ' 1 / / / / ' / / .-// \ \ \ \ с Одноосная q среда

Рис. 4. Все возможные комбинации первой и второй среды и типов падающего луча. Окружности и эллипсы на картинках - это срезы плоскостью П поверхности распространения фазы. Точки пересечения (жирные точки) перпендикуляра (тонкая вертикальная линия) с окружностями и эллипсами показывают векторы распространения фазы падающего и порожденных лучей. Вектор g - это проекция всех векторов распространения фазы на границу раздела сред Е. Оптические оси в одноосной среде обозначены c, cb c2. Подстрочные индексы i, o, e, r, t обозначают отношение вектора к изотропному, обыкновенному, необыкновенному, отраженному и преломленному лучу соответственно. Выражение вида mtoe обозначает «transmitted - ordinary - extraordinary», т. е. вектор распространения фазы преломленного необыкновенного луча, возникшего при падении обыкновенного луча на границу сред

лярный направлению распространения луча, т. е. вектору шу ру, то положим ш1 = и. Тогда система (28), (29) приобретает вид

Е, = А и + В [п,и ], (39)

И, = А [ш,и] + пВ^и. (40)

Далее система (31) и (32) сводится к виду

[(Ее +(ЕГ1 +... + Е„)-( Ея +... + Е,)) ц ] = 0,

И 0 +(И Г1 +... + И „)-( Ия +... + И ,) = 0.

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

([Е^]... + [Е^ц ])-([Е ^ ] +... + [Е „л]) = [Е0ц ], (41)

(Ид +... + И,)-( И Г1 +... + И „) = И 0. (42)

Для решения системы (41), (42) необходимо подставить представления векторов Е и И (14)-(16), (39), (40) в нее и умножить каждое из уравнений скалярно на два произвольных взаимно перпендикулярных вектора, лежащие в плоскости Е (мы используем векторы и и V).

Учитывая то, что V = [ци] и и = ^ц] (см. рис. 1), получаем систему из четырех уравнений:

(и +... + Е,и ) - (Ели +... + Е„ и ) = Е0и, (43)

( V +... + Е, V) - (ЕГ1V +... + Е„ V) = Е0 V, (44)

( Ияи +... + И ) - ( И Йи +... + И „и ) = И 0и, (45)

( И я V +... + И, V) - ( И Г1V +... + И „ V) = И 0 V. (46)

Поскольку значения в правых частях уравнений определяются падающей волной, будем считать их известными. Необходимо разрешить систему относительно величин Аг и Ву, кроме у = 0. Система оказывается однозначно разрешимой только в случае, если число свободных переменных равно числу уравнений. Если во взаимодействии с границей участвует изотропный луч (луч, распространяющийся в изотропной среде), то ему соответствуют две свободные переменные в уравнении: А, и В, в (39) и (40), а для плоско-поляризованного луча, распространяющегося в одноосной среде, переменная одна: Ао, Ае в (14)-(16). Это соответствует действительности, так как число порожденных при взаимодействии лучей и их поляризация должны быть такими, что число свободных переменных строго равно четырем. При падении луча на границу двух анизотропных сред образуется в общем случае четыре линейно-поляризованных луча, на границе двух изотропных сред образуется два изотропных луча, а на границе изотропной и одноосной среды образуется два плоско-поляризованных луча и один изотропный луч. Таким образом, в любом случае образуется четыре свободных переменных и уравнение однозначно разрешимо.

Рассмотрим пример составления такой системы уравнений. При падении луча из изотропной среды на границу с одноосной средой получаем систему уравнений:

((] +... + М- [Е^НЕ^], (47)

(Ио + Ие)-И г = И 0. (48)

Здесь индексы о, е, г, и 0 обозначают обыкновенный преломленный, необыкновенный преломленный, изотропный отраженный и изотропный падающий соответственно.

Пусть коэффициент преломления изотропной среды равен п. Подставляя соответствующие разложения векторов Е и И ((14)-(16) и (39), (40)), получаем:

(Ао М] + Ае ^ ]) - (Аг [иц ] + Вг [[п г и ]) = А N ] + В0 [[п0и ], (49)

( АоЬ о + АеЬ е) - ( Аг [ш г и ] - Вг (пи )) = А [шеи ] - В0 (пи ). (50)

Умножая оба уравнения на и и V, а затем решая полученную систему из четырех уравнений, получаем значения величин Ао, Ае, Аг, Вг. Далее мы можем вычислить векторы Е и И для всех отраженных и преломленных лучей (см. (14)-(16) и (39), (40)).

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

Замечание. Отметим, что в данной работе мы не рассматриваем изменение фазы луча между падающей волной и порожденными волнами.

Представление поляризации света

В данной работе для представления поляризации света мы используем подход, основанный на матрицах когерентности и их модификаторах [6; 7]. Рассмотрим кратко этот подход, необходимый для алгоритма расчета лучей.

Введем некоторую декартовую систему координат, связанную с лучом, с базисом у, z}.

Вектор z направлен по направлению распространения луча. Комплексным вектором Джонса называется вектор [7]:

Е =

е е ы+5,) Е е'(юг+5у)

здесь ю - круговая частота колебаний; 5, и 5у - начальные фазы колебаний; Е, и Еу - амплитуды проекции вектора E на оси x и у соответственно. Реальным значением вектора электрического поля E будет ЯеЕ. Матрицей когерентности называется комплексная матрица следующего вида [7]:

J =

(Е X Е

Е, (г)

Еу (г)

х(Е* (г) Е* (г ))) =

( ЕхЕх^ ( Е,Е*) ^ ЕуЕ*) ( ЕуЕу )

, ^'ух •^уу ,

Здесь (•) обозначает взятие комплексного сопряжения; ^ - математическое ожидание; (•У - транспонированное комплексное сопряжение. Матрица когерентности в совокупности с системой координат у, z}, связанной с лучом, полностью описывает состояние поляризации луча.

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

I = Jt = J„ + Jyy = ( Е,Е*) + {ЕуЕ*). (51)

Пусть при некотором взаимодействии электрическое поле луча изменилось в соответствии с законом

Е' = м Е,

где М - это комплексная матрица, описывающая изменение электрического поля луча, будем называть ее матрицей-модификатором матрицы когерентности. Новая матрица когерентности J' вычисляется как

J ' = ( Е х Е = ( м Е х ЕМ = м( Е х Е*)мt = мт ^ (52)

В данной работе для вычислений удобно использовать различные системы координат для каждого типа среды. Для одноосной среды луч поляризован линейно, и удобно использовать следующую систему координат (см. (14)-(16)):

ев Ьк

У =

х =

ъ =

Рр

|Рр|

(53)

|еР I 1иР1

где Ре {о,е}.

Для изотропной среды ввиду отсутствия в ней какого-либо выделенного направления и произвольности поляризации луча, распространяющегося в ней, такую систему координат мы можем определить только на границе сред (см. (39), (40)), а именно:

x = ^ у = [™ ], z = п.

(54)

В некоторых случаях необходимо поворачивать систему координат, связанную с лучом. Обозначим оси новой системы координат {х',у',г'}, тогда для пересчета матрицы когерентности можно использовать следующую матрицу-модификатор:

Mrotate

xx

УХ

Л

ХУ УУ у

(55)

Алгоритм расчета взаимодействия луча

с границей между двумя средами

Рассмотрим алгоритм расчета взаимодействия луча с границей двух изотропных и одноосных немагнитных прозрачных сред. Назовем лучом следующее множество:

Я = {г0, г,, т, 3, СБ ,ТР, X, ЕЫ}, (56)

здесь г0 - начальная точка луча; г^ - вектор направления распространения луча; т - вектор направления распространения фазы; X - длина волны луча в вакууме (она однозначно определяет круговую частоту ю); 3 - матрица когерентности луча; СБ - система координат, связанная с лучом; ТР - тип луча (изотропный, обыкновенный, необыкновенный) и булевский флаг ЕЫ (равен «истина», если луч Ф-нормализован).

Пусть луч Я падает на границу двух сред. Необходимо вычислить эти же параметры для всех отраженных и преломленных лучей. Обозначим величины в (56), соответствующие порожденному лучу надстрочным символом «*».

Алгоритм состоит из следующих шагов.

1. Ф-нормализация падающего луча. Этот шаг выполняется только, если ЕЫ равен «ложь»:

а) с помощью уравнений (22), (24) вычисляем направление вектора распространения фазы т', используя в качестве р вектор г,. Вектор т' будет параллелен т, но не Ф-нормализован;

б) для того чтобы Ф-нормализовать т', подставим кт' в качестве т в (11)—(13) для первой среды и падающего луча типа ТР. Полученное уравнение имеет несколько корней, но только один может быть положительным и вещественным, выбираем его в качестве значения к. Тогда т = кт';

в) вспоминая соотношение (20), вычисляем вектор распространения энергии р = га/тг, .

2. Преобразовать матрицу когерентности в систему координат, необходимую для вычислений. В зависимости от среды, типа луча ТР и границы, на которой мы рассчитываем взаимодействие, определяем нужную систему координат (53) или (54). Применяя матрицу-модификатор (55) по правилу (52), преобразовываем матрицу когерентности в новую систему координат.

3. Вычислить проекцию g вектора т падающего луча на плоскость границы раздела сред Е по (35).

4. Рассчитать векторы т для отраженных и преломленных лучей. Сначала необходимо вычислить такие значения п , которые удовлетворяют уравнениям нормалей для первой и второй сред (см. (36)-(38)). Затем по формуле (34) вычислить векторы т .

5. Вычислить типы порожденных лучей ТР в зависимости от уравнения, используемого на шаге 4.

6. Выбрать векторы т и вычислить векторы р :

а) для каждого вещественного вектора т по формулам (21), (23) вычисляем вектор р . Выбираем такие векторы т , для которых выполняется одно из условий:

1) для отраженных лучей р ц < 0, и вектор т получен из уравнения нормалей для первой

среды; * *

2) для преломленных лучей р ц > 0, и вектор т получен из уравнения нормалей для второй среды;

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

7. Составить амплитудное уравнение и разрешить его. Для каждого из векторов т , выбранных на шаге 6, вычисляем коэффициенты амплитудного уравнения (см. (43)-(46) и пояснения к ним). Для формирования уравнения можно использовать следующий алгоритм:

а) используя (14)-(16), для всех порожденных лучей, распространяющихся в одноосной среде, вычисляем один коэффициент-столбец вида

Г

Ь„

(х ^

V уру

р у

где ре {о,в, +,-}.

Используя (39) и (40), для всех порожденных лучей, распространяющихся в изотропной среде, вычисляем два коэффициента-столбца вида

(Х1 ^([ияП и (Х2^Мя?

V У1 У

,Ки].

и

V У2 У

V —Пи У

б) в соответствии с системой уравнений (41), (42) берем коэффициенты-столбцы для преломленных лучей с положительным знаком, а для отраженных с отрицательным знаком и формируем матрицу 2 х 4 (2 строки, 4 столбца), соответствующую левой части системы уравнений (41), (42);

в) аналогично вычисляем один или два коэффициента столбца для падающего луча и формируем правую часть уравнения;

г) в результате получаем матрицу 2 х 5, если первая среда - одноосная, и матрицу 2 х 6, если - изотропная;

д) умножая каждую из строк на и и V скалярно, получаем матрицу 4 х 5 или 4 х 6 амплитудного уравнения;

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

Продолжим рассмотрение примера, приведенного выше (47)-(50). Для преломленных обыкновенного и необыкновенного луча получаем следующие коэффициенты-столбцы, см. (14)-(16):

' Хо ^ (М П ( Хв ПМР

и

V Ув У

Для отраженного луча:

Для падающего луча:

( Х11_([ия]'

[т г и]

V

V 1г1 У

(х ^

Лг 2

V

V 1Г 2 У

[[п г и ]я ]'

-пи

( X 01 КИП ( Х02 Н[Мя ]'

V

V Х01 У

,[т ои ]

V

V *02 У V

-пи

Из коэффициентов-столбцов составляем матрицу 2 х 6:

Х0 Хв —Хг1 —Хг 2 Х01 Х02

V Уо

У —У —У

1в 1г1 1г 2

101 02 У

Умножая строки матрицы скалярно на векторы и и V, получаем следующую матрицу: (Хои Хви —Хяи — Х г 2и Х01и Х02и ^ Х о V Х в V — Х^ — Х^ XolV Х02 V У0и Уви —УЙи —Уг2и У01и У02и

\У0 V Ув V —УгlV —Уг2V УolV У02 V у Полученная матрица является формой системы уравнений (43)-(46):

(МММ4){А0ААВГ)т _ (м6)(()т,

где М7 - 7-я колонка матрицы.

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

а) обнуляем матрицу-модификатор Мв;

б) копируем элементы из 5-го и 6-го столбцов (если 6-й есть) матрицы амплитудного уравнения в столбцы матрицы-модификатора М0 со строки с номером равным номеру первого коэффициента-столбца в левой части уравнения (столбцы 1-4), соответствующего порожденному лучу. Количество копируемых строк равно количеству коэффициентов-столбцов в уравнении, соответствующих порожденному лучу (одна для луча, распространяющегося в одноосной среде, и две для изотропного луча);

в) применяем рассчитанную матрицу-модификатор М0 к матрице когерентности падающего луча по правилу (52), получая матрицу когерентности порожденного луча:

J ' = МВ1М1; (58)

г) для того чтобы формула (51) оставалась верной в любой среде, необходимо умножить матрицу когерентности J ' на коэффициент пропускания или отражения без учета амплитуд векторов E и ^ так как они уже учтены в преобразовании (58). Обозначим этот коэффициент / тогда по формулам (27), (30) получаем

Г = Ца

и '

где а,ре {0,1,о,е}, причем в соответствует падающему лучу, а а - порожденному;

д) тогда J* = /А ' - матрица когерентности порожденного луча.

Продолжим рассмотрение примера. Для границы между изотропной и одноосной средами после приведения левой части матрицы, полученной на шаге 7, к единичному виду получаем

(1 0 0 0 М15 М16 ^

0 1 0 0 М 25 М 26

0 0 1 0 М35 М36

V 0 0 0 1 М 45 М46 У

Это решение системы (57), поскольку

{лолелгвг )т = (ММ2М3М4 у1 (М5М6 )((в0 )т =

( М15 М16 ^

М 25 М 26 ( Л 1

М35 М36 V в0 у

V М 45 М46 у

Тогда матрицы-модификаторы для обыкновенного, необыкновенного и отраженного луча соответственно равны:

МВо =

( М„

М16 1

МВе =

( М 2

М26 1

МВг =

( М35 V М 45

М36 1

М46 У

Далее вычисляем значение / для каждого из лучей:

/о =

М о ]д

/е =

[^е ]д

m oq

/ =

m г q

= 1.

9. Вычислить и установить новые системы координат С$ для каждого из порожденных лучей, см. (53), (54).

10. Установить свойства луча для каждого порожденного луча: г0 = Q (точка пересечения луча с границей сред) и г* = p. Также устанавливаем значение ЕЛ* «истина», так как p* и m Ф-нормализованы.

11. Конец алгоритма.

Эксперименты

Лучшим тестом для проверки алгоритмов фотореалистической визуализации является сравнение рассчитанного изображения с фотографией. Именно такой подход использован, например, в работах [3; 8; 9]. Во всех случаях авторы столкнулись с рядом неизбежных проблем, связанных с совмещением параметров реальной и виртуальной сцены: оптических характеристик и геометрий объектов сцены, положений, геометрий и спектров источников света, характеристик и ракурса камеры. Тем не менее только этот подход может подтвердить корректность алгоритма фотореалистической визуализации кристаллов, поэтому такое сравнение было выполнено и в этой работе.

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

Детальное устройство реальной сцены и примерные размеры приведены на схеме на рис. 5 (относительные размеры не соблюдены), а на рис. 6 - фотография реальной сцены, включая расположение фотоаппарата, установленного на штативе рядом со сценой.

300 х 300 мм

120 мм

белая диффузная бумага клей

2 белые люминесцентные лампы 0 11 мм

100 мм

черная диффузная бумага

черная диффузная бумага

прозрачное цветные нет щелей стекло 2 мм чернила

87 х 87 мм

2 деревянные стойки

Ж -Ж- ^ 1

2 деревянные стойки

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

Виртуальная сцена состоит из прямоугольного текстурированного источника света и кристалла, расположенного параллельно ему на минимальном расстоянии, достаточном для отсутствия геометрических ошибок при вычислениях. Кристалл представлен в виде набора граней, форма каждой из них подобна форме грани реального кристалла. Для того чтобы добиться пропорциональности граней между реальным и виртуальным кристаллом были выполнены следующие шаги: 1) каждая из граней реального кристалла была отсканирована; 2) на каждом из отсканированных изображений были вручную выделены точки, соответствующие вершинам кристалла; 3) по разметке получена двухмерная геометрия граней; 4) восстановлена трехмерная геометрия кристалла. Оптические свойства виртуального кристалла соответствуют оптическим характеристикам кальцита. Для того чтобы определить направление оптической оси кристалла, перед сканированным на двух больших гранях кристалла были отмечены две точки таким образом, чтобы при просмотре через кристалл они сливались в одну точку (т. е. двоение не наблюдается). Оптическая ось кристалла проходит через эти две точки. В процессе восстановления геометрии рассчитываются трехмерные позиции этих двух точек, т. е. определяется направление оси.

В качестве виртуальной камеры используется камера из алгоритма классической обратной лучевой трассировки [1]. В качестве реальной камеры была выбрана Canon 450D с объективом EF-S18-55 мм f/3.5-5.6 IS. Чувствительность матрицы была выбрана минимальной (ISO100), чтобы уменьшить шумы, а значение диафрагмы - максимальным (F36), чтобы уменьшить апертуру объектива (апертура виртуальной камеры равна 0). Ракурс виртуальной камеры совмещался с ракурсом реальной камеры вручную. Выдержка должна быть существенно больше периода мигания люминесцентной лампы, для фотографии на рис. 7 она составляет 13 с. Фотографирование производится при отсутствии внешних источников освещения.

Существует ряд проблем при совмещении спектров текстур: 1) используемые лампы обладают неизвестным спектром излучения; 2) спектр пропускания бумаги и цветных чернил неизвестен; 3) объектив вносит неизвестные искажения в спектр проходящего через него света; 4) чувствительность матрицы фотоаппарата к отдельным компонентам спектра также неизвестна. Для того чтобы минимизировать различия между цветами текселей текстур на фотографии и рассчитанном изображении, цвета виртуальной текстуры выбираются по размытой фотографии реальной текстуры (без кристалла), сделанной с тем же ракурсом и такой же выдержкой. Цвета виртуальной текстуры хранятся в виде RGB-цветов, а для вычисления подходящего спектра по RGB-цвету был использован поход, предложенный в работе [10].

На рис. 7 приведена фотография кристалла, а на рис. 8 - рассчитанное изображение виртуальной сцены. Изображение было рассчитано алгоритмом обратной спектральной рекурсивной лучевой трассировки. При расчете использовались спектры с равномерно распределенными 21 сэмплами на видимом диапазоне с 380 до 780 нм.

Рис. 7. Фотография кристалла кальцита

Рис. 8. Рассчитанное изображение

r'S- 1

Г L Й 4 Й M A

1 1 1 4 HJ

I

и L_ra

Рис. 9. Попиксельная разница между фотографией и рассчитанным изображением

На рис. 9 приведена абсолютная пиксельная разница между фотографией и рассчитанным изображением. Видно, что разница между изображениями подтверждает корректность алгоритма. Имеющаяся разница объясняется несколькими причинами, прежде всего отличием ракурсов реальной и виртуальной камер, так как идеальное совмещение ракурсов практически невозможно без использования алгоритмов калибровки камеры, что хорошо видно по изображениям, приведенным в статьях [3; 8; 9]. В случае сравнения изображений прозрачных монокристаллов ручной подбор ракурса еще более сложен, так как малейшее изменение ракурса камеры приводит к сильным изменениям в изображении, что особенно хорошо заметно на вторичных гранях (видимых через другие грани). Во-вторых, различием между цветами текселей физической текстуры и текстуры на рассчитанном изображении, которое объясняется неравномерной структурой бумаги, что хорошо видно на разнице между изображениями. В-третьих, отношение размера текстуры к размеру кристалла немного отличается в реальной и виртуальной сценах.

Выводы

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

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

1. Whitted T. An Improved Illumination Model for Shaded Display // Commun. ACM. 1980. Vol. 23. No. 6. P. 343-349.

2. Guy S., Soler C. Graphics Gems Revisited // ACM Transactions on Graphics: Proc. of the SIGGRAPH Conference. 2004. Vol. 23. No. 3. P. 231-238.

3. Weidlich A., Wilkie A. A. Realistic Rendering of Birefringency in Uniaxial Crystals // ACM Transactions on Graphics: Proc. of the SIGGRAPH Conference). 2008. Vol. 27 (1). P. 6:1-6:12.

4. McClain S. C., Hillman L. W., Chipman R. A. Polarization Ray Tracing in Anisotropic Optically Active Media. 2. Theory and Physics // Journal of Optical Society of America. 1993. Vol. 10. No. 11. P. 2383-2393.

5. Федоров Ф. И., Филиппов В. В. Отражение и преломление света прозрачными кристаллами. Минск: Наука и техника, 1976.

6. БорнМ., Вольф Э. Основы Оптики. М.: Наука, 1973.

7. Tannenbaum D. C., Tannenbaum P., Wozny M. J. Polarization and Birefringency Considerations in Rendering // Proc. of SIGGRAPH'94 (Orlando, Florida, July 24-29). 1994. P. 221-222.

8. Goral C. M., Torrance K. E., Greenberg D. P., Battaile B. Modeling the Interaction of Light between Diffuse Surfaces // Proc. of the XI Annual Conference on Computer Graphics and Interactive Techniques. N. Y.: ACM, 1984. P. 213-222.

9. Meyer G. W., Rushmeier H. E., CohenM. F., Greenberg D. P., Torrance K. E. An Experimental Evaluation of Computer Graphics Imagery // ACM Transactions on Graphics. 1986. Vol. 5. P. 30-50.

10. Glassner A. S. How to Derive a Spectrum from an Rgb-Triplet // IEEE Comput. Graph. Appl. 1989. Vol. 9. Р. 95-99.

Материал поступил в редколлегию 22.06.2011

V. A. Debelov, D. S. Kozlov

A LOCAL MODEL OF LIGHT INTERACTION WITH ISOTROPIC AND UNIAXIAL TRANSPARENT MEDIA

The paper is devoted to the derivation of a local model of light interaction with a boundary of two transparent crystals (media) with different optical characteristics (isotropic or uniaxial). Refracting and reflecting properties of each medium are defined by the direction of the optical axis (if any) and refraction spectra. Local model describes a set of rules to calculate directions, intensities and polarization states of all outgoing rays by corresponding attributes of an incoming ray and optical properties of adjacent media. The developed model operates with the term of polarized ray, i.e. geometrical ray with linked intensity and polarization state attributes. The suggested algorithm for calculation of rays allows calculating properties of all reflected and refracted rays (up to 4 rays). In the paper, isotropic and uniaxial media are considered in a similar manner. The correctness of the model is validated by comparison of a photo of a real uniaxial crystal with the corresponding computed image.

Keywords: isotropic crystal, uniaxial crystal, optical axis, light polarization, optical dispersion, birefringence, photorealistic rendering, ray tracing.