CC BY
23
ДИФРАКЦИОННАЯ ОПТИКА, ОПТИЧЕСКИЕ ТЕХНОЛОГИИ
Метод расчёта функции эйконала и его применение для синтеза дифракционных оптических элементов для фокусировки в заданную область
Л.Л. Досколович1-2, А.А. Мингазов1-2, Е.В. Бызов1-2, Д.А. Быков1-2, Е.А. Безус12 1ИСОИ РАН - филиал ФНИЦ «Кристаллография и фотоника» РАН, 443001, Россия, г. Самара, ул. Молодогвардейская, д. 151; 2 Самарский национальный исследовательский университет имени академика С.П. Королёва, 443086, Россия, г. Самара, Московское шоссе, д. 34
Аннотация
Предложен метод расчёта функции эйконала светового поля (или фазовой функции), обеспечивающей формирование требуемого распределения освещённости в приближении геометрической оптики. В рамках метода задача расчёта функции эйконала сформулирована в полудискретной форме как задача максимизации вогнутой функции. Для решения задачи максимизации используется градиентный метод, при этом для градиента получены аналитические выражения. С использованием предложенного метода рассчитана функция эйконала, обеспечивающая формирование «разрывного» распределения освещённости в виде изображения гексаграммы. Показано, что использование полученного геометрооптиче-ского решения в качестве начального приближения в итерационном алгоритме Гершберга-Сакстона позволяет рассчитывать дифракционные оптические элементы с квазирегулярной структурой микрорельефа.
Ключевые слова: геометрическая оптика, обратная задача, эйконал, дифракционный оптический элемент, приближение Френеля, алгоритм Гершберга-Сакстона.
Цитирование: Досколович, Л.Л. Метод расчёта функции эйконала и его применение для синтеза дифракционных оптических элементов для фокусировки в заданную область / Л.Л. Досколович, А.А. Мингазов, Е.В. Бызов, Д.А. Быков, Е.А. Безус // Компьютерная оптика. -2022. - Т. 46, № 2. - С. 173-183. - DOI: I0.18287/2412-6179-C0-I029.
Citation: Doskolovich LL, Mingazov AA, Byzov EV, Bykov DA, Bezus EA. Method for calculating the eikonal function and its application to design of diffractive optical elements for optical beam shaping. Computer Optics 2022; 46(2): 173-183. DOI: 10.18287/2412-6179-C0-1029.
Введение
В последние годы задачи дизайна дифракционных оптических элементов (ДОЭ) в рамках скалярной теории дифракции вновь стали предметом интенсивных исследований. Основными причинами интереса к данной области исследований являются компактность, технологичность и эффективность применения ДОЭ при решении широкого класса задач по преобразованию и фокусировке оптического излучения [1 - 5]. Рассматриваемая в данной статье задача расчёта ДОЭ состоит в определении формы «фазового» дифракционного микрорельефа, обеспечивающего формирование требуемого распределения освещённости в некоторой заданной плоскости [6]. Поскольку высота микрорельефа ДОЭ пропорциональна фазовой функции светового поля, формируемого элементом, то задачу расчёта ДОЭ рассматривают как задачу расчёта «фазовой функции ДОЭ» или, более точно, как задачу восстановления фазы по двум известным распределениям освещённости. Для решения данной задачи в приближении Френеля и Фраунгофера предложены различные
итерационные алгоритмы, включающие «классический» алгоритм Гершберга-Сакстона (ГС-алгоритм) и широкий спектр его различных модификаций [6 - 11].
Несмотря на большое количество перспективных применений ДОЭ [1 - 11], включающих задачи преобразования и фокусировки лазерных пучков, формирование сложных изображений в системах дополненной реальности и голографических дисплеях, задачи структурного освещения, оптического захвата, микроскопии и т.п., существующие итерационные методы расчёта ДОЭ имеют ряд важных недостатков. В частности, микрорельеф ДОЭ, рассчитанный итерационным алгоритмом, как правило, имеет сложный и нерегулярный вид, зачастую напоминающий белый шум. Такой рельеф сложен в изготовлении, а погрешности изготовления приводят к неконтролируемому рассеянию света на микрорельефе и формированию спекл-структур, заметно ухудшающих качество формируемого распределения. Кроме того, несмотря на долгую историю развития (ГС-алгоритм был предложен в 1972 году), существующие итерационные алгоритмы всё ещё имеют существенные недостатки, обусловленные проблемами
медленной сходимости и стагнации итерационного процесса в локальных минимумах.
Для «борьбы» с указанными проблемами в последние годы были предложены «гибридные» методы, основанные на совместном использовании итерационных алгоритмов и геометрооптических методов [2, 12, 13]. Это связано с тем, что фазовая функция, рассчитанная в приближении геометрической оптики из условия фокусировки в заданную область, является гладкой или кусочно-гладкой. Использование геомет-рооптической фазовой функции в качестве начального приближения обеспечивает, как правило, быструю сходимость итерационного алгоритма. При этом не возникает проблем стагнации и зацикливания алгоритма в локальных минимумах. Важно отметить, что получаемая в итерационном алгоритме фазовая функция сохраняет квазирегулярный вид, что существенно упрощает изготовление ДОЭ. Более того, существуют модификации итерационных алгоритмов [2], которые позволяют контролировать «гладкость» получаемого решения.
Ключевой составляющей гибридных методов [2, 12, 13] является метод расчёта геометрооптической фазовой функции. Задача расчёта данной функции может быть сведена к решению нелинейного дифференциального уравнения (НДУ) эллиптического типа (уравнения Монжа-Ампера) [12 - 14]. Расчёт геомет-рооптической фазовой функции в работах [2, 12, 13] был основан на численном решении данного НДУ конечно-разностными методами [14 - 16]. Отметим, что эти методы по сути также являются итерационными, поскольку они сводят задачу к решению системы нелинейных уравнений относительно отсчётов фазовой функции на некоторой сетке, которая затем решается итерационным методом. Описанный подход к расчёту геометрооптической фазовой функции накладывает ряд существенных ограничений на класс формируемых распределений освещённости. Действительно, формулировка задачи расчёта геометро-оптической фазовой функции как задачи решения НДУ предполагает, что рассчитываемая фазовая функция является гладкой. Гладкая фазовая функция может реализовать только непрерывное лучевое отображение и поэтому не позволит сформировать распределение освещённости, определённое в несвязной области или в области с негладкими границами.
В данной работе мы предлагаем метод расчёта гео-метрооптической фазовой функции, обеспечивающей формирование требуемого распределения освещённости, в т. ч. в несвязной области или в области со сложными негладкими границами. В рамках метода задача расчёта фазовой функции формулируется в полудискретной форме как задача максимизации вогнутой функции, для решения которой используется градиентный метод. При этом выражения для градиента получены в аналитическом виде. Предлагаемый метод является простым в реализации и основан всего на не-
скольких базовых формулах. В качестве примера рассчитана геометрооптическая фазовая функция, обеспечивающая формирование «разрывного» распределения освещённости в виде изображения гексаграммы. Использование полученного геометрооптического решения в качестве начального приближения в итерационном алгоритме показывает возможность расчёта ДОЭ с квазирегулярной структурой микрорельефа в рамках скалярной теории дифракции.
Статья организована следующим образом. В параграфе 1 приведена постановка задачи расчёта ДОЭ в приближении геометрической оптики. Данная задача формулируется как задача расчёта функции эйконала (или геометрооптической фазовой функции) светового поля, заданной в некоторой исходной плоскости, из условия формирования заданного распределения освещённости в некоторой другой плоскости. В параграфе 2 предложен градиентный метод расчёта функции эйконала. В параграфе 3 рассмотрена задача расчёта фазовой функции ДОЭ в приближении Френеля скалярной теории дифракции. В параграфе 4 представлены результаты расчёта геометрооптической фазовой функции для формирования распределения освещённости в виде изображения гексаграммы, а также результаты оптимизации данной функции с использованием итерационных алгоритмов в рамках скалярной теории дифракции.
1. Постановка задачи в приближении геометрической оптики
Как отмечено выше, задача расчёта ДОЭ обычно формулируется как задача расчёта фазовой функции светового поля, заданной в некоторой исходной плоскости, из условия формирования заданного распределения освещённости в некоторой другой плоскости [6]. Однако в приближении геометрической оптики более правильно работать не с фазовой функцией, а с функцией эйконала, которая равна фазовой функции, делённой на волновое число k=2л/Х, где X - длина волны. Это связано с тем, что в приближении геометрической оптики не учитывается конечность длины волны падающего излучения.
Рассмотрим обратную задачу расчёта функции эйконала более подробно. Пусть в области G в плоскости г = 0 заданы распределение освещённости /о(и), и е G и функция эйконала Ф(и), и е G, где и = (иь и2) - декартовы координаты в плоскости г = 0 (рис. 1).
Функция эйконала Ф(и) определяет направление луча
е(и) = (УФ(и)^1 - (УФ(и))2),
где
ч ( дФ дФ УФ(и) = I-,-
I ди1 ди2
исходящего из точки и плоскости г = 0 в положительном направлении оси г (предполагается, что лучи распространяются в среде с единичным показателем преломления). Здесь и далее верхней стрелкой мы будем обозначать трёхмерные векторы. Лучи позволяют задать соответствие у : О ^ Б между точками исходной плоскости и точками некоторой области Б, расположенной в плоскости г =/> 0. Обозначим х = (XI, х2) декартовы координаты в плоскости г = / Тогда «лучевое» соответствие определяется следующим образом. Точка х = у(и) е Б определяется как точка пересечения плоскости г =f и прямой, проходящей через точку и = (щ, и2) с направляющим вектором е(и). В координатной форме отображение у : О ^ Б имеет вид
х = и + /УФ(и)^1 - (УФ(и))2 .
(1)
Рис. 1. Геометрия задачи расчёта функции эйконала
Распределение освещённости I (х), формируемое в плоскости г =f при функции эйконала Ф(и), определяется из закона сохранения светового потока, который в интегральной форме имеет следующий вид:
{1о(и)аи = {I (х)ах;
(2)
у-1 (В)
где В с Б - произвольное подмножество, а и = у-1(х) -обратное отображение к отображению (1).
Сформулируем обратную задачу расчёта функции эйконала из условия формирования заданного распределения освещённости. Пусть в области О плоскости г = 0 и в области Б плоскости г =f заданы распределения освещённости 10(и), иеО и I (х), хеБ. Требуется рассчитать функцию эйконала Ф(и), иеО, определяющую такое отображение у : О ^ Б, которое удовлетворяет закону сохранения светового потока (2) при заданных распределениях 10(и) и I (х).
Для решения данной обратной задачи удобно представить функцию эйконала в специальном виде. Назовём эйконалом линзы эйконал сходящегося сферического пучка с фокусом в точке х = (х1, х2) плоскости г =/, записанный в плоскости г = 0. Функция эйконала линзы реализует отображение области О в точку х = (х1, х2) и имеет вид
Ф кш (и; х, Т(х)) = Т (х) -р(и, х),
(3)
где
Р(и, х) = д//2 + (х - и)2
- расстояние между точками (и,0) и (х, /). Координаты точки фокусировки х = (х1, х2) в (3) мы рассматриваем как параметры. При этом константа Т(х) соответствует эйконалу в точке фокусировки. Далее мы представим искомую функцию эйконала Ф(и) в виде огибающей семейства функций Ф^ (и; х, Т(х)), обеспечивающих фокусировку в точки плоскости г =/, по параметрам х = (хь х2)еБ [17]. Огибающая семейства функций (3) определяется следующими уравнениями:
Ф (и) = Т (х) — р(и, х),
_д
дХ;
[Т(х) — р(и, х)] = 0, I = 1,2.
(4.1)
(4.2)
Отметим, что функция Т(х) = Ф(и) + р (и, х) совпадает с эйконалом, формируемым в плоскости г =/
Уравнения (4.2) определяют критическую точку функции Т(х) - р (и, х) по переменным х = (х1, х2) при фиксированном значении и = (м1, м2). Далее мы будем рассматривать частный случай огибающей, когда уравнение (4.2) соответствует не просто какой-то критической точке, а глобальному максимуму. В этом случае огибающая примет вид
Ф (и) = тах (Т(х) -р (и, х)). (5)
хеБ
Уравнение (5) позволяет определить лучевое отображение через функцию Т(х) в виде
у (и) = а^тах (Т (х ) -р(и, х)). (6)
хеБ
Несложно получить [17], что для биективного лучевого отображения у : О ^ Б функция Т(х) в (6) определяется следующим соотношением:
Т(х) = тт(Ф(и) + р(и, х)).
иеО
(7)
Уравнения (5), (6) позволяют рассматривать задачу расчёта функции эйконала как задачу расчёта такой функции Т(х), что определяемое ею лучевое отображение (6) удовлетворяет закону сохранения светового потока (2).
Далее мы будем предполагать, что функции Ф(и) и Т(х) являются непрерывными и кусочно-гладкими. Можно показать [17], что формулы (5 - 7) в этом случае также имеют место, разве что отображение у : О ^ Б оказывается многозначным на множестве меры нуль. Заметим, что многозначность отображения на множестве меры нуль не влияет на условие сохранения светового потока (2).
2. Градиентный метод расчёта функции эйконала
В данном параграфе мы рассмотрим градиентный метод для решения задачи расчёта функции эйконала. Пусть формируемое распределение освещённости
I(х), хеВ, аппроксимируется дискретным распределением I, , , = 1, N, заданным в N точках х,еВ (в узлах некоторой сетки, покрывающей область В). В этом случае функция Т(х) заменяется на набор из N чисел Т ,, , = 1, N. Их мы будем варьировать в процессе работы метода. Функция эйконала Ф(и) в этом случае является кусочно-гладкой функцией, состоящей из фрагментов эйконалов линз (3), обеспечивающих фокусировку в точках х,, , = 1, N, плоскости г =/, и, согласно (5), имеет вид
Ф(и) = тах(Т,. -р(и,х,)) = -тт(р(и,х,)-Т,). (8)
/ /
Для фиксированных значений Т,, , = 1, N, можно вычислить значения светового потока, фокусируемые в точках х,, , = 1, N. Согласно (6), (8), в точку х, фокусируется световой поток из области С (х,; Т ...,Тдт) с О, которая определяется следующим условием
c (x,; ¥ь..., ¥ n ) =
ГТ7) (9)
= {и е О | р(и,х,) -Т, < р(и,х,) -Т], V] = 1,N}.
Соответственно, световой поток, фокусирующийся в точку х,, можно вычислить как интеграл
Iest,i 1 est,i (¥b...,¥N) = J 10(u)du .
(10)
c (x, ;¥,,...,¥n )
Если значения Т,, , = 1, N выбраны таким образом, что значения , = 1, N, совпадают с требуемыми значениями I,=I (х,), то задача решена и функция эйконала, обеспечивающая формирование требуемого дискретного распределения, определяется формулой (8). Заметим, что значения Т,, , = 1, N, обеспечивающие выполнение равенств
I, = 1es, Д¥ь..., ¥ n ), , = 1,N ,
(11)
соответствуют стационарной точке следующей функции:
V (¥ь..., ¥ n ) =
= Х J [p(u,x,) - ¥,]lo(u)du + £¥,I,.
i=1 C (x, ;¥,,...,¥„ ) ,=1
(12)
Действительно, вычисление производных функции вида (12), причём при функции р (и, х) произвольного вида, рассмотрено в недавней работе [18]. Согласно [18], указанные производные имеют вид
dV (¥1,..., ¥n )
5¥,
Ii 1 est ,i (^Ь..^ ¥ N ) =
--It - J 10(u)du, i = 1, N.
(13)
C (x, ;¥,,...,¥N )
Более того, несложно показать, что функция К(ТЬ ...,Тл?) является вогнутой [18, 19]. Согласно (13), в критической точке функции К(Ть ...,Тдт) выполнены соотношения (11), причём в силу вогнутости этой
функции эта критическая точка является точкой максимума.
В настоящей работе для отыскания максимума функции V(Ть ...,Тдт) мы будем использовать градиентный метод. В этом случае значения Т,, , = 1, N, рассчитываются итерационно в соответствии с правилом
¥ = ¥ +
А i,n+1 А i,n 1
+д„
Ii - J Io(u)du
V С( ;¥,,„ ....,¥N,n)
i = 1, N,
(14)
где индекс n обозначает номер итерации, а Д„ > 0 -шаг градиентного метода. Важно отметить, что поскольку функция V(¥1, ...,¥N) является вогнутой, то любой локальный максимум будет также и глобальным максимумом. Поэтому для градиентного метода (14) нет проблемы зацикливания в локальном экстремуме, которая характерна для «обычных» градиентных методов.
Градиентный метод (14) имеет простую физическую интерпретацию. Действительно, из формулы (9) нетрудно понять, что при увеличении ¥i (при условии сохранения остальных ¥/ для j Ф i) область C(x,; ¥1, ...,¥n) увеличивается и, соответственно, увеличивается световой поток Iest,i, заданный формулой (10). Поэтому если в точке x, требуемый световой поток Ii > Iestii, то требуется увеличить значение ¥*, и уменьшить его в противном случае. В градиентном методе (14) величина увеличения (уменьшения) ¥* выбирается пропорциональной величине разности (Li- Lesti) между заданным и расчётным значениями светового потока.
Реализация градиентного метода требует вычисления значений Iestj в (10), выражающихся через интегралы по областям C (x,; ¥1, ...,¥N). Для вычисления значений Iest,i можно использовать метод трассировки лучей, который соответствует формуле прямоугольников для вычисления интегралов в (10). В этом случае в области G в плоскости z = 0 генерируется N r лучей, исходящих из узлов некоторой сетки u/, j = 1, Ntr, введённой в области G. Эти лучи «несут» значения светового потока 10/, j = 1, Ntr, пропорциональные значениям освещённости в узлах сетки. Согласно (6), (8), для каждого луча, исходящего из точки uj, находится индекс эйконала линзы
i(j) = argmax(¥( -p(u/,x,)),
i
которая в точке u/ принимает максимальное значение. Этот луч направляется в соответствующую точку xi. Далее в точках xieD плоскости z = f вычисляются световые потоки, равные суммам световых потоков лучей, пришедших в эти точки.
Выше мы рассмотрели градиентный метод расчёта функции эйконала Ф (u) (или геометрооптической
фазовой функции фго (и) = к Ф (и)) в общем непараксиальном случае. В задачах дизайна ДОЭ обычно предполагается выполненным условие параксиального приближения. В этом случае все формулы данного пункта сохраняют свой вид, но вместо расстояния
р(и, х) = ^/2 + (х - и)2
в вышеприведенных формулах следует использовать параксиальную формулу
(х - и)2
р(и, х) Ирраг (и, х) = / +у ' . (15)
2 /
3. Расчёт ДОЭ в приближении скалярной теории дифракции
Рассмотрим задачу формирования светового поля с заданной освещённостью в рамках скалярной теории дифракции. Обозначим ^(и) = ^(и) ехр (;ф (и)) комплексную амплитуду поля в плоскости г = 0, где ф (и), иеО - фазовая функция ДОЭ, Л(и) = -у/Ми) , иеО -амплитуда поля при г = 0. Тогда в приближении Френеля комплексная амплитуда поля в плоскости г =/ будет представляться следующим интегралом:
ж (х ) = —к— ("Г ж, (и )ехр ¡1— (х - и )2 [■ d2u = ^ > XI/ ((oJ ^ 2/У > ] (16)
= Р/ [Ж0 (и)],
где Р/ [•] - оператор распространения света от плоскости г = 0 к плоскости г =/ Отметим, что данный оператор несложно выразить через оператор преобразования Фурье [6].
В приближении Френеля задача расчёта ДОЭ формулируется как задача расчёта фазовой функции ф (и), обеспечивающей формирование в плоскости г =/ светового поля с заданной освещённостью I (х) = (х)|2, хеБ. Отметим, что при вычислении интеграла (16) с использованием приближённого метода стационарной фазы [20] освещённость формируемого поля I(х) = |ж (х)|2 совпадёт с распределением освещённости, рассчитанным в приближении геометрической оптики (в параксиальном приближении) [12, 13]. При этом уравнение на стационарную точку совпадает с параксиальным уравнением луча. Метод стационарной фазы при вычислении интеграла (16) имеет ошибку о (1/к). Соответственно, геометрооптическую фазовую функцию фго(и), рассчитанную в приближении к ^ ж (X ^ 0), можно рассматривать как приближённое решение задачи в рамках скалярной теории дифракции.
Для строгого решения задачи широко используется алгоритм Гершберга-Сакстона, а также его различные модификации [6 - 11]. В ГС-алгоритме фазовая функция ф (и) рассчитывается итерационно, и каждая итерация включает 4 следующих простых шага: 1) Пересчёт поля из плоскости г = 0 в плоскость г =/
wn (x) = Pf [Wo,„_i (u)] = Pf [Ao (u)exp(/'9„_i(u))],
где фп-i (u) - фазовая функция, полученная на предыдущей итерации.
2) Замена амплитуды поля в плоскости z = f на требуемую:
wn (x)^^/7(X)exp (' arg ((x))).
3) Обратный пересчёт поля из плоскости z =f в плоскость z = 0:
wo,n (u ) = P- f [wn (x)] .
4) Замена амплитуды поля в плоскости z = 0 на требуемую:
Wo,n (u) = yl 1o(x) exp ( (u)),
где фп (u) = arg(w0n (u)) - фазовая функция, полученная на текущей n-й итерации. ГС-алгоритм обладает свойством уменьшения (а точнее, неувеличения) ошибки [8]
6n =jj(|wn(x)|-%//СЮ)2dx , (17)
представляющей отличие формируемого и требуемого распределений на n-й итерации. Сходимость ГС-алгоритма существенно зависит от выбора начальной фазы фо (u), используемой на первом шаге алгоритма. Использование в качестве начального приближения геометрооптической фазовой функции фо (u) = фго (u) обеспечивает, как правило, быструю сходимость алгоритма. К сожалению, фазовая функция в процессе итераций не остаётся гладкой (кусочно-гладкой), а приобретает нерегулярную высокочастотную составляющую, которая усложняет процесс изготовления ДОЭ. Для контроля гладкости получаемой фазы можно использовать модификацию ГС-алгоритма со сглаживанием фазы, предложенную в недавней работе [2]. В этом случае фазовая функция фп (u) в плоскости ДОЭ сглаживается на каждой итерации за счёт свёртки с Гауссовой функцией
G(u; а) = (-£)
в соответствии со следующей формулой: фп (u) = ф8 o(u) +
+ mod2^[arg (w o,n (u)) (u)]* G (u; ст).
Параметр ст Гауссовой функции определяет степень сглаживания. В пределе при ст ^ 0 функция G (u; ст) стремится к дельта-функции 5 (u), и тогда фаза (i8) примет вид
фп(u) = ф^(u) + mod2^ [arg(i?o,n(u)) -фg»(u)] . (19)
В силу 2л-периодичности функции ехр(/ф) замена (19) эквивалентна замене фазы фп (u) = arg(w0n(u)), осуществляемой в ГС-алгоритме. При итерационном расчёте ДОЭ в следующем пункте пересчёт фазы в плоскости z = 0 будет осуществляться с использованием формул (18) и (19).
4. Расчётные примеры
В данном параграфе представлены результаты расчёта геометрооптической фазовой функции для формирования «разрывного» распределения освещённости, заданного в несвязной области в виде гексаграммы (рис. 26), а также результаты исследования гибридного подхода к расчёту ДОЭ, в рамках которого геометрооптическое решение используется в качестве начального приближения для итерационного алгоритма в приближении Френеля.
Для получения геометрооптической фазовой функции была рассчитана функция эйконала, заданная в круглой области G в плоскости z = 0, из условия формирования постоянного распределения освещённости в области D в виде гексаграммы (рис. 26) в плоскости z = f. Для расчёта функции эйконала использовался предложенный градиентный метод из параграфа 2 при постоянных распределениях освещённости I0 (u), ueG и I (x), xeD. Расчёт производился при следующих параметрах: радиус области G (апертуры ДОЭ) R = 5 мм, расстояние до плоскости фокусировки f= 100 мм, размер сторон внутренних и внешних треугольников, образующих гексаграмму, d1 = 10 мм и d2 = 15 мм, толщина линий гексаграммы w = 0,45 мм.
В соответствии с методом расчёта требуемое распределение освещённости в виде гексаграммы было аппроксимировано дискретным распределением, содержащим N = 16768 точек, расположенных в узлах квадратной сетки (шаг сетки D = 0,05 мм), покрывающей область гексаграммы. Расчёт эйконала проводился итерационно. В соответствии с описанием градиентного метода на каждой итерации методом трассировки лучей осуществлялся расчёт значений светового потока в точках дискретного распределения, и затем выполнялась коррекция значений ¥j, i = 1, N, по формуле (14). В методе трассировки лучей падающий пучок аппроксимировался набором из N = 4N лучей, заданных на квадратной сетке в плоскости z = 0. Было выполнено около 20000 итераций. Общее время расчёта функции эйконала на персональном компьютере (Core i9-7940X CPU) составило около 2,5 часов.
Рассчитанная функция эйконала Ф^) и формируемое при данной функции распределение освещённости приведены на рис. 2а,в соответственно. Формируемое распределение освещённости было рассчитано методом трассировки лучей и с хорошей точностью совпадает с требуемым распределением: нормированное среднеквадратичное отклонение (СКО) сформированного распределения (рис. 2в) от заданно-
го распределения (рис. 2б) составляет всего 4,5 %. Следует отметить, что полученная функция эйконала является непрерывной и кусочно-гладкой. Изломы (разрывы в производных) у функции эйконала хорошо видны на рис. 2а и обусловлены разрывным видом формируемого распределения. Важно отметить, что при распределениях такого вида методы расчёта геометрооптической фазовой функции в [2, 12, 13], основанные на численном решении уравнения Мон-жа-Ампера, оказываются неприменимыми или являются численно неустойчивыми. Также отметим, что в недавней работе [2], в которой говорится о перспективности использования в качестве начального приближения для итерационного алгоритма геометрооп-тической фазовой функции и представлены результаты расчёта элемента микрооптики, формирующего разрывное распределение в виде изображения робота, расчёт геометрооптической фазовой функции в действительности не проводился. В качестве начального приближения для итерационного алгоритма в [2] использовалась всего лишь квадратичная начальная фаза, обеспечивающая формирование расфокусированного светового пятна, примерно совпадающего с размерами формируемого распределения в виде робота.
Рис. 2. (а) Функция эйконала, рассчитанная из условия формирования постоянного распределения освещённости
в области в виде гексаграммы (б); (в) нормированное распределение освещённости, формируемое в плоскости фокусировки при функции эйконала (а)
На рис. 3а показана геометрооптическая фазовая функция фго (и) = mod2ям [кФ(и)] при длине волны X = 0,55 мкм (к = 2л / X = 11,424 мкм -1). Для удобства визуализации фазовая функция показана по модулю 8л (при М = 4). В силу непрерывности функции эйконала границы зон фазовой функции являются непрерывными, но содержат изломы, обусловленные изло-
мами функции эйконала. На рис. 3б показано распределение освещённости I(х) = (х)|2, формируемое при геометрооптической фазовой функции в приближении Френеля. Распределение освещённости было рассчитано численно с использованием процедуры быстрого преобразования Фурье. Рис. 3б демонстри-
рует достаточно хорошее качество сформированного распределения, хотя и с проявлением дифракционных эффектов (см. врезку). Нормированное СКО сформированного распределения (рис. 3б) от заданного постоянного распределения в виде гексаграммы (рис. 2б) составляет 17,5 %.
У
I — \
1, у1 о
V
ж
V
I __ \
I 1
Рис. 3. (а),(в),(д): Геометрооптическая фазовая функция по модулю 8 Ж (а) при длине волны Х=0,55 мкм и фазовые функции по модулю 8ж, полученные в ГС-алгоритме (в) и в ГС-алгоритме со сглаживанием фазы (д); (б),(г),(е): нормированные распределения освещённости, формируемые при фазовых функциях (а), (в) и (д) соответственно
Рис. 3в, г представляют результат оптимизации геометрооптической фазовой функции в алгоритме Гершберга-Сакстона (с заменой фазы по формуле (19)) и полученное распределение освещённости за 30 итераций (данное число итераций оказалось достаточным для сходимости алгоритма). Распределение освещенности на рис. 3г близко к «идеальному» (рис. 2б). Нормированное СКО сформированного распределения (рис. 3г) от заданного постоянного распределения в виде гексаграммы (рис. 2б) составляет всего 2,1 %. В то же время полученная фазовая функция (рис. 3в) хотя и сохранила регулярный вид, но приобрела высокочастотную составляющую (сравни врезки к рис. 3 а и 3в). В частности, на врезке к рис. 3в видно, что границы зон стали «изрезанными», а внутри зон появились характерные модуляции в виде решёток. Наличие указанных нерегулярностей усложняет процесс изготовления ДОЭ.
На рис. 3д, е показаны результат оптимизации геометрооптической фазовой функции в алгоритме Гершберга-Сакстона со сглаживанием фазы (18) (при ст = 7,4 мкм) и полученное распределение освещённости. Нормированное СКО сформированного распределения на рис. 3е от заданного распределения в виде гексаграммы (рис. 2б) составляет 7,4 %. Это СКО больше, чем в предыдущем случае. В то же время полученная фазовая функция (рис. 3д и врезка) является кусочно-гладкой. Отметим, что за счёт выбора параметра ст можно обеспечить требуемый компромисс между среднеквадратичной ошибкой формируемого распределения и степенью гладкости получаемой фазы. Использованное значение ст = 7,4 мкм было выбрано в соответствии с технологическими возможностями оборудования, которое имеется в распоряжении авторов статьи и которое в последующем планируется использовать для изготовления ДОЭ.
Таким образом, использование геометрооптиче-ской фазовой функции в итерационном алгоритме со сглаживанием (а = 7,4 мкм) позволило получить кусочно-гладкую фазовую функцию, обеспечивающую формирование требуемого разрывного распределения в виде гексаграммы со среднеквадратической ошибкой менее 8 %.
Для демонстрации преимуществ использования точного геометрооптического решения в качестве начального приближения в итерационных алгоритмах был проведён итерационный расчёт фазовой функции при следующей начальной фазе:
фо(и) = -к ^ - £). (20)
Фазовая функция (20) соответствует параксиальной фазовой функции линзы с фокусом
и в приближении геометрической оптики (в параксиальном приближении) обеспечивает в плоскости 2 = / формирование круга радиуса г с постоянной освещенностью. При выбранных параметрах и г = 9 мм гексаграмма оказывается вписанной в данный круг. Отметим, что аналогичная начальная фаза, обеспечивающая формирование светового пятна, совпадающего с размерами формируемого распределения, использовалась в [2] при расчёте элемента, формирующего разрывное изображение в виде робота.
На рис. 4а показана фазовая функция, полученная в алгоритме Гершберга-Сакстона при квадратичной начальной фазе (20) (при г = 9 мм). Данная фазовая функция также обеспечивает формирование заданного распределения освещённости (рис. 4б), но имеет крайне нерегулярный вид. Отметим, что при использовании в ГС-алгоритме в качестве начального приближения случайной или постоянной фазы получаются ещё более нерегулярные фазовые функции, напоминающие по виду реализации некоррелированного случайного процесса типа «белый шум». Важно также отметить, что при начальной фазе (20) или при случайном начальном приближении ГС-алгоритм требует на порядок большего числа итераций. В частности, для получения распределения на рис. 4б было сделано 300 итераций (вместо 30 итераций при геометро-оптической начальной фазе (см. рис. 3в, г)). При этом среднеквадратическая ошибка полученного распределения (рис. 4б) составляет 9,3 %, что даже больше, чем в ГС-алгоритме со сглаживанием фазы при гео-метрооптической начальной фазе (рис. 4г).
В рассмотренном выше примере при длине волны X = 0,55 мкм эффекты дифракции были относительно слабыми. Представляет интерес исследование гибридного подхода в случае более сильных дифракционных эффектов. На рис. 5а показана геометрооптиче-
ская фазовая функция ф^о (и) по модулю 2 л при в 10 раз большей длине волны X=5,5 мкм. В данном случае распределение интенсивности, формируемое в приближении Френеля-Кирхгофа, является существенно более размытым (рис. 5б). Нормированное СКО сформированного распределения (рис. 5б) от заданного постоянного распределения в виде гексаграммы (рис. 2б) составляет 26,4 %.
Рис. 4. Фазовая функция по модулю 8л (при Х = 0,55 мкм), полученная в ГС-алгоритме при квадратичной начальной фазе (20) (а); нормированное распределение освещённости, формируемое в плоскости фокусировки (б)
На рис. 5в, д представлены результаты оптимизации геометрооптической фазовой функции в алгоритме Гершберга-Сакстона и в ГС-алгоритме со сглаживанием фазы при а = 20 мкм. Распределения освещённости, формируемые при данных фазах, показаны на рис. 5г, е соответственно. Нормированные СКО сформированных распределений от заданного постоянного распределения в виде гексаграммы (рис. 2б) составляют 5,6 % и 10,5 %.
В данном случае фазовая функция, полученная в ГС-алгоритме (рис. 5в), имеет существенно более нерегулярный вид, чем исходная геометрооптическая фаза (рис. 5а). В то же время ГС-алгоритм со сглажи-
ванием (18) позволил получить существенно более гладкую фазу (рис. 5д) при относительно небольшом
увеличении среднеквадратической ошибки формируемого распределения.
Рис. 5. (а),(в),(д): Геометрооптическая фазовая функция по модулю 2Ж (а) при длине волны Х = 5,5 мкм и фазовые функции
по модулю 2ж, полученные в ГС-алгоритме (в) и в ГС-алгоритме со сглаживанием фазы (д); (б),(г),(е): нормированные распределения освещённости, формируемые в плоскости фокусировки при фазовых функциях (а), (в) и (д) соответственно
По мнению авторов статьи, предложенный градиентный метод расчёта геометрооптической фазовой функции расширяет возможности гибридных методов синтеза ДОЭ, в особенности в задаче формирования разрывных распределений интенсивности. Получаемые в результате применения итерационного алгоритма распределения фазы оказываются кусочно-гладкими или близкими к ним. Это упрощает технологию изготовления ДОЭ, а также уменьшает негативные эффекты рассеяния света на нерегулярном микрорельефе.
Заключение
В работе предложен геометрооптический метод расчёта функции эйконала (или геометрооптической фазовой функции), обеспечивающей формирование требуемого распределения освещённости, в т.ч. в несвязной области или в области со сложными негладкими границами. В рамках метода задача расчёта функции эйконала сформулирована в полудискретном виде как задача максимизации вогнутой функции. Для решения задачи максимизации был использован градиентный метод, при этом для градиента получены аналитические выражения. Предложенный метод является простым в реализации и основан всего на нескольких базовых формулах.
С использованием предложенного метода была рассчитана геометрооптическая фазовая функция,
обеспечивающая формирование «разрывного» распределения освещённости в виде изображения гексаграммы. Показано, что использование полученного геометрооптического решения в качестве начального приближения в итерационном алгоритме Гершберга-Сакстона позволяет рассчитывать ДОЭ с квазирегулярной структурой микрорельефа. В частности, использование геометрооптической фазовой функции в ГС-алгоритме со сглаживанием (ст = 7,4 мкм) позволило получить кусочно-гладкую фазовую функцию, обеспечивающую формирование требуемого распределения в виде гексаграммы со среднеквадратической ошибкой менее 8 %. При этом полученная фазовая функция имеет значительно более простой вид (с точки зрения технологии изготовления), чем фазовая функция, полученная в ГС-алгоритме при квадратичной или случайной начальной фазе.
Благодарности
Работа выполнена при поддержке гранта Российского научного фонда (проект № 18-19-00326) в части разработки градиентного метода расчёта функции эйконала и Министерства науки и высшего образования РФ в рамках выполнения работ по Государственному заданию ФНИЦ «Кристаллография и фотоника» РАН (соглашение № 007-ГЗ/Ч3363/26) в части дизайна ДОЭ в рамках скалярной теории дифракции.
References
[1] Zhang J, Pegard N, Zhong J, Adesnik H, Waller L. 3D computer-generated holography by non-convex optimization. Optica 2017; 4(10): 1306-1313. DOI: 10.1364/OPTICA.4.001306.
[2] Schmidt S, Thiele S, Toulouse A, Bösel C, Tiess T, Herkommer A, Gross H, Giessen H. Tailored micro-optical freeform holograms for integrated complex beam shaping. Optica 2020; 7(10): 1279-1286. DOI: 10.1364/0PTICA.395177.
[3] Banerji S, Meem M, Majumder A, Vasquez FG, Sensale-Rodriguez B, Menon R. Imaging with flat optics: metalenses or diffractive lenses? Optica 2019; 6(6): 805810. DOI: 10.1364/OPTICA.6.000805.
[4] Banerji S, Sensale-Rodriguez B. A computational design framework for efficient, fabrication error-tolerant, planar THz diffractive optical elements. Sci Rep 2019; 9: 5801. DOI: 10.1038/s41598-019-42243-5.
[5] Pegard NC, Mardinly AR, Oldenburg IA, Sridharan S, Waller L, Adesnik H. Three-dimensional scanless holographic opto-genetics with temporal focusing (3D-SHOT). Nat Commun 2017; 8: 1228. DOI: 10.1038/s41467-017-01031-3.
[6] Soifer VA, Kotlyar VV, Doskolovich LL. Iterative methods for diffractive optical elements computation. London: CRC Press; 1997. ISBN: 978-0-7484-0634-0.
[7] Gerchberg RW, Saxton WO. A practical algorithm for the determination of phase from image and dif-fraction plane pictures. Optik 1972; 35(2): 237-246.
[8] Fienup JR. Phase retrieval algorithms: a comparison. Appl Opt 1982; 21(15): 2758-2769. DOI: 10.1364/AO.21.002758.
[9] Shechtman Y, Eldar YC, Cohen O, Chapman HN, Miao JW, Segev M. Phase retrieval with application to optical imaging. IEEE Signal Process Mag 2015; 32(3): 87-109. DOI: 10.1109/MSP.2014.2352673.
[10] Latychevskaia T. Iterative phase retrieval in coherent dif-fractive imaging: practical issues. Appl Opt 2018; 57(25): 7187-7197. DOI: 10.1364/AO.57.007187.
[11] Ripoll O, Kettunen V, Herzig HP. Review of iterative Fourier transform algorithms for beam shaping applications. Opt Eng 2004; 43(11): 2549-2556. DOI: 10.1117/1.1804543.
[12] Feng Z, Froese BD, Liang R. Composite method for precise freeform optical beam shaping. Appl Opt 2015; 54(31): 9364-9369. DOI: 10.1364/AO.54.009364.
[13] Yang L, Badar I, Hellmann C, Wyrowski F. Light-shaping design by a fourier pair synthesis: the homeomorphic case. Opt Express 2021; 29(3): 3621-3630. DOI: 10.1364/OE.415649.
[14] Bösel C, Gross H. Ray mapping approach for the efficient design of continuous freeform surfaces. Opt Express 2016; 24(13): 14271-14282. DOI: 10.1364/OE.24.014271.
[15] Benamou JD, Froese BD, Oberman AM. Numerical solu-tionof the optimal transportation problem using the Monge-Ampere equation. J Comput Phys 2014; 260: 107126. DOI: 10.1016/j.jcp.2013.12.015.
[16] Prins C, Beltman R, ten Thije Boonkkamp J, Ijzerman W, Tukker T. A least-squares method for optimal transport using the Monge-Ampere equation. SIAM J Sci Comput 2015; 37(6): B937-B961. DOI: 10.1137/140986414.
[17] Doskolovich LL, Mingazov AA, Bykov DA, Andreev ES, Bezus EA. Variational approach to calculation of light field eikonal function for illuminating a prescribed region. Opt Express 2017; 25(22): 26378-26392. DOI: 10.1364/OE.25.026378.
[18] Mingazov AA, Bykov DA, Bezus EA, Doskolovich LL. On the use of the supporting quadric method in the problem of designing double freeform surfaces for collimated beam shaping. Opt Express 2020; 28(15): 22642-22657. DOI: 10.1364/OE.398990.
[19] Merigot Q. A multiscale approach to optimal transport. Comput Graph Forum 2011; 30(5): 1583-1592. DOI: 10.1111/j.1467-8659.2011.02032.x.
[20] Bleistein N, Handelsman RA. Asymptotic expansions of integrals. New York: Dover Publications Inc; 1986. ISBN: 0-486-65082-0.
Сведения об авторах
Досколович Леонид Леонидович в 1989 году с отличием окончил Куйбышевский авиационный институт (КуАИ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королева) по специальности «Прикладная математика». Доктор физико-математических наук (2001 год), профессор, работает заведующим лабораторией дифракционной оптики ИСОИ РАН - филиала ФНИЦ «Кристаллография и фотоника» РАН, профессором кафедры технической кибернетики Самарского университета. Специалист в области дифракционной оптики, лазерных информационных технологий, нанофотоники. E-mail: leonid@ipsiras.ru .
Мингазов Альберт Айдарович, в 2010 году окончил Самарский государственный университет (СамГУ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Математика». В 2012 году окончил магистратуру НИУ ВШЭ по направлению «Математика», в 2015 окончил аспирантуру Санкт-Петербургского отделения математического института имени В. А. Стеклова. Кандидат физико-математических наук (2015), научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН - филиала ФНИЦ «Кристаллография и фотоника» РАН. E-mail: mineazov88@smail.com .
Бызов Егор Владимирович, 1988 года рождения. В 2014 году с отличием окончил обучение в магистратуре Самарского государственного аэрокосмического университета имени академика С.П. Королёва (СГАУ, ныне Самарский национальный исследовательский университет имени академика С.П. Королева) по направлению «Прикладные математика и физика». Область научных интересов: методы расчетов формирующей неизоб-ражающей оптики для светодиодов. E-mail: egor.byzov@smail.com .
Быков Дмитрий Александрович в 2009 году с отличием окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва (СГАУ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Прикладная математика и информатика». Доктор физико-математических наук (2017 г.), старший научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН - филиала ФНИЦ «Кристаллография и фотоника» РАН, доцент кафедры технической кибернетики Самарского университета. Области научных интересов: нанофото-ника, магнитооптика, плазмоника, электромагнитная теория дифракции. E-mail: bvkovd@gmail.com .
Безус Евгений Анатольевич в 2009 году с отличием окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва (СГАУ) по специальности «Прикладная математика и информатика». Кандидат физико-математических наук (2012 г.), научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН - филиала ФНИЦ «Кристаллография и фотоника» РАН, доцент кафедры технической кибернетики Самарского университета. Области научных интересов: нанофото-ника, плазмоника, электромагнитная теория дифракции. E-mail: evgeni.bezus@gmail.com .
ГРНТИ: 26.31.29
Поступила в редакцию 1 октября 2021 г. Окончательный вариант - 10 ноября 2021 г.
Method for calculating the eikonal function and its application to design of diffractive optical elements for optical beam shaping
L.L. Doskolovich 2, A.A. Mingazov 2, E.V. Byzov1-2, D.A. Bykov 2, E.A. Bezus12 1IPSIRAS - Branch of the FSRC "Crystallography and Photonics" RAS, 443001, Samara, Russia, Molodogvardeyskaya 151;
2 Samara National Research University, 443086, Samara, Russia, Moskovskoye Shosse 34
Abstract
We develop a method for calculating the eikonal function (or the phase function) of the light field, ensuring the formation of a prescribed irradiance distribution in the geometrical optics approximation. In the proposed method, the problem being solved is formulated in a semi-discrete form as a problem of the maximization of a concave function. For finding the solution to the latter problem, a gradient method is used, with analytical expressions obtained for the gradient. Using the developed method, we calculate an eikonal function that provides the formation of a "discontinuous" hexagram-shaped irradiance distribution. We demonstrate that the use of the solution obtained in the framework of the geometrical optics as an initial approximation in iterative Fourier transform algorithms allows one to calculate diffractive optical elements having a quasi-regular microrelief.
Keywords: geometrical optics, inverse problem, eikonal, diffractive optical element, Fresnel approximation, Gerchberg-Saxton algorithm.
Citation: Doskolovich LL, Mingazov AA, Byzov EV, Bykov DA, Bezus EA. Method for calculating the eikonal function and its application to design of diffractive optical elements for optical beam shaping. Computer Optics 2022; 46(2): 173-183. DOI: 10.18287/2412-6179-CC>-1029.
Acknowledgements: The development of the gradient method and the design of the geometrical-optics phase function was funded by the Russian Science Foundation under project 18-1900326 and the design of DOEs in the framework of the scalar diffraction theory was funded by the Russian Federation Ministry of Science and Higher Education (State contract with the "Crystallography and Photonics" Research Center of the RAS under agreement 007-GZ/Ch3363/26).
Authors' information
Leonid Leonidovich Doskolovich graduated with honors (1989) from S.P. Korolyov Kuibyshev Aviation Institute (presently, Samara National Research University), majoring in Applied Mathematics. He received his Doctor in Physics & Maths (2001) degree from Samara National Research University. Head of laboratory of the IPSI RAS - Branch of the FSRC "Crystallography and Photonics" RAS, professor at Samara University's Technical Cybernetics sub-department. Current research interests include diffractive optics, laser information technologies, nanophotonics. E-mail: leonid@ipsiras.ru .
Albert Aidarovich Mingazov graduated from Samara State University (presently, Samara National Reseach University) in 2010, majoring in Mathematics. In 2012 he got the master degree in Mathematics at National Research University Higher School of Economics, in 2015 he graduated from the post-graduate course of the St. Peterburg department of Steklov Mathematical Institute of Russian Academy of Sciences. Candidate of Physics and Mathematics (2015). Currently he is a researcher at Diffractive Optics laboratory of the Image Processing Systems Institute RAS -Branch of the FSRC "Crystallography and Photonics" RAS. E-mail: mingazov88@gmail.com .
Egor Vladimirovich Byzov (b. 1988) graduated with honors (2014) from Samara State Aerospace University named after S.P. Korolyov (now - Samara National Research University named after academician S.P. Korolyov), majoring in Applied Mathematics and Physics. Research interests: design methods of nonimaging optics for LEDs. E-mail: egor.byzov@gmail.com .
Dmitry Alexandrovich Bykov graduated with honors (2009) from Samara State Aerospace University (SSAU), majoring in Applied Mathematics and Computer Science. Doctor of Physics and Mathematics (2017). Senior researcher at Samara University and at the Diffractive Optics Laboratory of the Image Processing Systems Institute RAS - Branch of the FSRC "Crystallography and Photonics" RAS. His current research interests include optics of resonant diffractive structures, magneto-optics of nanostructured materials, and electromagnetic diffraction theory. E-mail: bykovd@gmail.com.
Evgeni Anatolievich Bezus graduated with honors (2009) from the Samara State Aerospace University named after academician S.P. Korolyov (SSAU), majoring in Applied Mathematics and Computer Science. Candidate in Physics and Mathematics (2012). Currently he is a researcher at the Diffractive Optics Laboratory of the Image Processing Systems Institute (IPSI RAS - Branch of the FSRC "Crystallography and Photonics RAS") and an associate professor at Technical Cybernetics department of Samara University. His current research interests include nanophotonics, plasmon-ics and electromagnetic diffraction theory. E-mail: evgeni.bezus@smail.com .
Received october 01, 2021. The final version - November 10, 2021.
