CC BY
16
Метод согласованных квадрик для коллимированных пучков
А.А. Мингазов1, Л.Л. Досколович1,2, Д.А. Быков1-2 1ИСОИ РАН - филиал ФНИЦ «Кристаллография и фотоника» РАН, 443001, Россия, г. Самара, ул. Молодогвардейская, д. 151, 2 Самарский национальный исследовательский университет имени академика С.П. Королёва, 443086, Россия, г. Самара, Московское шоссе, д. 34
Аннотация
В статье рассматривается задача расчета преломляющего элемента с двумя поверхностями, формирующего плоский фронт и заданное распределение освещенности. Формулируется метод согласованных квадрик для расчета данного оптического элемента, и показывается, что данный метод совпадает с градиентным методом для некоторого функционала, связанного с задачей перемещения масс Монжа-Канторовича. Это дает возможность адаптивного выбора шага в методе согласованных квадрик. В конце статьи приводится расчетный пример.
Ключевые слова: геометрическая оптика, неизображающая оптика, обратная задача, задача Монжа-Канторовича о перемещении масс.
Цитирование: Мингазов, А.А. Метод согласованных квадрик для коллимированных пучков / А.А. Мингазов, Л. Л. Досколович, Д. А. Быков // Компьютерная оптика. - 2021. -Т. 45, № 1. - С. 29-37. - DOI: 10.18287/2412-6179-CO-783.
Citation: Mingazov AA, Doskolovich LL, Bykov DA. Supporting quadric method for colli-mated beams. Computer Optics 2021; 45(1): 29-37. DOI: 10.18287/2412-6179-CO-783.
Введение
Задача расчета преломляющей или отражающей оптической поверхности из условия формирования заданного распределения освещенности в некоторой области относится к классу обратных задач неизоб-ражающей оптики и является крайне сложной. В приближении геометрической оптики данная задача сводится к решению нелинейного дифференциального уравнения (НДУ) эллиптического типа. Решение этого НДУ является сложной теоретической и вычислительной задачей.
Одним из подходов к данной задаче является метод согласованных квадрик [1 - 3]. Он состоит в следующем. Требуемое распределение освещенности приближается распределением, сосредоточенным в конечном числе точек некоторой сетки. Оптический элемент представляется в виде кусочно-гладкой поверхности, состоящей из фрагментов оптических элементов, каждый из которых фокусирует распределение освещенности в одну точку. В зависимости от задачи эти поверхности могут быть параболоидами вращения, эллипсоидами, гиперболоидами или более сложными поверхностями (например, [3]). Далее параметры фрагментов подбираются итерационно в зависимости от того, насколько отличается количество энергии, приходящее в каждую точку, от требуемого количества.
Мы описываем метод согласованных квадрик для задачи расчета оптического элемента из двух преломляющих поверхностей, формирующих световой пучок с плоским волновым фронтом и заданным распределением освещенности. Далее мы пользуемся
формулировкой данной задачи как задачи перемещения масс Монжа-Канторовича [4], чтобы описать метод согласованных квадрик как градиентный спуск функционала, связанного с функционалом Лагранжа для задачи перемещения масс.
В статьях [1 - 3] в качестве критерия качества используется среднеквадратичное отклонение, которое в процессе работы алгоритма меняется нетривиальным образом. Результат данной статьи позволяет использовать другой критерий качества, который в процессе работы алгоритма согласованных квадрик изменяется монотонно (до тех пор пока значение шага метода не станет слишком большим). Это позволяет адаптивно выбирать шаг, что может ускорять сходимость алгоритма. Хотя в статье рассматривается задача расчета преломляющего элемента, формирующего изображение с плоским фронтом, результаты могут быть легко распространены на любую задачу неизображающей оптики, которая может быть переформулирована как задача перемещения масс Мон-жа-Канторовича.
1. Постановка задачи
Рассмотрим трехмерное пространство Е3 с координатами (х1, х2, г). Пусть в области О плоскости г = 0 задана функция 1(х), определяющая распределение интенсивности коллимированного светового пучка, распространяющегося в направлении к = (0, 0, 1). Геометрия задачи приведена на рис. 1.
Зададим оптический элемент с двумя поверхностями свободной формы. Пусть первая преломляющая поверхность Я1 задается функцией м1(х), х = (х1, Х2) £ О. Тогда
Rl = {(x,щ1(x))|x е G}.
Зафиксируем некоторую плоскость г=I, декартовы координаты на этой плоскости будем обозначать (у\, у). Вторая преломляющая поверхность R2 будет задаваться с помощью функции и2(у) следующим образом
R2 = {(у, г) е Е31г = I - и2(у)}.
Рис. 1. Геометрия задачи
Мы предполагаем, что показатель преломления среды при г < щ(х) равен п0, при щ (х) < г < I - и2(х) показатель преломления равен П1, при I - щ2 (х) < г показатель преломления снова равен по.
Луч, выходящий из точки х е G, преломляясь на двух поверхностях R1, R2, попадает в некоторую точку у области Б в плоскости г = I. Таким образом, определяется лучевое отображение у : G ^ Б. Отображение у формирует в плоскости г = I распределение освещенности, задаваемое функцией Ь(х), которая определяется условием сохранения светового потока:
J I(x)dx = J L(y)dy,
у-1 (В )
где В с Б - произвольное борелевское подмножество.
Задача, которую мы рассматриваем, в классической постановке формулируется следующим образом. Пусть задано распределение интенсивности источника 1(х), х е G, и требуемое распределение освещенности Ь(х) в области Б плоскости г = I. Требуется найти такие дифференцируемые функции (и1( х ), и2 (у )), чтобы выполнялись два условия:
1) Сформированное распределение освещенности имеет плоский фронт, параллельный плоскости г = I.
2) Отображение у : G ^ Б взаимно однозначно и удовлетворяет условию сохранения светового потока:
J I(x)dx = J L(y)dy,
y-1( в) в
(1)
где B с Б — произвольное борелевское подмножество.
Лемма 1. Условие сохранения энергии для отображения у : G ^ Б равносильно любому из условий:
1) I(x) = L(y(x))Jy(x), а Jy (x) - якобиан отображения у, вычисленный в точке x;
2) для любой непрерывной на Б функции h выполнено
J h(y) L(y)dy = J h( у (x)) I (x)dx.
Б G
Proof. Доказано в [5, lemma 4.11] или [6, VI.1.1]. □
В силу леммы 1.1 классическое решение данной проблемы - это решение дифференциального уравнения Монжа-Ампера для функции отображения у специального вида. Мы сформулируем понятие слабого решения для данной задачи, не предполагающее гладкости рефракторов.
2. Функция стоимости
Пусть луч из точки x плоскости z = 0, преломляясь на двух поверхностях R1, R2, приходит в точку y плоскости z = l. Длина оптического пути F между этими точками не зависит от x и y. Число F является параметром оптического элемента, который требуется выбрать. Он неявно связан с толщиной оптического элемента (чем больше F, тем толще оптический элемент). При слишком маленьком параметре F, например, при F < l решения рассматриваемой задачи не существует. Итак, вычисляя длину оптического пути, получим соотношение
(2)
n0u1(x) + n^II y - x II2 +(l - U1(x) - u2(y))2 +
+noU2( y) = F. Обозначим K(x, y) = Ui (x) + u2(y).
Тогда соотношение (2) является уравнением на
K (x, y):
no K (x, y) + n^J | y - x |2 +(l - K (x, y))2 = F. Решая уравнение (2) относительно K, получим
n2l - n0 F n1
(3)
K (x, y) = ■
n2 - no2
n2 - no2
<7(F - nol)2 - (n12 - no2) I y - x I2.
В [7] показано, что оптическому элементу соответствует только решение
K ( х, y) =
ni2/ - n0 F ni
ni - По2
I ni2 - По2 I
<у!(F - По/)2 - (ni2 - По2) н У - х н2.
Итак, имеет место следующая теорема. Теорема 1. Если система из преломляющих поверхностей (м1(х), и2(у)) отображает луч, выходящий из точки х в точку у, то выполнено соотношение
и (х) + и 2 (у) = К (х, у), где
П2/ - п0 Е п1
K ( х, y) =
ni - по I ni - по I
х^(F - По/)2 - (ni2 - По2) Il y - х I
Нам понадобится представление дифференцируемых рефракторов Я1, Я2 в виде огибающих семейств поверхностей специального вида. Пусть лучевое соответствие у : О ^ Б переводит точку х е О в точку у е Б. Для х е О луч, приходящий в точку (х, и1(х)) преломляющей поверхности Я1, далее приходит в точку (у, I - и2 (у)) второй поверхности. Тогда поверхность Я1 может быть представлена в виде огибающей семейства поверхностей Яу, фокусирующих лучи в точки (у, I - и2 (у)) второй поверхности Я2.
Уравнение рефрактора Яу, фокусирующего в точку (у, I - и2 (у)), нетрудно получить, записав условие равенства оптического пути из точки плоскости г = 0 до точки (у, I - и2 (у)) величине Ь - и2 (у). Заметим, что в теореме 1 проведено вычисление именно этого оптического пути. В результате функция и{ (х), задающая рефрактор Яу, удовлетворяет уравнению:
и{ (х) + и2 (у) = К(х, у), то есть
иу (х) = К (х, у) - и2( у).
Итак, функция и1(х) удовлетворяет уравнениям огибающей семейства поверхностей и? (х), то есть выполнено:
ui( х) = K ( х, y) - u2( y),
д —
—[K ( х, y) - u2( y)] = о, i = i, 2.
.dy-
(4)
Второй рефрактор Я2 аналогично можно рассматривать как огибающую семейства поверхностей, фокусирующих при обратном распространении все лучи в одной точке (х, щ(х)).
Щ(у) = К(х,у) -х),
д -
—[К(х, у) - и (х)] = 0,, = 1,2.
дх,
Мы будем рассматривать частный случай огибающей (4), когда второе условие соответствует не просто некоторой экстремальной точке, а глобальному минимуму. В этом случае соотношения примут вид
Ui ( х) = min( K ( х, y) - U2 (y )).
yeD
(5)
Лемма 1. Пусть для системы рефракторов (mi, м2), индуцирующих взаимно однозначное лучевое соответствие, имеет место представление (5). Тогда для функции м2 выполнено
м2 (у) = min(K(х, y) - Ml (х)).
ХЕ.О
Proof. Из (5) сразу следует, что
Mj(х) < K(х, у) - м2(у) Vx e G, у е
Значит, выполнено
м2 (у) < K(х, у) - Mi (х) Vx е G, у е
Отсюда, переходя к минимуму по х, получаем неравенство
м2 (у) < min( K (х, у) - м1 (х)).
Предположим теперь, что существует уо е D, для которого
м2( у о) < min( K (х, уо) - Mi( х)).
xeG
Значит,
M2 (уо) < K(х,уо) - Mi (х) Vx е G, или
м1 (х) + м2 (уо) < K(х, уо) Vx е G.
Последнее равенство означает, что ни для какой точки х е G оптический путь до уо не совпадает с K (х, уо), то есть никакая точка не отображается в уо, что противоречит биективности у. □
3. Слабое решение и метод согласованных квадрик
Мы сформулируем понятие слабого решения, не требующее дифференцируемости функций (м1(х), м2(у)). Это понятие, в частности, определяет, что мы будем понимать под решением в случае дискретного распределения освещенности, задаваемого не непрерывной функцией L, а набором точек {у е D} с заданным в каждой из них световым потоком Ly. В этом случае функция м2 является дискретным набором значений, а м1 является кусочно-гладкой.
Определение 1. K-выпуклой парой (м1(х), м2(у)) будем называть пару непрерывных функций м1 ( х ) на G, м2 (у ) на D таких, что
M1( х) = min(K (х, у) - M2( у));
yeD
м2 (у) = min( K (х, у) - м1 (х)).
xeG
х
Лемма 1. Пусть (u^x), u2(y)) - K-выпуклая пара. Тогда функции u1, u2 липшицевы.
Proof. Аналогично доказательству [8, lemma 3.2]. □
По теореме Радемахера [6, XI.4.2] липшицева функция является почти всюду дифференцируемой. Поэтому, если (ui(x), щ(у)) - K-выпуклая пара, то ui(x) и u2(y) почти всюду дифференцируемы. Значит, они определяют почти всюду дифференцируемые преломляющие поверхности. Поэтому определено (почти всюду) лучевое соответствие у : G ^ D, как отображение преломления на двух поверхностях R1, R2 (см. раздел 1).
Лемма 2. Пусть (u1(x), u2(y)) - K-выпуклая пара, у : G ^ D — лучевое соответствие, определенное почти всюду. Тогда минимум в формуле
ui (x) = min( K (x, y) - u2 (y))
y£D
достигается в точке y = y(x) (там, где отображение у определено). Аналогично для функции u2(y).
Proof. Пусть xo £ G - точка, в которой лучевое соответствие у определено, y0 - точка, в которой достигается минимум. Тогда в точке x0 функция u1(x) касается функции K (x, y0) - u2(y0), которая соответствует поверхности, фокусирующей весь световой пучок в точке (y0, l - u2( y0)). А значит, y(x0) = y0. □
Поскольку, как уже было отмечено, функции u1(x), u2(y) почти всюду дифференцируемы, многозначное отображение
У(x) = {y £ D|u1 (x) + u2 (y) = K(x, y)}
почти всюду однозначно. Следовательно, можно говорить о том, удовлетворяет ли оно условию сохранения энергии или нет. Многозначное отображение у для K-выпуклой пары можно записать также следующим образом
у( x) = argmin(K (x, y) - u2( y)). (6)
yeD
Определение 2. Слабым решением называется K-выпуклая пара функций (u1(x), u2(y))£C (G) x C(D), для которой обобщенное отображение преломления
у (x) = arg min(K (x, y) - u2 (y))
y£D
удовлетворяет условию сохранения энергии: для любого борелевского B с D вполнено
j L(y)dy = j I(x)dx. (7)
B у-1 (B)
Замечание 1. Описанное понятие слабого решения 2 имеет смысл в случае, когда требуемое распределение освещенности L (y) является дискретным, то есть имеет вид
X Ly8 y,
y£D
где 5у - дельта-функция, сосредоточенная в точке у. Но в этом случае функция и2 задана на дискретном множестве Б, а поверхность R2 представляет собой дискретный набор точек. «Преломление» на системе понимается в данном случае в смысле обобщенного отображения преломления (6).
Мы опишем геометрически метод расчета преломляющих поверхностей поверхностей R1, R2, являющийся аналогом метода согласованных квадрик, описанного в статьях [1 - 3] для других задач неизоб-ражающей оптики. Он позволяет приблизиться к слабому решению в случае, когда формируемое распределение освещенности дискретно. В дальнейшем мы докажем, что данный метод является градиентным, а потому сходится при достаточно малом шаге.
Пусть формируемое распределение освещенности представлено дискретным распределением
ХА 8 у,
1=1
у, е Б. Тогда функцию и2(у), как уже было отмечено в замечании 1, заменяет набор из к чисел и2,, где , = 1, к. Их мы будем варьировать в процессе работы алгоритма. Рефрактор Rl представляет собой кусочно-гладкую поверхность, состоящую из фрагментов, фокусирующих исходящий световой пучок в точки (у ,, I—Щ2 ,). Более точно, из определения слабого решения функция и1(х), задающая поверхность Rl предста-вима в виде
Щ1(х) = тт( К (х, уI) - и2,-). (8)
/
Теперь для фиксированных значений и2,, , = 1, к, можно вычислить дискретное распределение, формируемое парой (и1, и2). В точку и 2, фокусируется световой поток с области С(у,) с G, которая определяется следующим условием
С (у,) = {х е 01К (х, у,) - и2, < К (х, у) - и^ ,У/ = 1, к}.
Соответственно, световой поток, фокусирующийся в точку и2 , можно вычислить как интеграл
| I (х)й?х.
С (у,)
Нетрудно понять, что при увеличении и2, (при условии сохранения остальных и2, для ] Ф,) область С( у,) увеличивается, соответственно, увеличивается световой поток.
Нам известно, что требуемый световой поток в точке и2, равен А,. Соответственно, если
А > | I (х)й?х,
С (у,)
требуется увеличить и2 и уменьшить в противном случае. Метод согласованных квадрик предполагает
вычисление следующего приближения иП+' для второй поверхности по предыдущему и2п таким образом
= un + л • (Li - J I(х)ёх), i = i, к,
(9)
C ( y, )
где "л — шаг алгоритма, выбираемый отдельно.
В дальнейшем мы покажем, что метод согласованных квадрик является градиентным спуском некоторого функционала, связанного с задачей перемещения масс.
4. Задача перемещения масс
Задача о нахождении слабого решения оказывается эквивалентна задаче перемещения масс Монжа-Канторовича. Данная теорема аналогична классическим результатам [5, 8, 9], также схожий результат для функции эйконала используется для расчета оптического элемента, формирующего изображение в ближней зоне. Переформулировка рассматриваемой задачи в виде задачи перемещения масс была получена в статье [4].
Точная формулировка приведена в следующей теореме. Мы приводим ее, поскольку доказательство является отправной точкой для дальнейшего рассуждения.
Теорема 1.
Рассмотрим функционал
F(T) = J K(х,Т(х))I(х)ёх,
(1о)
заданный на пространстве отображений Т: О^Б, удовлетворяющих условию сохранения светового потока
I ( х) = L(T ( х)) • Jt ( х).
(11)
Пусть отображение у является точкой минимума данного функционала. Тогда данное отображение является лучевым соответствием для пары рефракторов R.1, Я2.
Proof. Рассмотрим функционал Лагранжа для рассматриваемой вариационной задачи. Возьмем непрерывную на D функцию ц, и пусть
G(T, ц) = J [ K ( х,Т ( х)) I ( х) + ц( х)( L(T ( х)) )
(12)
xJT ( х) -1 ( х))]й?х.
Так как ц - произвольная функция, в силу биек-тивности Т, вместо нее можно записать X (Т (х)), где X - произвольная непрерывная функция на О. Функционал О мы переобозначим через
Н (Т ,Х) = | [К (х,Т (х)) I (х) + Х(Т (х))( Ь(Т (х)) х
О
хЗТ (х) -1 (х))]й?х.
Разбивая на два интеграла и делая замену, получим
н (Т, X) = | [ К (х,Т (х)) - Х(Т (х))]1 (х)ёх +
О (13)
+| Ь( у)Х( у)4у.
Минимум функционала (10) при ограничениях (11) соответствует некоторой критической точке функционала (12). Вариация функционала О легко вычисляется. Условие 5хН = 0, как нетрудно понять из формулы (13), равносильно условию сохранения светового потока. Условие равенства нулю вариации по Т 8Т Н = 0, равносильно
_д_
dy-
[K(х, y) -Х(y)] |y=T(х) = о, i = 1,2.
(14)
Как нетрудно увидеть, это в точности условие огибающей (4), то есть отображение Т задается системой рефракторов, и функция X определяет поверхность второго рефрактора. □
5. Модификация функционала Лагранжа и градиентный метод
В доказательстве теоремы 1 мы видели, что если (Т, X) является критической точкой функционала Лагранжа Н (Т, X), то отображение Т : О ^ Б удовлетворяет уравнению Монжа-Ампера и задается системой рефракторов (и:, и2), при этом X совпадает с и2. То есть в критической точке отображение Т: О ^ Б определяется функцией X условием (14), а именно
ТХ ( хо) = Уо ;
д
dyi
[K(хо,y)-Х(у)]|у=у0 = о, i = 1,2.
(15)
Аналогично тому, как мы ограничивались рассмотрением огибающих специального вида (5), мы будем рассматривать критические точки специального вида, а именно: те, для которых второе соотношение не просто определяет критическую точку, а точку глобального минимума. В этом случае отображение Т можно записать аналогично тому, как описывалось лучевое соответствие ранее (6):
ТX (х) = а^шт(К (х, у) у)).
уеБ
В результате мы можем ограничиться рассмотрением функционала
Н (X) = | [ К (х,Т X (х)) - X(T X (х))]1 (х)ёх +
bJ L( у)Х( y)dy.
(16)
Рассмотрим теперь дискретизацию этого функционала. До этого мы рассматривали только биективные отображения из О в Б. В описанной постановке задачи это являлось разумным, поскольку для распределений освещенности, задаваемых обычными функциями, искомое лучевое соответствие является именно таким. Но для вычислений мы будем аппроксимировать непрерывное распределение интенсивности Ь(у) дискретным распределением
X Ly5 y .
yeD
В этом случае функционал (16) сохраняет смысл, Б - конечное множество точек. Тогда функция X определяется набором чисел Xy, y е Б. Отображение T X в этом случае является отображением области G в дискретное множество точек.
Определение 1. Для фиксированного набора чисел Xy, y е Б, обобщенными клетками Вороного будем называть области C X(y) с G, y е Б, определяемые следующим условием. Точка x е G лежит в клетке C X( yo) тогда и только тогда, когда
K(x,yo)-Xyo <K(x,y)-Xy Vy еБ.
Замечание 1. Клетка C X( y) совпадает с прообразом (T X) -1(y). Соответственно, разбиение на обобщенные клетки Вороного - это разбиение G на прообразы точек y е Б при отображении TX. Поэтому разбиение области G на клетки Вороного для разных весов можно воспринимать как комбинаторно-геометрический способ перебирать отображения вида T X.
Теперь функционал H( X) становится функцией на Rk:
h(X) = X J [K(x, y) -Xy ]I(x)dx + X XyLy, (17)
yеБ CX (y) yеБ
где X - набор весов обобщенных клеток Вороного.
Аналогичный функционал для случая, когда функция стоимости K является квадратичной, используется для расчета оптического элемента в параксиальном приближении в статье [10], связь этого функционала с задачей перемещения масс с квадратичной функцией стоимости была освещена в статьях [11, 12].
Предложение 1. Функция h ( X) является вогнутой.
Proof. Доказательство аналогично доказательству соответствующего утверждения для квадратичной стоимости [12]. Достаточно показать вогнутость первого слагаемого, которое обозначим
Y(X) = X J [K(x, y) -Xy ]I(x)dx =
yеБ C X (y)
= J [ K (x, T X (x)) - X(T X (x))]I (x)dx.
G
Для произвольного отображения T : G ^ Б и любого xеC X (y) выполнено
K (x, y) -X( y) < K (x,T (x)) - X(T (x)),
значит, для любого xе G верно
K(x,TX (x)) - X(TX (x)) < K(x, T(x)) - X(T(x)).
Интегрируя по I(x)dx, получим соотношение
Y (X) <Yt (X),
где
YT (X) = J [ K (x,T (x)) - X(T (x))]I (x)dx.
G
Поскольку для T = T X величины Y ( X) и YT ( X) совпадают, мы можем записать
Y(X) = inf YT (X).
T
Функции YT ( X) являются линейными по X, а значит, вогнутыми. Нетрудно показать, что инфимум вогнутых функций является вогнутой функцией. Соответственно, функция Y является вогнутой, следовательно, h является вогнутой, как сумма двух вогнутых. □
В силу вогнутости функции h она имеет не более одного экстремума, который может быть только ее глобальным максимумом. Вопрос существования этого максимума остается при этом нетривиальным, мы не будем его касаться.
Далее мы собираемся использовать градиентный метод для поиска максимума функции h. В следующей теореме вычисляется градиент функции h( X).
Теорема 1. Частная производная функции h( X) по переменной Xy имеет вид
■^X" = Ly - J I (x)dx.
y CX ( y)
Proof. Для доказательства достаточно вычислить производную функции Y, введенной ранее в доказательстве утверждения (1):
Y(X) = X J [K(x,y)-Xy]I(x)dx.
¡п=Б CX (y)
Для удобства введем обозначение
M (x,X) = min(K (x, y) - Xy).
yеБ
Тогда, пользуясь определением клеток Вороного, можно записать
Y(X) = X J [K(x,y)-Xy]I(x)dx =
yеБ CX (y)
= J M (x,X) I (x)dx.
G
Пусть 5y - вектор (0,_, 0, 5, 0,...,0), имеющий ненулевую компоненту только на y-м месте. Рассмотрим разность
Y(X + 5y) -Y(X) =
= J (M (x, X + 5 y) - M (x,X)) I (x)dx.
G
Предположим для определенности, что 5 < 0 (для 5 >0 можно провести аналогичное рассуждение, но области разбиения будут несколько другие). Тогда клетка Вороного C X (y), вес которой мы варьируем, уменьшается. Значит, имеется цепочка из двух вложений
С^у (у) с СX (у) с О.
Разобьем последний интеграл на три, продолжив равенство:
Т^+бу ) =
= | (М (хД + 5 у) - М (хД)) I (х)ёх +
С ^у ( у )
+ | (М (хД+ 5 у)-M(x,X))I(x)dx +
СX (у)\С"-+5у (у)
+ | (М (хД + 5 у) - М (хД)) I (х)ёх.
О\СX (у)
Рассмотрим отдельно каждый интеграл, начиная с первого. Так как С^ (у) с СX (у), оба минимума достигаются в точке у:
М (хД + 5у) = К (х, у) -Ду +5); М (хД) = К (х, у) -7,у.
Тогда мы получаем
| (М (хД + 5у) - М (хД)) I (х)ёх =
| I (х)ёх.
сХ+5у (у)
= -5-
С"у (у)
Для второго интеграла мы уже не можем сказать, в какой точке достигаются оба минимума. Тем не менее, для любой точки х е СX (у) \ С ^ (у) взвешенное расстояние К (х, у) -Xу до произвольной точки у е Б либо не изменилось, либо изменилось на |5|. Соответственно, минимум этих расстояний изменился не более, чем на |5|. То есть с учетом знаков можно записать
0 < М (х, X + 5 у) - М (хД) < -5.
В результате для второго интеграла получаем оценку
0 < | (М (хД + 5 у) -
С11 (у)\С ^ (у)
-М(хД))I(х)ёх <-5- | I(х)ёх.
С(у)\С"-+5у (у)
Заметим, что при 5 ^ 0 разность ^(у) / С+^у стремится к пустому множеству, а значит, и интеграл
| I(х)ёх
СX (у)\СМу (у)
также стремится к нулю. То есть мы установили, что
Л
lim
5^о
Л
(M( х, Х+5у ) - M ( х, Х))1 (х)сСх
СХ (у)\СХ+5у (у)
= о.
Третий интеграл равен нулю, поскольку оба минимума достигаются в других клетках, а там значение функций под знаком минимума одинаково. Теперь легко вычислим производную:
дТ = У(Х + 5у ) -Т(Х)
дХ " 5^о 5
= - J I (х)Сх.
С ( у )
Зная градиент функции к ( X), можем искать максимум, пользуясь градиентным методом. Более точно, следующее приближение Xn + 1 вычисляется по предыдущему следующим образом:
Xп+1 у (Ьу - | I(х)ёх), (18)
Сxn (у)
где ^ — шаг градиентного метода. В качестве критерия качества имеет смысл рассматривать функцию к ( X). К сожалению, предсказать значение к ( X) в точке максимума невозможно, тем не менее, разумно уменьшать шаг градиентного метода, когда значение к ( X) перестает увеличиваться при операции (18).
Поскольку функция X, как мы видели в доказательстве теоремы 1, является поверхностью второго рефрактора, алгоритм градиентного спуска в точности совпадает с методом согласованных квадрик, описанным ранее (9).
6. Расчетный пример
Рассмотрим расчет оптического элемента, преобразующего пучок с круглым сечением и равномерным распределением
I (х) = Чт, х е G = {(х1, х2 ) | х? + х2 < R2}
nR2
(19)
в плоский пучок с распределением в виде квадрата. Более точно, формируемое распределение определяется формулой
L(у) = -1-, у e D,
w2
(2о)
где D = {(хь х2) | max |{| х1 |, | х2 |} < w}, i = 1,2 - два квадрат со стороной w. Расчет проводился по следующим параметрам: расстояние l = 1о мм; оптическая длина пути F = 11 мм; радиус падающего луча R = 1 мм; размер формируемого квадрата w = 1 мм; индексы преломления материала оптического элемента и окружающей среды n1 = 1,5 и по = 1.
Для вычисления оптического элемента требуемое распределение освещенности L (у) было аппроксимировано дискретным распределением, определенным на квадратной сетке из N = 1оо2 точек. Для вычисления второй поверхности (то есть значений m2, i, i = 1,..., N) используется метод согласованных квадрик (формула 9), описанный в параграфе 5. После вычисления значений m2, ,, i = 1,., N вторая поверхность аппроксимируется, используя метод наименьших квадратов, двумерным 5-сплайном. После вычисления второй поверхности первая поверхность восстанавливается с помощью соотношения (8). Для расчета приведенного примера проделано 2оооо
□
итераций (9). Время работы на персональном компьютере (Intel Core i9-7940X CPU) 4,5 минуты, вос-
становление поверхностей рефрактора заняло около 3 минут.
Рис. 2. Формирование плоского пучка с распределением освещенности в виде квадрата
Полученные поверхности показаны на рис. 1а. На рис. 16 - г приведены нормализованные распределения освещенности, сформированные оптическим элементом в плоскостях z = l = 10 мм, z = 20 мм и z = 30 мм соответственно.
Представленные результаты моделирования демонстрируют формирование прямоугольной рамки с необходимыми размерами и практически равномерной освещенностью. Нормализованные среднеквадратичные отклонения полученных распределений освещенности от равномерного составляют 7%, 10 % и 15 % при z = 10 мм, z = 20 мм и z = 30 мм соответственно. Отметим, что размеры квадрата на разных расстояниях от оптического элемента остаются неизменными, поэтому волновой фронт генерируемого выходного луча является плоским.
Благодарности
Работа выполнена при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию ФНИЦ «Кристаллография и фотоника» РАН в части численной реализации алгоритма расчета, а также грантов РФФИ (№ 18-29-03067, 18-07-00982) в части формулировки метода согласованных квадрик и доказательства совпадения с градиентным методом для соответствующего функционала.
Литература
1. Kochengin, S.A. Computational algorithms for constructing reflectors / S.A. Kochengin, V.I. Oliker // Computing and Visualization in Science. - 2003. - Vol. 6. - P. 15-21. -DOI: 10.1007/s00791-003-0103-2.
2. Досколович, Л.Л. О применении метода согласованных квадрик к расчёту дифракционных оптических элементов / Л.Л. Досколович, М.А. Моисеев, Н.Л. Казанский // Компьютерная оптика. - 2015. - Т. 39, № 3. - С. 339-346. -DOI: 10.18287/0134-2452-2015-39-3-339-346.
3. Андреева, К. В. Метод расчёта оптических элементов с поверхностью свободной формы, работающей по принципу полного внутреннего отражения / К.В. Андреева, М.А. Моисеев, С.В. Кравченко, Л.Л. Досколович // Компьютерная оптика. - 2016. - Т. 40, № 4. - С. 467474. - DOI: 10.18287/2412-6179-2016-40-4-467-474.
4. Rubinstein, J. Intensity control with a free-form lens / J. Rubinstein, G. Wolansky // Journal of the Optical Society of America A. - 2007. - Vol. 24, Issue 2. - P. 463-469. -DOI: 10.1364/JOSAA.24.000463.
5. Glimm, T. Optical design of two-reflector systems, the Monge-Kantorovich mass transfer problem and Fermat's principle / T. Glimm, V.I. Oliker // Indiana University Mathematics Journal. - 2004. - Vol. 53, Issue 5. - P. 12551277. - DOI: 10.1512/iumj.2004.53.2455.
6. Макаров, Б.М. Лекции по вещественному анализу / Б.М. Макаров, А.Н. Подкорытов. - Санкт-Петербург: БХВ-Петербург, 2011. - 688 с.
7. Yadav, N.K Monge-Ampere problems with non-quadratic cost function: application to freeform optics / N.K. Yadav. -Eindhoven: Technische Universiteit Eindhoven, 2018. -152 p.
8. Glimm, T. Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem / T. Glimm, V.I. Oliker // Journal of Mathematical Sciences. - 2003. -Vol. 117. - P. 4096-4108. - DOI: 10.1023/A:1024856201493.
9. Wang, X.-J. On design of a reflector antenna II / X.-J. Wang // Calculus of Variations and Partial Differential Equations. - 2004. - Vol. 20. - P. 329-341. - DOI: 10.1007/s00526-003-0239-4.
10. Schwartzburg, Y. High-contrast computational caustic design / Y. Schwartzburg, R. Testuz, A. Tagliasacchi, M. Pauly // ACM Transactions on Graphics. - 2014. -Vol. 33, Issue 4. - 74. - DOI: 10.1145/2601097.2601200.
11. Aurenhammer, F. Minkowski-type theorems and least-squares clustering / F. Aurenhammer, F. Hoffmann, B. Aronov // Algorithmica. - 1998. - Vol. 20, Issue 1. -P. 61-76. - DOI: 10.1007/PL00009187.
12. Merigot, Q. A multiscale approach to optimal transport / Q. Merigot // Computer Graphics Forum. - 2011. - Vol. 30, Issue 5. - P. 1584-1592. - DOI: 10.1111/j. 1467-8659.2011.02032.x.
Сведения об авторах
Мингазов Альберт Айдарович, в 2010 году окончил Самарский государственный университет (СамГУ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Математика». В 2012 году окончил магистратуру НИУ ВШЭ по направлению «Математика», в 2015
окончил аспирантуру Санкт-Петербургского отделения математического института имени В. А. Стеклова. Кандидат физико-математических наук (2015), научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН - филиала ФНИЦ «Кристаллография и фотоника» РАН. E-mail: mingazov88@,gmail.com .
Досколович Леонид Леонидович, в 1989 году с отличием окончил Куйбышевский авиационный институт (КуАИ, ныне — Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Прикладная математика». Доктор физико-математических наук (2001 год), профессор, главный научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН (ИС-ОИ РАН) — филиала ФНИЦ «Кристаллография и фотоника» РАН, профессор кафедры технической кибернетики Самарского университета и ведущий научный сотрудник научно-исследовательской лаборатории прорывных технологий дистанционного зондирования Земли Самарского университета. Специалист в области дифракционной оптики, лазерных информационных технологий, нанофотоники. E-mail: leonid@ipsiras.ru .
Быков Дмитрий Александрович, в 2009 году с отличием окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва (СГАУ, ныне — Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Прикладная математика и информатика». Доктор физико-математических наук (2017 г.), старший научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН (ИСОИ РАН) - филиала ФНИЦ «Кристаллография и фотоника» РАН и доцент Самарского университета. Области научных интересов: оптика резонансных дифракционных структур, магнитооптика, электромагнитная теория дифракции. E-mail: bykovd@gmail.com .
ГРНТИ: 29.31.01
Поступила в редакцию 18 июля 2020 г. Окончательный вариант - 14 октября 2020 г.
Supporting quadric method for collimated beams
A.A. Mingazov1, L.L. Doskolovich1,2, D.A. Bykov1,2 'IPSIRAS - Branch of the FSRC "Crystallography and Photonics" RAS, 443001, Samara, Russia, Molodogvardeyskaya 151, 2Samara National Research University, 443086, Samara, Russia, Moskovskoye Shosse 34
Abstract
We consider the problem of calculating a refractive element with two surfaces, forming a flat front and a given distribution of illumination. The supporting quadrics method is formulated for calculating a given optical element and it is shown that this method coincides with the gradient method for some functional related to the problem of the Monge-Kantorovich mass transfer problem. This enables adaptive selection of the step in the supporting quadric method. At the end of the article a design example is given .
Keywords: geometric optics, non-imaging optics, inverse problem, Monge-Kantorovich mass transfer problem.
Citation: Mingazov AA, Doskolovich LL, Bykov DA. Supporting quadric method for collimated beams. Computer Optics 2021; 45(1): 29-37. DOI: 10.18287/2412-6179-CO-783.
Acknowledgements: This work was supported by Ministry of Science and Higher Education of the Russian Federation (FSRC "Crystallography and Photonics" RAS) in part of parts of the numerical implementation of the calculation algorithm and Russian Foundation for Basic Research (18-29-03067, 18-07-00982) in part of of formulating the supporting quadric method and proving the coincidence with the gradient method for the corresponding functional.
References
[1] Kochengin SA, Oliker VI. Computational algorithms for constructing reflectors. Comput Vis Sci 2003; 6: 15-21. DOI: 10.1007/s00791-003-0103-2.
[2] Doskolovich LL, Moiseev MA, Kazanskiy NL. On using a supporting quadric method to design diffractive optical elements. Computer Optics 2015; 39(3): 339-346. DOI: 10.18287/0134-2452-2015-39-3-339-346.
[3] Andreeva KV, Moiseev MA, Kravchenko SV, Doskolovich LL. Design of optical elements with TIR freeform surface. Computer Optics 2016; 40(4): 467-474. DOI: 10.18287/2412-6179-2016-40-4-467-474.
[4] Rubinstein J, Wolansky G. Intensity control with a freeform lens. J Opt Soc Am A 2007; 24(2): 463-469. DOI: 10.1364/JOSAA.24.000463.
[5] Glimm T, Oliker VI. Optical design of two-reflector systems, the Monge-Kantorovich mass transfer problem and Fermat's principle. Indiana Univ Math J 2004; 53(5): 1255-1277. DOI:10.1512/iumj.2004.53.2455.
[6] Makarov B, Podkorytov A. Real analysis: Measures, Integrals and applications: Measures, integrals and applications. London: Springer-Verlag; 2013.
[7] Yadav NK. Monge-Ampere problems with non-quadratic cost function: application to freeform optics. Eindhoven: Technische Universiteit Eindhoven; 2018.
[8] Glimm T, Oliker VI. Optical design of single reflector systems and the Monge-Kantorovich mass transfer problem. J Math Sci 2003; 117: 4096-4108. DOI: 10.1023/A:1024856201493.
[9] Wang X-J. On design of a reflector antenna II. Calc Var Partial Differ Equ 2004; 20: 329-341. DOI: 10.1007/s00526-003-0239-4.
[10] Schwartzburg Y, Testuz R, Tagliasacchi A, Pauly M. High-contrast computational caustic design. ACM Trans Graph 2014; 33(4): 74. DOI: 10.1145/2601097.2601200.
[11] Aurenhammer F, Hoffmann F, Aronov B. Minkowski-type theorems and least-squares clustering. Algorithmica 1998; 20(1): 61-76. DOI: 10.1007/PL00009187.
[12] Merigot Q. A multiscale approach to optimal transport. Comput Graph Forum 2011; 30(5): 1584-1592. DOI: 10.1111/j.1467-8659.2011.02032.x.
Authors' information
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 Phisics 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 .
Leonid Leonidovich Doskolovich graduated with honours (1989) from S.P. Korolyov Kuibyshev Aviation Institute (presently, Samara National Reseach University), majoring in Applied Mathematics. He received his Doctor in Physics & Maths (2001) degree from Samara State Aerospace University. Main reseacher of the Image Processing Systems Institute RAS - Branch of the FSRC "Crystallography and Photonics" RAS, professor at Technical Cybernetics depart-
ment of National Research University, senior researcher at the Breakthrough Technologies for Earth's Remote Sensing laboratory at SSAU. His leading research interests include diffractive optics, laser information technologies, nanopho-tonics. E-mail: leonid@smr.ru .
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 the 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 ofnanostructured materials, and electromagnetic diffraction theory. E-mail: bykovd@gmail.com .
Received Jule 18, 2020. The final version - October 14, 2020.
