CC BY
18
Вычисление параметров объемного рассеяния на основе измерений рассеяния света
Поздняков С.Г., Ершов С.В., Волобой А.Г., ИПМ им. М.В. Келдыша РАН pozdnyakov@gin.keldysh.ru, sergey_65@mail.ru, voloboy@gin.keldysh.ru
Аннотация
Статья посвящена описанию процедуры определения параметров, описывающих распространение видимого света в материалах с объемным рассеянием. Показано, как с помощью несложных вычислений такая процедура реализуется на основе простой лучевой трассировки. В качестве иллюстрации применения в статье подробно описаны расчеты для двух существенно различных материалов с объемным рассеянием.
1 Введение
В современной оптической индустрии широко используются устройства и элементы, произведенные из материалов с объемным рассеянием света. К ним можно отнести различные диффузоры, элементы подсветки шкал, выходные линзы светодиодов и т.п. Для проектирования таких устройств необходимо с достаточной точностью знать оптико-физические свойства материалов с объемным рассеянием. Как правило, при фабрикации таких материалов параметры, описывающие объемное рассеяние не могут быть заданы априори. Вместо этого производители материалов изготавливают достаточно широкий набор таких материалов, из которых заказчик может выбрать наиболее подходящие.
Следует отметить, что подбор материалов с объемным рассеянием методом проб и ошибок представляет собой достаточно дорогую и трудоемкую процедуру, особенно в случае сложной геометрии проектируемых оптических элементов. Таким образом, определение параметров объемного рассеяния, исходя из результатов оптических измерений для заданного материала, представляет собой важную и актуальную задачу.
Эта задача была решена для использования в программном комплексе моделирования освещенности и синтеза реалистичных изображений Ьишюер^ разработанном в Институте прикладной математики им. М.В. Келдыша РАН [Жданов, Потемин, Галактионов, Барла-дян, Востряков, Шапиро, 2011].
2 Методы получения параметров объемного рассеяния
Непосредственные измерения параметров объемного рассеяния достаточно сложны и трудоемки и, поэтому, не всегда возможны. Часть информации о параметрах может быть получена при измерении рассеяния света достаточно тонкой (чтобы избежать более чем однократного рассеяния) пластинкой материала. Однако, в таком случае "видна" не вся фазовая функция (угловое распределение) из-за эффекта полного внутреннего отражения: свет, рассеянный под большими углами, не может покинуть пластинку из-за этого эффекта. В принципе можно выполнить измерения, при которых влияние полного внутреннего отражения исключается. В этом случае к обеим поверхностям пластинки иммерсионно приклеиваются две прозрачные полусферы, через которые проходит и освещающий световой пучок, и рассеянный свет. Однако такие измерения достаточно сложны и могут быть непосредственно выполнены не на всякой измерительной установке. Кроме этого, изготовление тонких пластинок с толщиной порядка длины свободного пробега (часто меньше чем -0.05^0.1 мм) представляет собой дополнительную проблему.
Вычисление параметров объемного рассеяния, исходя из свойств компонентов материала, также практически невозможно из-за очень большой сложности процедуры фабрикации и небольшого количества точно решаемых моделей рассеяния (по сути дела пригодно только решение Mie для рассеяния на диэлектрическом шаре [Mie Scattering]).
Альтернативой служит процедура измерения рассеяния света для одной или нескольких пластин (разной толщины) исследуемого материала, а затем вычисление параметров рассеяния с помощью решения обратной задачи рассеяния. В этом случае производится расчет рассеяния света для пластин таких же толщин с параметрами рассеяния, которые варьируются таким образом, чтобы достичь наилучшего согласия между измеренными и
вычисленными зависимостями. В литературе опубликовано несколько статей, использовавших такой подход. Практически во всех из них исследуются различные биологические ткани или эмульсии [см., например Тучин]. В большинство этих статей в качестве фазовой функции используется хорошо известная индикатриса Неупеу-Огее^ет или ей подобные. Анализ результатов, полученных в этих статьях, позволил сделать вывод о достаточной применимости такого подхода. Однако следует отметить, что результаты, полученные в некоторых из этих статей, представляются не бесспорными. В частности, использование суммы двух индикатрис Неупеу-Огееп81ет для описания прямого и обратного рассеяния представляется не совсем корректным.
В дополнение ко всему сказанному следует добавить, что ряд ведущих мировых оптических фирм (например, ОРТК [ОРТК]) производят определение параметров объемного рассеяния для предоставляемых материалов. Впоследствии эти параметры используются при компьютерном моделировании и проектировании оптических элементов и устройств в их программном продукте. Вероятно, что в этом случае используется процедура, близкая к описанной выше. Описание процедур извлечения параметров, используемых этими фирмами, в литературе не опубликовано, так как, видимо, составляет коммерческую тайну.
3 Измерения и анализ измеренных данных
В нашем случае использовался подход, основанный на "извлечении" параметров объемного рассеяния, исходя из результатов измерения рассеяния света пластинами исследуемого материала. Измерения проводились на существующей светоизмерительной установке [Волобой, Галактионов, Ершов, Летунов, Потемин, 2006]. При этом предполагалось, что измерения будут выполнены без значительных отклонений от штатного режима. Единственное отличие состояло в том, что были выполнены дополнительные измерения интегральных коэффициентов рассеяния и пропускания, которые позволили несколько повысить точность.
Измерения проводились для двух материалов А и В, изготовленных из молочного органического стекла (коэффициент преломления объемлющей среды равен 1.49), которые бы-
ли представлены в виде трех пластин различной толщины для каждого материала. Толщины пластин для материала А были равны 1.53мм, 1.98мм и 3.15мм, для материала В -2.05мм, 3.06мм и 3.63мм. Для всех 6 пластин были измерены двунаправленные функции отражения (ДФО) для 5 углов падения (10°, 20°, 30°, 45°, 60°) и двунаправленные функции пропускания (ДФП) для 6 углов падения (0°, 10°, 20°, 30°, 45°, 60°) при освещении с обеих сторон каждой пластины. Измерения проводились в спектральном представлении в видимом диапазоне с шагом по длине волны в 10нм. Интегральные коэффициенты отражения и пропускания были уточнены с помощью дополнительных измерений.
Прежде всего, в случае пластин, изготовленных из материала А, была обнаружена флюоресценция для длин волн £ 420 нм. Вследствие этого все расчеты для этого материала выполнялись для длин волн, начиная с 430нм.
Поверхности пластин не были идеально отполированы, а были шероховатыми. При визуальном анализе не были обнаружены ни зеркальное отражение, ни прямое пропускание. Измерения ДФО показали, что микрорельефы поверхностей различных пластин и даже поверхностей одной пластины несколько отличались. Для того чтобы выяснить роль микрорельефов при формировании ДФО были выполнены дополнительные иммерсионные измерения для части пластин и меньшего числа углов падения и длин волн. Результаты этих измерений позволили с определенностью сделать вывод о том, что узкий пик у части ДФО и ДФП около направлений зеркального отражения и прямого пропускания возникает за счет влияния шероховатых поверхностей пластин и не связан с объемным рассеянием.
Анализ измеренных ДФП показал, что для пластин из материала А даже в случае наиболее тонкой пластины пропускание практически диффузное. Из этого был сделан вывод, что оптическая толщина (отношение толщины к длине свободного пробега) всех пластин А велика.
Поведение ДФП пластин из материала В заметно отличалось от ДФП пластин из материала А. ДФП наиболее тонкой пластины имела заметный пик в области, близкой к направлению прямого пропускания, особенно в длинноволновой части спектра. При увеличении толщины пластины пик ДФП заметно
уменьшался. В случае наиболее толстой пластины пик практически отсутствовал. Исходя из этого, был сделан вывод о заметно меньших величинах оптических толщин пластин, изготовленных из материала B.
Кроме этого, анализ также показал не слишком сильную спектральную зависимость измеренных распределений и величин в случае пластин из материала A. При этом форма измеренных ДФО и ДФП достаточно слабо зависит от длины волны освещения.
В случае пластин из материала B ситуация заметно другая. Спектральные зависимости интегральных коэффициентов отражения и пропускания заметно сильнее, чем для материала A. В то же время, спектральная зависимость формы ДФП (за исключением области вблизи направления прямого пропускания) не намного сильнее, чем в случае пластин из материала A.
4 Процедура вычисления параметров объемного рассеяния
В качестве физической модели рассеяния света пластинами была выбрана модель, основанная на феноменологическом Уравнении Переноса Излучения (УПИ) [Radiative Transfer]. Предел применимости такой модели прост: необходимо чтобы характерное расстояние между неоднородностями среды, на которых происходит рассеяние света, было бы много больше длины волны излучения. Как правило, для всех реально применяющихся материалов это условие выполнено с достаточным запасом.
Особо отметим, что УПИ является феноменологическим способом описания распространения света в средах с объемным рассеянием, при этом реальное распределение показателя преломления среды (объемлющая среда с неоднородностями) не используется. Такой подход оперирует со световыми лучами и для описания их распространения в рассеивающей среде использует следующие величины: среднюю длину свободного пробега, поглощение на единицу длины пробега и фазовую функцию (угловое распределение рассеяния). Естественно, все эти величины должны быть заданы в спектральном виде. Одной из самых простых вычислительных реализаций такого подхода является лучевая трассировка методом Монте-Карло.
Вычисление параметров объемного рассеяния может быть сформулировано в виде об-
ратной задачи, т.е. искомые величины или функции находятся из условия совпадения результатов измерений и данных, вычисленных с помощью искомых распределений и величин.
Во многих случаях решение обратной задачи может быть не единственно. Кроме этого, в процессе решения обратной задачи возможно появление различных неустойчивостей. Универсального рецепта преодоления возможных трудностей не существует.
Наиболее часто приближенное решение обратной задачи находится с помощью процедуры минимизации некоторой целевой функции, которая определяется невязками между измеренными и вычисленными величинами. Устойчивость такого подхода и качество процедуры минимизации в значительной степени зависит от адекватности физической модели, удачного выбора параметризации искомых распределений и удачного выбора начального набора варьируемых параметров. Если физическая модель адекватна, то, как правило, число необходимых параметров невелико, а сама процедура минимизации устойчива.
В данном случае применялся именно такой подход. Учитывая, что точность измеренных ДФО/ДФП и интегральных величин невелика и составляет ~5 - 10%, было принято решение ограничиться небольшим числом варьируемых параметров.
Учитывая не очень сильную спектральную зависимость, предварительные вычисления проводились для величин, усредненных по всем измеренным длинам волн (черно-белый режим). В результате этих вычислений удалось для обоих материалов приближенно найти длины свободного пробега и коэффициенты поглощения на единицу длины пробега луча. Методом проб и ошибок, удалось определить аналитический вид адекватной фазовой функции. Было найдено, что фазовая функция подобная фазовой индикатрисе Неупеу-Огееп81ет позволила получить удовлетворительные результаты:
£>(0) =-А-
[1 + g2 - 2g 008(0)]' здесь А — нормировочная постоянная, g — так называемый параметр анизотропии, 0 - угол рассеяния. Использованная фазовая функция при целых и полу-целых значениях у позволяет обращать соответствующую кумулятив-
ную функцию аналитически также как и оригинальная индикатриса Неупеу-Огее^ет (7 = 1.5), что заметно ускоряет трассировочные вычисления.
Кроме этого, было найдено, что углы наклонов поверхностных шероховатостей невелики. Поэтому даже для самых тонких пластин влияние поверхностных рельефов на интегральные коэффициенты отражения и пропускания весьма мало. А это значит, что можно учитывать наличие поверхностного рельефа приближенно, с помощью распределения плотности вероятности наклона нормалей. Такое неявное задание рельефов слабо влияет на интегральные коэффициенты отражения и пропускания.
Плотности вероятности наклонов нормалей поверхностей были приближенно вычислены как средние по всем поверхностям пластин данного материала, используя коэффициент преломления РММА и измеренные ДФО для угла падения 10°.
Далее было установлено, что сильнее всего от искомых параметров зависят формы ДФП наиболее тонких пластин. Слабее - для пластин из материала А и сильнее - для пластин из В. Исходя из этого факта, были найдены величины 7 = 2.5 для материала А и 7 = 10 для материала В. Поиск величин 7 был выполнен полуавтоматически: для нескольких значений 7 выполнялась процедура минимизации для всех 3 пластин данного материала совокупно. При этом целевая функция задавалась как сумма по всем трем пластинам данного материала, где каждое слагаемое было суммой квадратов отклонений отношений вычисленных и измеренных интегральных коэффициентов отражения и пропускания от единицы:
3
F(m, р, g )=
1=1
Эта целевая функция зависит от трех параметров: длины свободного пробега m, коэффициента поглощения р и параметра анизотропии g. Здесь о - угол падения освещающего светового пучка, R(s) и Rexp (о) -вычисленный и измеренный интегральные
коэффициенты отражения, Т (о) и Техр (о) -вычисленный и измеренный интегральный коэффициент пропускания, Ж (о) - вес, с которым учитывался вклад для данного угла падения. Учитывая точность измерений, вклад для угла падения 45° брался с меньшим весом, чем для всех остальных углов падения, а угол падения 60° в целевой функции совсем не учитывался.
Вычисление целевой функции проводилось с помощью лучевой трассировки для всех трех пластин данного материала, используя необходимое число падающих лучей.
Минимизация целевой функции Ф(т, р, ё) осуществлялась с помощью алгоритма нелинейного симплекс - метода или алгоритма Нелдера-Мида. Использованный алгоритм для настоящей целевой функции, по сути, является единственно практически возможным. Действительно, целевая функция вычисляется с помощью лучевой трассировки, что влечет неизбежный стохастический шум. Метод Нелдера-Мида позволяет работать с "зашум-ленными" данными при условии небольшого числа варьируемых переменных, и когда уровень шума не слишком велик. Приемлемый уровень шума достигается трассировкой достаточно большого числа падающих лучей.
Применение более эффективных градиентных методов из-за шума невозможно. Все же другие методы минимизации, не требующие вычисления градиентов, в данном случае малоэффективны.
После завершения процедуры минимизации при фиксированном значении 7 и найденных соответствующих величинах т, р и ё с помощью лучевой трассировки вычислялось угловое распределение прошедшего света для угла падения о = 00 (сечение ДФП при о = 00 ) для наиболее тонких пластин каждого из материалов.
Как уже было упомянуто выше, форма ДФП самых тонких пластин наиболее чувствительна к параметрам объемного рассеяния. Используя вычисленные сечения ДФП
при о = 00 и соответствующие измеренные сечения, усредненные по спектру, были найдены оптимальные величины 7, при которых различия между вычисленными и измеренными распределениями достаточно малы. В процессе предварительных вычислений было установлено, что из-за ограниченной
S
о
S
' r(s) ;
Rexp (о)
t (о) ;
Texp (о)
• W (о) +
• W (о)
2
2
о
точности измерений различия между измеренными и вычисленными сечениями ДФП
при С = 00 достаточно слабо зависят от небольшого изменения величины у в окрестности около найденных значений. Поэтому, как уже было сказано выше, предпочтение было сделано в пользу полу-целых и целых величин с целью ускорить вычисления.
После того, как были определены величины у для обоих материалов А и В вычисление параметров рассеяния для всех измеренных длин волн проводилось аналогично. При фиксированном значении у проводилась минимизация точно такой же целевой функции, в которой интегральные коэффициенты отражения и пропускания брались для заданной длины волны. Естественно, начальная точка (тройка величин Д, р и g) для какой-либо длины волны выбиралась с учетом минимизаций для спектрально усредненных данных и ранее выполненных минимизаций для других длин волн.
Следует сделать ряд замечаний обо всей процедуре поиска параметров объемного рассеяния. На первый взгляд, примененная вычислительная схема выглядит переусложненной и не слишком оптимальной. Более последовательным была бы минимизация по всем четырем параметрам объемного рассеяния: у, т, р и g с целевой функцией, которая бы непосредственно учитывала различия между измеренными и вычисленными ДФО/ДФП.
Однако применение такой вычислительной процедуры сопряжено с очень большими трудностями. Прежде всего, следует еще раз отметить, что решение задачи "извлечения" параметров объемного рассеяния может быть неединственным. В случае решения обратных задач очень большое значение имеет точность экспериментальных данных. Часто при недостаточной точности исходных данных увеличение числа варьируемых переменных практически не улучшает решение задачи минимизации. Более того, возможна ситуация, когда в зависимости от стартовой точки будет найдено несколько решений примерно одинакового качества. При этом найденные наборы варьируемых параметров будут заметно отличаться друг от друга.
В нашем случае относительная точность измерений ДФО/ДФП составляла ~10%, а точность измерения интегральных коэффициентов отражения и пропускания была в 1.5 - 2
раза выше. Вследствие этого было принято решение в первую очередь использовать интегральные коэффициенты, а только во вторую - учитывать форму ДФП. С помощью предварительных вычислений было определено, что для каждой из измеренных длин волн можно ограничиться тремя варьируемыми параметрами, указанными выше. А величину у можно фиксировать для всех пластин, изготовленных из одного и того же материала.
Применение алгоритма Нелдера-Мида в случае зашумленных данных требует большого количества вычислений целевой функции. Более того, количество вычислений целевой функции заметно возрастает при увеличении размерности процедуры минимизации. Кроме этого, вычисление ДФО/ДФП с помощью лучевой трассировки на двумерной угловой сетке адекватного разрешения при достаточно небольшом уровне шума требует существенно большего количества трассируемых лучей, чем вычисление коэффициентов интегрального отражения и пропускания с таким же уровнем шума.
5 Результаты вычислений
На Рис. 1-3 представлены графики для вычисленных параметров объемного рассеяния в спектральном виде. Напомним, что из-за флюоресценции параметры для материала А были определены, начиная с длины волны 430нм.
Результаты вычислений параметров объемного рассеяния для обоих материалов А и В, представленные на графиках выглядят вполне разумно. Порядки всех величин не противоречат физической модели объемного рассеяния. Естественно, найденные параметры были проверены путем сравнения измеренных и вычисленных с помощью найденных параметров ДФО/ДФП. Такие вычисления были проведены для достаточно большого количества длин волн и углов падения освещающего пучка света.
Нужно отметить, что было достигнуто, в целом, хорошее согласие для вычисленных и измеренных распределений. Естественно, согласие несколько хуже в небольшой области вблизи направлений зеркального отражения и прямого пропускания. Кроме этого, заметные отличия наблюдаются для направлений света, рассеянного вблизи поверхности пластин.
Параметр анизотропии g
-•- Материал A * Материал B
Длина волны нм
Рис. 1. Параметр анизотропии g
Длина свободного пробега
-•- Материал A * Материал B
Рис. 2. Длина свободного пробегаm
Коэффициент абсорбции 1/мм
-•- Материал A Материал B
\
Рис. 3. Коэффициент абсорбции р
В обоих случаях эти отклонения вызваны приближенным учетом рассеяния на поверхностных рельефах и не имеют непосредственного отношения к самому объемному рассеянию.
6 Результаты вычислений
В настоящей статье описана не слишком сложная процедура определения параметров объемного рассеяния, исходя из измерений рассеяния света пластин, изготовленных из материалов с таким рассеянием. Предложенная процедура не слишком сложна и дает возможность получить неплохие результаты.
При необходимости качество описанного подхода может быть заметно улучшено, если применить данные, измеренные с большей точностью.
Кроме этого, трудоемкость расчетов может быть значительно снижена, а их быстродействие увеличено, в случае применения более мощных математических и вычислительных методов (например, метод удвоения или параллельные вычисления).
Благодарности
Работа выполнена при частичной поддержке РФФИ, грант № 17-01-00363.
Список литературы
Mie Scattering
URL : https://en.wikipedia.org/wiki/Mie_scattering (дата обращения: 19.02.2018)
Тучин В.В. // Исследование биотканей методом светорассеяния "УФН", 1997, том 167, № 5.
OPTIS
URL: https://www.optis-world.com (дата обращения: 19.02.2018)
Radiative Transfer
URL : https: //en. wikipedia. org/wiki/Radiative_trans fer (дата обращения: 19.02.2018)
Волобой А.Г., Галактионов В.А., Ершов С.В., Летунов А. А., Потемин И.С. // Аппаратно-программный комплекс для измерения свето-рассеивающих свойств поверхностей. «Информационные технологии и вычислительные системы», № 4, 2006, с. 24-39.
Жданов Д.Д., Потемин И.С., Галактионов В.А., Барладян Б.Х., Востряков К.А., Шапиро Л.З. // Спектральная трассировка лучей в задачах построения фотореалистичных изображений. "Программирование", 2011, № 5, с. 13-26.
2 0.3
=т о
=г о
