CC BY
13
ДИФРАКЦИОННАЯ ОПТИКА, ОПТИЧЕСКИЕ ТЕХНОЛОГИИ
Метод опорных квадрик в задачах неизображающей оптики, допускающих переформулировку в виде задачи перемещения масс
А.А. Мингазов1, Л.Л. Досколович ',2, Д.А. Быков 12, Е.В. Бызов1 1ИСОИ РАН - филиал ФНИЦ «Кристаллография и фотоника» РАН, 443001, Россия, г. Самара, ул. Молодогвардейская, д. 151; 2 Самарский национальный исследовательский университет имени академика С.П. Королёва, 443086, Россия, г. Самара, Московское шоссе, д. 34
Аннотация
В статье рассматриваются обратные задачи расчёта оптических элементов для формирования заданных распределений освещённости, которые могут быть переформулированы как задачи перемещения масс Монжа-Канторовича. Для всех задач такого типа мы единообразно формулируем метод опорных квадрик и показываем, что он совпадает с градиентным методом для нахождения максимума некоторой вогнутой функции.
Ключевые слова: неизображающая оптика, геометрическая оптика, обратная задача, задача Монжа-Канторовича, метод опорных квадрик.
Цитирование: Мингазов, А.А. Метод опорных квадрик в задачах неизображающей оптики, допускающих переформулировку в виде задачи перемещения масс / А.А. Мингазов, Л.Л. Досколович, Д.А. Быков, Е.В. Бызов // Компьютерная оптика. - 2022. - Т. 46, № 3. -С. 353-365. - DOI: I0.18287/2412-6179-C0-I055.
Citation: Mingazov AA, Doskolovich LL, Bykov DA, Byzov EV. Support quadric method in non-imaging optics problems that can be reformulated as a mass transfer problem. Computer Optics 2022; 46(3): 353-365. DOI: 10.18287/2412-6179-C0-1055.
Введение
Задача расчёта преломляющей или отражающей оптической поверхности из условия формирования заданного распределения освещённости в некоторой области относится к классу обратных задач неизображающей оптики и является крайне сложной. В приближении геометрической оптики данная задача сводится к решению нелинейного дифференциального уравнения (НДУ) эллиптического типа. Решение этого НДУ является сложной теоретической и вычислительной задачей. Несмотря на то, что для решения данного НДУ были предложены различные конечно-разностные методы [1 - 7], расчёт оптических элементов в рамках такого подхода является сложным и имеет ограничения. В частности, формулировка задачи расчёта оптического элемента в виде НДУ предполагает, что рассчитываемые поверхности оптического элемента являются гладкими. Требование гладкости ограничивает класс распределений освещённости, которые могут быть сформированы оптическим элементом. Например, оптический элемент с гладкими поверхностями не позволяет сформировать распределение освещённости, определённое в несвязной области или в области с негладкими границами.
Одним из наиболее универсальных подходов к данной задаче является метод опорных квадрик [8 -10]. Он состоит в следующем. Требуемое распределение освещённости приближается дискретным рас-
пределением, сосредоточенным в конечном числе точек. Оптический элемент представляется в виде кусочно-гладкой поверхности, состоящей из фрагментов оптических элементов, каждый из которых фокусирует падающий на элемент пучок в одну точку. В зависимости от задачи эти поверхности могут быть параболоидами вращения, эллипсоидами, гиперболоидами или более сложными поверхностями (например, [10]). Далее параметры фрагментов подбираются итерационно в зависимости от того, насколько отличается количество энергии, приходящее в каждую точку, от требуемого значения.
Мы рассматриваем задачи неизображающей оптики, формулируемые в виде задачи перемещения масс Монжа-Канторовича. В первом параграфе данной статьи мы приводим описание всех оптических элементов, для которых известно к настоящему моменту, что задача расчёта элемента из условия формирования заданного распределения освещённости может быть представлена в виде задачи перемещения масс. Для всех этих задач мы единообразно формулируем метод опорных квадрик (МОК), доказывая его совпадение с градиентным методом для нахождения максимума некоторой вогнутой функции. Это переносит результаты статьи [11] на широкий класс задач.
В статьях [8, 9] в качестве критерия качества в МОК используется среднеквадратичное отклонение между формируемым и требуемым распределениями, которое в процессе работы алгоритма меняется не-
тривиальным образом. Результат данной статьи позволяет использовать другой критерий качества для описанного ряда задач, который имеет некоторые преимущества по сравнению со среднеквадратичным отклонением, в частности, в процессе работы алгоритма опорных квадрик изменяется монотонно.
Статья организована следующим образом. В первом параграфе описываются оптические элементы, для которых задача расчёта из условия формирования заданного распределения освещённости может быть переформулирована как задача перемещения масс Монжа-Канторовича, а также их представления в виде огибающих специального вида. Во втором параграфе мы единообразно формулируем задачу расчёта оптического элемента из условия формирования заданного распределения освещённости. В третьем параграфе мы формулируем итерационный метод опорных квадрик для расчёта описанных оптических элементов. Затем в четвёртом параграфе показываем, что МОК является градиентным методом для поиска максимума некоторой выпуклой функции. В последнем параграфе мы приводим два примера расчёта преломляющих элементов, формирующих заданные распределения освещённости в дальней зоне. На первом из них мы иллюстрируем преимущества рассмотрения в качестве критерия сходимости метода опорных квадрик предложенной выпуклой функции перед средне-квадратическим отклонением. В статье использована единая сквозная нумерация определений, теорем, лемм и утверждений.
1. Описания рассматриваемых задач и представление поверхностей оптических элементов в виде огибающих
В данном параграфе мы кратко опишем несколько оптических элементов, для которых задачи формирования заданных распределений освещённости допускают переформулировку в виде задачи перемещения масс (ЗПМ) Монжа-Канторовича с некоторой функцией стоимости К. Данные оптические элементы рассмотрены в отдельных работах [11 - 23]. Для каждого элемента мы введём понятие прямой ф и т.н. двойственной функции у и лучевого отображения Т. В следующем параграфе мы, используя введённые обозначения, приводим общую для всех оптических элементов формулировку проблемы в терминах функций у, T и функций I, L, представляющих интенсивность источника и требуемое распределение освещённости (или интенсивности).
1.1. Преломляющий (отражающий) элемент, формирующий заданное распределение интенсивности при точечном источнике
Пусть имеется точечный источник света, расположенный в начале координат в трёхмерном пространстве с координатами (х1, х2, х3). Рассмотрим еди-
ничную сферу § с центром в начале координат. Положение точки
х = (хьх2Л/ГХГХГ) е §
на полусфере xз > 0 будем задавать координатами x = (х1, x2). Источник света будем характеризовать функцией интенсивности I' (х1, x2), x е G, равной световому потоку, излучаемому в единицу телесного угла. Обозначим ^ элемент телесного угла, тогда световой поток через область ю с § вычисляется как
Ф
(ю) = {I'(х^ст = |т
ю ю V 1
I'(х)
Х1 X2
dx1dx2.
Далее нам будем удобно характеризовать излучение источника функцией
Iм=- I'(х)
, x е G.
(1)
X22
Пусть излучение от точечного источника падает на преломляющую поверхность W. Будем считать, что показатель преломления среды под поверхностью W (в области расположения источника) равен щ, а в области над поверхностью - п2.
Рассмотрим некоторый луч, вышедший из начала координат (из источника) с направлением х . Данный луч преломляется на поверхности W и приобретает направление единичного вектора
У = (УьУ2,\Д-У2 -У22 ) е § .
Это определяет «лучевое отображение» T : G ^ §. Отображение T задаёт на сфере § функцию распределения интенсивности L' (у) для преломлённого пучка. Функция L'(у) однозначно определяется условием сохранения светового потока, который в интегральной форме имеет вид:
Ф, (Т-1(ю)) = Ф4Я (ю) Уюс §,
(2)
где
Ф,
,(ю) = |Г( y)dст = ¡1=
ю ю V 1
_ г L'(y1,У2)
-У12 - У22
Ф -
световой поток преломлённого пучка через область ю с §. Далее мы будем характеризовать распределение интенсивности преломлённого пучка функцией
L( у) =
L'( У)
л/1 - У12 - у2 '
(3)
Преломляющую поверхность W удобно задать функцией ф(х) = - 1п р (х), где р (х) - расстояние от начала координат до поверхности в направлении
= (хь Х2 , \Д - х12 - Х2 ) .
Рассматриваемая в данном пункте задача состоит в расчёте функции ф (х), задающей преломляющую поверхность, из условия формирования заданного распределения интенсивности Ь (у), у е Б.
Для решения данной задачи искомую поверхность удобно представить в виде огибающей семейства эллипсоидов фр-у (х) = 1п(1 - п(х, у» - 1п р(у), где п = П2 / «1, р (у) - фокальный параметр [14], по параметрам у = (у:,у2) е Б. Эллипсоид ф р,у (х) преобразует падающий сферический пучок от источника в коллимированный пучок с направлением у = (у\, у2). Уравнение огибающей семейства эллипсоидов имеет вид [14]:
ф( х) = К( х, у) -у( у),
д —
—[К( х, у)-у(у)] = 0, г = 1,2,
дУг
(4)
где у (у) = 1пр (у) - так называемый двойственный параметр, а функция К (х, у) имеет вид
К( х, у) = 1п(1 - «<х, у».
(5)
В [12 - 15] показано, что задача расчёта преломляющей поверхности из условия формирования заданного распределения интенсивности эквивалента ЗПМ Монжа-Канторовича с функцией стоимости К (х, у), заданной уравнением (5). Формулировка ЗПМ приведена в теореме 4 следующего параграфа.
Формально подставив в (5) п = -1, мы получим верные формулы для случая отражающей оптической поверхности. Подробно данная задача рассмотрена в [12 - 15].
1.2. Преломляющий (отражающий) элемент, формирующий заданное распределение интенсивности при освещающем пучке с плоским волновым фронтом
Пусть в трёхмерном пространстве с координатами (х1, х2, х3) в направлении к = (0, 0, 1) распространяется световой пучок с плоским волновым фронтом, имеющий в плоскости П, задаваемой уравнением г = 0, распределение освещённости 1(х), х = (х1, х2) е О. Данный пучок падает на преломляющую поверхность Ж, задаваемую уравнением г = ф (х). Будем считать, что в области г < ф (х) показатель преломления среды равен «1, а в области г > ф (х) равен «2. Луч, вышедший из точки х = (х1, х2) плоскости г = 0, преломляется на поверхности Ж и приобретает направление единичного вектора
У = (у1,у2,л/1-у2 - У2 ) е §.
Это определяет лучевое отображение Т : О ^ §. Данное отображение задаёт на сфере § функцию распределения интенсивности Ь'(у) для преломлённого пучка. Аналогично предыдущему случаю распределение интенсивности преломлённого пучка будем за-
давать функцией Ь (у), определённой формулой (3). Функция Ь ( у) определяется условием сохранения светового потока
Фп(Т-1(ю)) = ФШ(ю) V®
(6)
где
ф ш (ю) = |ь( у)ау -
световой поток преломленного пучка через область на сфере §, а
Фп (Т-1(ю))= | I(х)ах -
Т->0)
световой поток падающего пучка через соответствующую область плоскости г = 0.
Рассматриваемая в данном пункте задача состоит в расчёте преломляющей поверхности Ж из условия формирования заданного распределения Ь (у), у е Б. Аналогично предыдущему случаю представим поверхность г = ф (х) в виде огибающей семейства поверхностей, преобразующих падающий пучок в пучок с плоским волновым фронтом и заданным направлением у = (у1, у2) е Б. В данной задаче такими поверхностями являются плоскости [16] вида
г = К(х, у) -у( у),
где
К( х, у) =
пх1 У1
пх2 у2
1 - Пу] 1 - У12 - у| 1 - 1 - У2 -
у2
(7)
(8)
где п = п2 / щ. Двойственный параметр у (у) в данной задаче имеет геометрический смысл смещения плоскости (7) относительно начала координат. В результате, так же как и в предыдущем случае, преломляющая поверхность г = ф (х) может быть представлена в виде огибающей семейства плоскостей, задаваемой уравнениями (4).
Данная задача также эквивалента ЗПМ Монжа-Канторовича с функцией стоимости (8). Формулировка ЗПМ приведена в теореме 4 следующего параграфа. Заменой координат можно добиться того, что функция стоимости станет квадратичной. Подробнее данная задача рассмотрена в [15, 16].
Формально подставив п = -1 в выражение (8) для функции стоимости К, можно получить формулы для случая отражающего оптического элемента, формирующего заданное распределение интенсивности.
1.3. Элемент из двух отражающих поверхностей, формирующий заданное распределение освещённости и плоский волновой фронт при освещающем пучке с плоским волновым фронтом
Рассмотрим трёхмерное пространство с координатами (хь х2, г). Пусть в области О плоскости г = 0 задана функция I(х), х = (хь х2) е О, определяющая рас-
пределение освещённости коллимированного светового пучка, распространяющегося в направлении к = (0, 0, 1). Пучок падает на оптический элемент с двумя отражающими поверхностями. Пусть первая поверхность элемента Я1 задаётся функцией г = ф (х), х е G. Зафиксируем некоторую «выходную» плоскость г = I, в которой будет рассматриваться преобразованный элементом пучок, и обозначим у = ( у1, у2) декартовы координаты на этой плоскости. Вторую отражающую поверхность Я2 зададим относительно выходной плоскости с помощью функции у (у) следующим образом г = I - у( у).
Мы предполагаем, что луч, выходящий из точки х е G плоскости г = 0 с направлением к = (0, 0, 1), последовательно отражается от двух поверхностей Я1, Я2 и затем попадает в некоторую точку у области Б в выходной плоскости г = I, сохраняя при этом направление к = (0, 0, 1). Это определяет лучевое отображение Т : G ^ Б. В результате на плоскости г = I формируется некоторое распределение освещённости L (у), определяемое условием сохранения светового потока, которое в данной задаче принимает вид:
J I(x)dx = Jl(y)dy, VracD.
(9)
T->0)
Задача состоит в расчёте функций ф (х) и у ( у), задающих поверхности элемента, из условия формирования выходного пучка, имеющего в плоскости г = I заданное распределение освещённости L (у) и плоский волновой фронт. В данной задаче мы будем рассматривать первую поверхность Я1 как огибающую семейства поверхностей, фокусирующих в точки второй поверхности Я2 [15, 17]. В этом случае поверхность семейства должна отражать падающий колли-мированный пучок в некоторую точку с координатами (у, I - у (у)) на второй поверхности. Нетрудно понять, что такими поверхностями являются параболоиды вращения с фокусами в точках второй поверхности, задаваемые уравнением
Фy,-и(y)(x) = Д(x,y) -у(y),
где
Д(x,y) = 62 -12- II Х-У I'2 ,
2( L -1 )
(10)
где Q - оптическая длина пути лучей между плоскостями г = 0 и г = I [17]. Отметим, что оптическая длина пути лучей между плоскостями г = 0 и г = I является постоянной в силу того, что волновые фронты падающего и выходного являются плоскими. Уравнения огибающей аналогично записываются в виде (4). Данная задача также эквивалента ЗПМ с функцией стоимости К [17]. Формулировка ЗПМ приведена в следующем параграфе. Подробное описание данной задачи дано в [15, 17].
1.4. Элемент из двух преломляющих поверхностей, формирующий заданное распределение освещённости и плоский волновой фронт при освещающем пучке с плоским волновым фронтом
Рассмотрим трёхмерное пространство с координатами (хь х2, г). Пусть в области G плоскости г = 0 задана функция I (х), х = (хь х2) е G, определяющая распределение освещённости коллимированного светового пучка, распространяющегося в направлении к = (0, 0, 1). Пучок падает на оптический элемент с двумя преломляющими поверхностями. Первая поверхность Яь задаётся функцией г = ф(х), х е G. Как и в предыдущей задаче, вторую поверхность Я2 зададим относительно некоторой выходной плоскости г = I через функцию у(у) в виде г = I- у (у), где У = (Уь, У2) - декартовы координаты на выходной плоскости. Будем считать, что показатель преломления среды при г < и (х) равен щ, при ф(х) < г < I- у (х) равен п2, а при I - у (х) < г снова равен щ.
Мы предполагаем, что луч, выходящий из точки х е G с направлением к = (0, 0, 1), последовательно преломляется на двух поверхностях Яь, Я2 и затем попадает в некоторую точку у области Б в выходной плоскости г = I, сохраняя направление к = (0, 0, 1). Это определяет лучевое соответствие Т : G ^ Б. В результате на плоскости г = I формируется некоторое распределение освещённости L (у) выходного пучка, определяемое условием сохранения светового потока (9).
Задача состоит в расчёте функций ф (х) и у ( у), задающих поверхности элемента, из условия формирования выходного пучка, имеющего в плоскости г = I заданное распределение освещённости L (у) и плоский волновой фронт.
Как и для случая двух отражающих поверхностей, мы будем рассматривать первую поверхность Я1 как огибающую семейства поверхностей, фокусирующих в точки второй поверхности Я2 [11, 19]. В этом случае имеют место уравнения огибающей (4) с функцией
) n2l - niQ Д( x, y) = ■ ^
П2
ni - n2
| П22 - П12
(11)
(Q - nil)2 - (П22 - П12) I y - x I
где Q - оптическая длина пути лучей между плоскостями z = 0 и z = l. Данная задача также эквивалента ЗПМ Монжа-Канторовича с функцией стоимости (11) [11, 15, 18]. Формулировка ЗПМ приведена в следующем параграфе. Подробное описание данной задачи приведено в [11, 15, 18, 19].
1.5. Элемент из двух отражающих поверхностей, формирующий заданное распределение освещённости и плоский фронт при точечном источнике
Рассмотрим трёхмерное пространство с координатами (x1, x2, z). Пусть в начале координат расположен точечный источник света. Обозначим через § сферу единичного радиуса с центром в источнике. Интен-
сивность источника, как и в случае первого элемента, описанного в пункте 1.1, будем характеризовать функцией I(хь х2), х = (х^ х2)еО, которая связана с фактической интенсивностью I' (х) соотношением (1).
Сферический пучок от точечного источника падает на оптический элемент с двумя отражающими поверхностями. Пусть первая отражающая поверхность Я1 задаётся функцией и (х), х е О:
Я1 = (и; (х) х | х е О},
где
х (х, ■ 1 х2 х2).
Зафиксируем некоторую «выходную» плоскость г = I, в которой будет рассматриваться преобразованный элементом пучок, и обозначим у = (у1, у2) декартовы координаты на этой плоскости. Вторую поверхность Я2 зададим функцией и2 (у) следующим образом
Я2=((у,г) е Е3 |г = 1 - ^(у)}.
Мы предполагаем, что луч, исходящий из источника в направлении х, последовательно отражается от поверхностей Я\, Я2 и затем попадает в точку у некоторой области Б на плоскости г = I, имея при этом направление к = (0, 0, 1). Это определяет лучевое соответствие Т : О ^ Б. Оно определяет в области Б плоскости г = I распределение освещённости Ь (у) выходного пучка, определяемое условием сохранения светового потока
I I(х)йх = \Ь(у)ду, V® с Б.
(12)
Т->0)
Задача состоит в расчёте отражающих поверхностей Я1, Я2, задаваемых функциями и1, и2, из условия формирования выходного пучка, имеющего в плоскости г = I заданное распределение освещённости Ь (у) и плоский волновой фронт.
Для получения соотношения вида (4) в данной задаче поверхности Я1, Я2 следует задавать другими функциями ф (х) и у (у), выражающимися через и1 (х), и2 (у), достаточно сложным образом [20, 21]. Пусть Q - оптическая длина пути лучей от источника до выходной плоскости. Отметим, что данная оптическая длина пути является постоянной, поскольку волновой фронт выходного пучка является плоским. В рассматриваемой задаче прямая и двойственная функции вводятся следующим образом [21]:
ф( х) = 1п
1
1
2и,(х)(1 -<х, к») 2^ -I) ( .. ! А
У(У) = 1п
и2( У)
1
(13)
Q2-I2-У УII2 2(Q-I)
В терминах функций ф и у снова можно записать соотношение (4), рассматривая ф как огибающую се-
мейства функций ф,у<у) (х) = К (х, у) - у (у), соответствующих поверхностям, фокусирующим падающий пучок в точки второй поверхности. В данной задаче функция стоимости имеет вид [20, 21]
К(х, у) =
(
= 1п
1
Q -<х, Л»
4(Q-1)2 2(Q-1)(Q2 -12-1у ||2)(1 -<х,к»)
(14)
где л = (У, I). Условие огибающей снова может быть записано в виде (4). Данная задача эквивалента [20, 21] ЗПМ Монжа-Канторовича с функцией стоимости (14). Формулировка ЗПМ приведена в следующем параграфе. Подробнее данная задача рассмотрена в [20, 21].
1.6. Элемент из двух отражающих поверхностей, формирующий заданное распределение освещённости и сферический фронт при точечном источнике
Рассмотрим трёхмерное пространство с координатами (х1, х2, г). Пусть в начале координат расположен точечный источник света. Обозначим через §1 сферу единичного радиуса с центром в источнике. Интенсивность источника света, как и для элемента в пункте 1.1, будем характеризовать функцией !(х1, х2), х = (х1, х2) е О, заданной соотношением (1).
Сферический пучок от точечного источника падает на оптический элемент с двумя отражающими поверхностями. Пусть первая поверхность Я1 задаётся функцией и1 (х), х е О:
Я1 = (и;(х) х | х е О},
где
(х1, х2 , х х2) .
Рассмотрим также точку Т с координатами (0, 0, Р) и единичную сферу §2 с центром в этой точке. Координаты на §2 обозначим у = (у1, У2), эти координаты определяют точку Т + у сферы §2, где
У = (У\,У2,д/1 У12 -У22).
Вторую отражающую поверхность зададим функцией и2 ( у) следующим образом
Я2=(Т + щ(у)у |у е Б}.
Мы предполагаем, что луч, исходящий из источника в направлении х, последовательно отражается от поверхностей Я\, Я2 и затем приходит в точку Т, пересекая сферу §2 в некоторой точке с координатами (У1, У2). Это определяет лучевое отображение Т : О ^ Б. Данное отображение задаёт на сфере §2 функцию распределения интенсивности Ь'( у) для преобразованного элементом пучка. Как и для элемента в пункте 1.1, формируемое на сфере §2 распределение интенсивности мы будем задавать функцией
Ь (у), уеБ, которая связана с фактическим распределением интенсивности Ь' (у) формулой (3).
Задача состоит в расчёте отражающих поверхностей Я^, Я2, задаваемых функциями щ, щ, из условия формирования выходного пучка, имеющего заданное распределение интенсивности Ь (у) и сферический волновой фронт с центром в точке Т.
Как и в предыдущем случае, для получения соотношения вида (4) в данной задаче поверхности Я\, Я2 следует задавать другими функциями ф (х) и у (у), выражающимися через щ, щ достаточно сложным образом [15].
Введём прямую и двойственную параметризацию системы отражающих поверхностей следующим образом [15]:
ф( х) = 1п
у(у) = 1п
1
1
2^(х)(© -F<к,х» Q2 -F2 1 1
2щ(уШ + F<к,у» Q2 -F2
(15)
где Q - оптическая длина пути лучей от источника до точки Т.
Представляя одну из отражающих поверхностей в виде огибающей семейства поверхностей, фокусирующих в точки другой поверхности, можно получить представление (4) с функцией стоимости [15]
К( х, у) = 1п
1
(Ь2 -F2)2 1 -<х, у)
2(Q - F<к, х))(© + F<к, у))(©2 - F2)
(16)
Данная задача эквивалента ЗПМ Монжа-Канторовича с функцией стоимости (16) [15]. Формулировки ЗПМ приведена в следующем параграфе. Подробно данная задача рассмотрена в [15].
1.7. Функция эйконала на поверхности, обеспечивающая формирование заданного распределения освещённости на другой поверхности
Пусть в трёхмерном пространстве с координатами (х1, Х2,2) расположены две поверхности: Щ\, заданная уравнением 2 =/(х), и Щ2, заданная уравнением 2 = g(х). На поверхности Щ1 заданы функция эйконала ф (х), х = (х1, х2) е О и распределение освещённости I'(х), хеО. Световой поток в области ю с Щвычисляется как
Фщ» = ]/'(х)№,
где - элемент площади поверхности. Нам будет удобно характеризовать распределение освещённости функцией
I (х) = I'(х^1+/х2+/2
где
Л=/ , у = 1,2.
ох-
В этом случае световой поток ФЩ1 вычисляется следующим образом
фщ (ю) = ]/(х)дх.
Функция эйконала ф (х), х = (х1, х2) е О определяет направления лучей, исходящих из точек первой поверхности. Будем считать, что поверхности расположены в среде с единичным показателем преломления. В этом случае направления лучей, распространяющихся в положительном направлении оси 2, задаются единичным вектором
p(х) = (Мх),p2(),
который определяется из следующей системы уравнений [23]:
1т = Pl+/х, -V 1 - Pl2 - P22; 1т- = P2 + /х2 V1 - Pl2 - P22.
ох2
(17)
Функция эйконала задаёт лучевое отображение Т между точками поверхностей Щ1, Щ2. Данное отображение определяется следующим образом. Образом точки хеО является точка у = (у1,у2) е Б, соответствующая пересечению поверхности щ2 и прямой, выходящей из точки (х,/(х)) поверхности Щ1, с единичным направляющим вектором p = ^^,^1- p12 - ), определяемым из системы уравнений (17).
Отображение Т задаёт на поверхности щ2 некоторое распределение освещённости. Аналогично поверхности щ1 будем характеризовать распределение освещённости функцией Ь (у), у е Б, которая связана с фактическим распределением освещённости Ь (у) следующим соотношением:
Ь( у) = Ь'( + .
Тогда результирующий световой поток для ю с Б вычисляется
Фщ2(ю) = |Ь( у )Ау,
и функция Ь (у) определяется условием сохранения светового потока
| I(х)6х = \Ь(у)Ау, Уюс Б.
(18)
Т-1(ю)
Задача состоит в расчёте функции эйконала ф (х), х е О, из условия формирования заданного распределения освещённости Ь (у), у е Б, на второй поверхности щ2.
Аналогично предыдущим задачам, мы представим функцию эйконала ф (х) в виде огибающей семейства функций эйконалов сферических пучков с фокусами в точках (у, g (у)) второй поверхности Ж2. Данные функции имеют вид [22, 23]
фу.у(у)(х) = Дх, у) _у (у), (19)
где
Д (х,У) = ^Ц х_у II2 +(/(х) _g(у))2. (20)
Огибающая семейства функций (19), (20) по параметрам у = (у1,у2) е Б определяется соотношением (4). При этом функция - у (у) совпадает с формируемым на поверхности Ж2 эйконалом.
Данная задача эквивалента задаче перемещения масс Монжа-Канторовича с функцией стоимости (20) [23]. Формулировка задачи перемещения масс приведена в следующем параграфе. Подробнее данная задача рассмотрена в [22, 23].
2. Формулировка обратной задачи
Сформулируем более явно роль представления поверхности оптического элемента в виде огибающей. Будем говорить, что произвольное лучевое отображение Т : О ^ Б задаётся соответствующим оптическим элементом, если существует некоторая функция у (у), для которой выполнено:
д —
—[Д(х,у) _у(у)] |у=Т(х)= 0, г = 1,2. (21)
ду,-
Прямой и двойственный параметры оптического элемента при этом связаны соотношением:
ф( х) + у (Т (х)) = Д( х, Т (х)). (22)
Определение 1. Будем называть отображение Т : О ^ Б интегрируемым, если оно удовлетворяет условию (21).
Определение 2. Будем называть отображение Т : О ^ Б удовлетворяющим условию сохранения светового потока, если для любого борелевского подмножества ю с Б выполнено
| I(х)6х = |цу)ёу. (23)
Т_ (ю) ю
Условие сохранения светового потока может быть переформулировано следующим образом.
Лемма 3. Условие сохранения энергии для обратимого отображения Т : О ^ Б равносильно любому из условий:
1) для любого ю с Б выполнено
I I (х)йх = |цу)6.у;
Т_ (ю) ю
2) для любой функции И е С(Б) выполнено
jh( y) L( y)dy = jh(T (X)) I (x)dx;
D G
3) I(x)=L (T(x)) • Jt(x), где Jt(x) - якобиан отображения T в точке x.
Доказательство. Доказано в [17, lemma 4.11] или [24, VI.1.1].
Задачи расчёта всех рассмотренных в предыдущем параграфе оптических элементов можно сформулировать в одном общем виде следующим образом. Для заданных распределений I(x) и L (у) требуется найти лучевое отображение T : G ^ D, которое удовлетворяет условию сохранения светового потока (23) и условию интегрируемости (21), а также восстановить прямой и двойственный параметры оптического элемента.
Данная задача может быть сформулирована как задача перемещения масс Монжа-Канторовича.
Теорема 4. Рассмотрим функционал
Т (T ) = jK( x,T (x)) I (x)dx, (24)
G
заданный на пространстве отображений T : G ^ D, удовлетворяющих условию сохранения светового потока (23).
Пусть отображение To является точкой минимума данного функционала. Тогда оно является интегрируемым отображением.
Доказательство. Аналогично доказательству [11, Theorem 1].
3. Метод опорных квадрик
Мы будем рассматривать решения специального вида, для которых условие интегрируемости (21) задаёт не просто критическую точку функции K(x, у) - у (у) по переменной у, а точку минимума. Геометрически это соответствует тому, что поверхность оптического элемента, задаваемая функцией ф, расположена по одну сторону от всех касающихся её поверхностей фу,у (y)(x) = К (x, у) - у (у). В этом случае условия (21), (22) можно записать в виде
ф( x) = min (К( x, у)-у(у)). (25)
yeD
При этом экстремум достигается в точке у = T (x), то есть соответствующее отображение T = T у : G ^ D имеет вид
T у (x) = argmin(K( x, у)-у(у)). (26)
yeD
При такой формулировке возможно рассмотрение дискретного варианта задачи, когда формируемое распределение освещённости представлено обобщённой функцией
EL 5 у,,
i=1
у, е Б. В этом случае функция у представляет собой набор чисел, а функция ф является непрерывной и кусочно-гладкой. Далее мы будем рассматривать такую постановку задачи, которая соответствует аппроксимации непрерывного распределения Ь ( у) дискретным, заданным в конечном числе точек {у, }к=1.
Определение 5. Пусть Б = {у, }к=1, у- = у (у,). Обобщёнными взвешенными клетками Вороного с весами у = (уь ... ук) будем называть области Су (у,) равные прообразам (Т у)-1( у,). Из представления (26) следует, что клетки Су (у,) могут быть определены следующим условием: точка х е О лежит в клетке С у(у,) тогда и только тогда, когда
Д (х, у-) -у, < Д(х, у;) - у] V; Ф,'.
Обозначим световой поток, фокусирующийся в точку у, через
О (у)= | I (х)йх. (27)
Су (у,)
Теперь условие сохранения светового потока (23) может быть записано так:
О1(у) = Ьь
■ ... (28)
Оп (у) = Ьп.
Система (28) является системой нелинейных уравнений относительно значений у = (у1, ... ук). Перепишем её следующим образом:
у1 = у1 +е-(А -О1(у)),
■ ... (29)
у п = у п +е- (Ьп - Оп (у)),
где е > 0 - некоторая константа. Рассматривая метод простой итерации для решения системы (29), мы получим следующий итерационный алгоритм.
у п+1 = уп + е • (Ь - О' (у)), ' = ГТп. (30)
Представленный алгоритм имеет понятный физический смысл. Пусть требуемый световой поток в точке у, равен Ь,, а текущий сформированный поток О, (у) для функции у меньше Ь, то есть Ь,> О, (у). Тогда нужно увеличить размер клетки Су (у,), а значит, нужно увеличить вес у,-. Если Ь, < О, (у), то размер клетки Су (у,) нужно уменьшить за счёт уменьшения веса у,. Алгоритм (30) увеличивает или уменьшает у, на величину, пропорциональную разнице требуемого и формируемого световых потоков.
В следующем параграфе мы покажем, что метод (30) является градиентным методом для нахождения максимума функции
^у) = £ | [Д(х, у)-у,- V (х)йх + £у Ь. (31)
'=1 Су (у,) ,-=1
В силу вогнутости функции h (у) она имеет не более одного экстремума, который может быть только её глобальным максимумом. Совпадение метода (30) с градиентным гарантирует его сходимость при достаточно малом значении шага е, а также даёт возможность адаптивного выбора шага метода.
Функция h (у) тесно связана с функционалом Ла-гранжа для задачи перемещения масс Монжа-Канторовича. В данной статье мы ограничимся лишь вычислением её градиента и доказательством вогнутости. Происхождение данной функции описано в статье [11].
4. Вогнутость и вычисление градиента
В этом параграфе мы покажем, что функция (31) является вогнутой, а также вычислим её частные производные.
Утверждение 6. Функция h (у) является вогнутой.
Доказательство. Доказательство полностью аналогично доказательству соответствующего утверждения для квадратичной стоимости [25]. Мы приводим доказательство для полноты изложения.
Для этого сначала покажем вогнутость первого слагаемого, которое обозначим
Я (у ) = Ё | [Д( х, у,)-у,- У (х)йх =
,-=1 С у (у,)
= \[К(х, Ту (х)) - у (Ту (х))У (х) йх.
О
Для произвольного отображения Т: О ^ Б и любого х е Су (у,) выполнено
Д(х, у,) - у,' < Д(х,Т (х)) - у(Т (х)), значит, для любого х е О верно
Д(х,Т у (х)) -у(Ту (х)) < Д(х,Т (х)) - у(Т (х)).
Интегрируя по I (х)йх, получим соотношение
Я (у) < Ят (у),
где
Ят (у) = |[Д(х, Т(х)) -у(Т (х)Ж (х)йх.
О
Поскольку для Т = Ту величины Я (у) и ЯТ (у) совпадают, мы можем записать
Я(у) = ЯТ (у).
Функции ЯТ (у) являются линейными по у, а значит, вогнутыми. Нетрудно показать, что инфимум вогнутых функций является вогнутой функцией. Соответственно, функция Я является вогнутой, следовательно, h (у) также является вогнутой, как сумма двух вогнутых функций.
Теорема 7. Частная производная функции h (у) по переменной у, имеет вид
-д— = Ь _ | I (х)дх.
ду' СУ (у,)
Доказательство. Для доказательства достаточно вычислить производную функции Я (у), введённой ранее в доказательстве утверждения 6. Для удобства введём обозначение:
М (х, у) = шш(Д( х, у,) _ у,).
г
Тогда, пользуясь определением клеток Вороного, можно записать
п
Я(у) = X I [Д х, у,-) _у, ]1(х)^х = |м (х, у)1(х)ах.
'=1 Су (у,) О
Пусть 5, - вектор (0, ..., 0, 5, 0, ..., 0), имеющий ненулевую компоненту только на ,-м месте. Рассмотрим разность
Я(у + 5,) _Я(у) = |(М (х, у + 5,) _М (х, у)) I (х)йх.
О
Предположим для определённости, что 5 > 0 (случай 5 < 0 рассматривается аналогично). Тогда клетка Вороного Су (у,), вес которой мы варьируем, увеличивается, то есть Су (у,) с Су+5- (у,) с О. Разобьём последний интеграл на три:
Я(у + 5,) _Я(у) =
= | (М(х,у + 5,) _М(х,у))!(х^ +
Су (у,)
+ I (М (х, у + 5,) _ М (х, у)) I (х)^ +
Су+5, (у, )\Су (у,)
+ I (М (х, у + 5,-) _ М (х, у)) I (х)^.
О\Су+5, (у,)
Теперь рассмотрим отдельно каждый интеграл, начиная с первого. Так как Су (у,) с Су+5- (у,), то оба минимума достигаются в точке у,, и поэтому
I (М (х, у + 5,-) _ М (х, у))I (х)^ =
)
5- I I(х)йх.
Су (у,)
Су (у,)
Рассмотрим второй интеграл. Для любой точки х е С у+5- (у,) \ Су (у,) взвешенное расстояние Д (х, у,) -у, до произвольной точки у, е Б либо не изменилось, либо изменилось на 5. Соответственно, минимум этих расстояний изменился не более, чем на 5. То есть
| М (х, у + 5,-) _М (х, у) | <5.
В результате для второго интеграла получаем оценку
I (М (х, у + 5,,) _ М (х, у))! (х)^
Су+5, (у, )\Су (у,)
<5- I I (х)йх.
Су+5, (у, )\Су (у,)
Заметим, что при 5 ^ 0 разность Су+5- (у,) \ Су (у,) стремится к пустому множеству, а значит, второй интеграл имеет порядок малости о (5).
Третий интеграл равен нулю, поскольку оба минимума достигаются в других клетках, а там значение функций под знаком минимума одинаково. Переходя к пределу при 5 ^ 0, мы получаем утверждение теоремы.
5. Расчётные примеры
Теоретические положения данной статьи мы проиллюстрируем на примере расчёта преломляющего оптического элемента, формирующего заданное распределение интенсивности при точечном источнике. Данный оптический элемент описан в параграфе 1.1. Будут рассмотрены два расчётных примера. Первый из них имеет небольшую расчётную сложность, на нём хорошо видно преимущество критерия сходимости И (у) (31) по сравнению со среднеквадратичным отклонением. Во втором примере демонстрируется возможность рассчитывать с помощью данного метода оптические элементы, формирующие сложные распределения освещённости.
5.1. Рефракционная поверхность, формирующая постоянное распределение освещённости на контуре квадрата
В качестве первого модельного примера рассмотрим расчёт преломляющей поверхности Ж, формирующей при точечном источнике излучения постоянное распределение освещённости на кривой в виде контура квадрата с размером стороны ё в удалённой плоскости г =/ Преломляющая поверхность разделяет среды с показателями преломления П1 = 1,49 (по-лиметилметакрилат) и п2 = 1. Понятие «удалённая плоскость» предполагает, что размерами преломляющей поверхности по сравнению с расстоянием до плоскости г =/ можно пренебречь. Данную задачу можно переформулировать [16] как задачу расчёта преломляющей поверхности из условия формирования заданного распределения интенсивности (см. подпараграф 1.1). Пусть (ть т2) - декартовы координаты на плоскости г = / а требуемое распределение задаётся функцией Е (т). Точка на этой плоскости однозначно соответствует единичному вектору
8 (т) =
(ть т2, /)
который мы отождествляем с точкой на сфере. Первые две координаты вектора 8 (т) соответствуют введённым в подпараграфе 1.1 координатам
(у1, у2) =
(ть т2)
^ т 2 + т2 + /2
Функция Ь (у (т)), с помощью которой мы характеризуем формируемое распределение интенсивности, связана с Е (т) следующим образом [16]:
L(y(m)) = E(m)
f2 + m2 + m| f
В качестве источника излучения будем использовать точечный ламбертовский источник, излучающий в телесный угол, соответствующий конусу с вершиной в начале координат и углом при вершине 26о. В этом случае интенсивность ламбертовского источника имеет вид
I(x1,x2) = -<J 1 -Xj2 - x|, x = (xbx2) e G ,
где
G =
j(xb x2 ^xj2
+ xf < sin2
■ }•
Расчёт поверхности производился при следующих параметрах: угол при вершине «конуса излучения» 260 = 90°, расстояние от источника до удалённой плоскости / = 1000 мм, размер стороны квадрата й=Дап (60°) = 1376 мм.
В соответствии с описанием метода опорных квадрик (МОК) в параграфе 3, требуемое распределение освещённости в виде контура квадрата было аппроксимировано дискретным распределением, содержащим N = 240 точек {у,}240. Точки были расположены равномерно на контуре с угловым шагом в 1°, то есть все требуемые распределения освещённости Ь, одинаковы. В этом случае преломляющая поверхность Щ представляет собой кусочно-гладкую поверхность из фрагментов N =240 квадрик (эллипсоидов вращения), фокусирующих весь световой поток в точки у,, задающихся уравнением в полярных координатах
р(х) - P
1 - n
{e (x), У i}
Поверхность W тогда задается в полярных координатах функцией
p(x) = max ,
1 - n
{e (x), У,)'
(32)
После применения 1п с учётом обозначений из параграфа 1 данное соотношение эквивалентно формуле
(25) параграфа 3. Поверхность зависит только от набора фокальных параметров p^. Расчёт фокальных параметров эллипсоидов проводился итерационно градиентным методом (30)
у П+1 = уП
+ B-(L -G, (у)),
где у, = lnp, а Gi (у) - световой поток, фокусирующийся оптическим элементом в точку у. Это соответствует такому преобразованию значений p,.
pn+1 = pnes'(Li-G" (у®
При этом на каждой итерации методом трассировки лучей [11] осуществлялся расчёт значений светового потока Gi (у) в точках дискретного распределения, и затем выполнялась коррекция значений у,, i = 1, N. При расчёте было выполнено 1600 итераций (этого оказалось достаточно для сходимости метода). Общее время расчёта преломляющей поверхности на персональном компьютере (Core i9-7940X CPU) составило около 3 минут.
Рассчитанная преломляющая поверхность W приведена на рис. 1а, где чёрной точкой отмечено положение источника излучения. В данном примере расстояние от источника до вершины поверхности составляет 3,1 мм. На рис. 16 приведено формируемое распределение освещённости в плоскости z =f= 1000 мм. Данное распределение было рассчитано методом трассировки лучей [11] и с хорошей точностью совпадает с требуемым распределением. нормированное среднеквадратичное отклонение (СКО) сформированного распределения от заданного постоянного распределения на контуре квадрата составляет всего 2,5 %. Следует отметить, что полученная рефракционная поверхность является непрерывной и кусочно-гладкой. Изломы (разрывы в производных) хорошо видны на рис. 1а и обусловлены тем, что формируемый контур квадрата не является гладкой кривой. Для проверки правильности выполненных расчётов было проведено моделирование работы рассчитанной поверхности в программе для светотехнических расчётов TracePro [26]. Результаты моделирования приведены на рис. 1е и также показывают формирование равномерно освещённого контура.
6) в) I
Рис. 1. Рефракционная поверхность. рассчитанная из условия формирования постоянного распределения освещенности
на контуре квадрата (а); нормированное распределение освещенности. формируемое в плоскости 2 =/= 1000 мм (б); результат моделирования работы преломляющей поверхности в программе для светотехнических расчетов TracePro (в)
"орГ" ю-3
ю-4
ю-5
Ю-6
Ю-7
Ю-8
ККМ8Е(Ье!1ГЬ) *10~2
100
ю-'
ю-2 _ Итерации,
О
4
10 12 14 16
а) 0 0,25 0,50 0,75 1,00 б) Рис. 2. График функции к (у) (сверху) и график нормированного СКО между расчетным и требуемым распределением в виде контура квадрата (снизу) на отрезке (1 - Г)у1 + Гу2, Г е [0,1] (а); графики сходимости МОК при расчете поверхности на рис. 1а: зависимость ошибки \ к (у) - кор\(сверху) и СКО-критерий (снизу) от номера итерации (б)
Как показано в параграфе 4, использованный для расчёта преломляющей поверхности МОК (30) соответствует градиентному методу максимизации вогнутой функции к (у) (31). Для иллюстрации вогнутости функции к (у) мы выбрали два случайных вектора у1 и у2, представляющих параметры эллипсоидов в уравнении для преломляющей поверхности (35), и построили функцию к (у) на отрезке (1 - Г) у1 + Г у2, Г е [0, 1] соединяющем эти точки. Полученная функция к (1 - Г) у1 + Г у2 показана на верхнем графике на рис. 2а и является вогнутой. Для сравнения на нижнем графике на рис. 2а мы показали график нормированной среднеквадратической ошибки (СКО) между требуемыми и расчётными значениями освещённости на контуре квадрата при параметрах эллипсоидов, лежащих на том же самом отрезке (1 - Г) у: + Г у2. На графике СКО хорошо видно наличие нескольких экстремумов. Это показывает, что СКО-критерий, широко используемый для контроля сходимости в МОК [1 - 3], является менее подходящим по сравнению с функцией к (у). Это также иллюстрируется графиками сходимости МОК на рис. 2б, полученными при расчёте преломляющей поверхности на рис. 1а. На верхнем графике рис. 2б приведена величина корГ- к (у), где корГ - значение функции, полученное в конце расчёта (на последней итерации). Из рис. 2б видно монотонное уменьшение величины корГ- к (у). Для сравнения в нижней части рис. 2б показан график СКО-критерия. Из сравнения графиков на рис. 2б видно, что СКО, по сравнению с корГ- к(у), является более плохим критерием для контроля сходимости итерационного процесса. В частности, можно видеть, что СКО-критерий после примерно тысячи итераций имеет зашумлённый и осциллирующий вид.
5.2. Рефракционная поверхность, формирующая распределение освёщенности в виде символа «Триглав»
В качестве более сложного примера, иллюстрирующего возможности предложенного МОК в зада-
чах формирования сложных распределений освещённости, была рассчитана преломляющая поверхность, формирующая постоянное распределение освещённости в области в виде символа «Триглав» (см. врезку к рис. 3в). Как и в примере 1, в качестве источника излучения был использован точечный ламбертовский источник, излучающий в телесный угол, соответствующий конусу с углом при вершине 26о = 90°. При этом угловой размер формируемого распределения (угловой размер квадрата, описанного вокруг символа «Триглав») составляет 100°.
При расчёте требуемое распределение освещённости в виде символа «Триглав» было аппроксимировано дискретным распределением, содержащим N = 20890 точек. Точки были расположены в узлах квадратной угловой сетки (шаг сетки 0,26°), покрывающей область символа. Расчёт преломляющей поверхности (параметров квадрик) проводился градиентным методом (30). Как и в примере 1, на каждой итерации методом трассировки лучей [11] осуществлялся расчёт значений светового потока в точках дискретного распределения и затем выполнялась коррекция значений у,-,, = 1,N по формуле (30). При расчёте было выполнено 25000 итераций. Общее время расчёта преломляющей поверхности составило около 3,5 часов.
Рассчитанная преломляющая поверхность приведена на рис. 3а, где чёрной точкой отмечено положение источника излучения (расстояние от источника до вершины поверхности составляет 3,0 мм). На рис. 3б приведено формируемое распределение освещённости в удалённой плоскости г =/= 1000 мм. Данное распределение было рассчитано методом трассировки лучей [11] при числе лучей Nг = 106. СКО сформированного распределения от заданного постоянного распределения в области символа «Триглав» составляет 3,6 %. Можно видеть, что рассчитанная рефракционная поверхность является непрерывной и кусочно-гладкой. Изломы у рассчитанной поверхно-
сти обусловлены «разрывным» видом формируемого распределения, заданного в многосвязной области. Важно отметить, что при распределениях такого вида широко используемые методы расчёта оптических поверхностей [1 - 7], основанные на численном решении нелинейного дифференциального уравнения эллиптического типа, оказываются неприменимыми
или являются численно неустойчивыми. Для проверки правильности выполненных расчётов, было проведено моделирование работы рассчитанной поверхности в программе ТгасеРго. Результаты моделирования приведены на рис. 3в и показывают формирование примерно постоянного распределения освещённости в области символа.
1,0
0,8
0,6
б) Е^нв^а
Рис. 3. Рефракционная поверхность. рассчитанная из условия формирования постоянного распределения освещенности в области в виде символа «Триглав» (а); нормированное распределение освещенности. формируемое в плоскости 2 =/= 1000 мм (б); результат моделирования работы преломляющей поверхности в программе для светотехнических расчетов TracePro (в). На врезке показано требуемое распределение
Благодарности
Работа выполнена при поддержке Министерства науки и высшего образования в рамках выполнения работ по Государственному заданию ФНИЦ «Кристаллография и фотоника» РАН в части численной реализации алгоритма расчёта, а также гранта РНФ18-19-00326 в части доказательства совпадения метода опорных квадрик с градиентным методом для соответствующего функционала.
References
[1] Chang S, Wu R, Li A, Zheng Z. Design beam shapers with double freeform surfaces to form a desired wavefront with prescribed illumination pattern by solving a Monge-Ampere type equation. J Opt 2016; 18(12): 125602. DOI: 10.1088/2040-8978/18/12/125602.
[2] Feng Z, Froese BD, Huang C-Y, Ma D, Liang R. Creating unconventional geometric beams with large depth of field using double freeform-surface optics. Appl Opt 2015; 54(20): 6277-6281. DOI: 10.1364/A0.54.006277.
[3] Feng Z, Froese BD, Liang R, Cheng D, Wang Y. Simplified freeform optics design for complicated laser beam shaping. Appl Opt 2017; 56(33): 9308-9314. DOI: 10.1364/A0.56.009308.
[4] Bösel C, Worku NG, Gross H. Ray-mapping approach in double freeform surface design for collimated beam shaping beyond the paraxial approximation. Appl Opt 2017; 56(13): 3679-3688. DOI: 10.1364/AO.56.003679.
[5] Bösel C, Gross, H. Double freeform illumination design for prescribed wavefronts and irradiances. J Opt Soc Am A 2018; 35(2): 236-243. DOI: 10.1364/JOSAA.35.000236.
[6] Mao X, Li J, Wang F, Gao R, Li X, Xie Y. Fast design method of smooth freeform lens with an arbitrary aperture for collimated beam shaping. OSA Technical Digest 2019; JT5A.2. Washington: DC United States; 2019. ISBN: 9781-943580-60-6.
[7] Wei S, Zhu Z, Fan Z, Yan Y, Ma D. Double freeform surfaces design for beam shaping with non-planar wavefront using an integrable ray mapping method. Opt Express 2019; 27(19): 26757-26771. DOI: 10.1364ЮЕ.27.026757.
[8]
[9]
Kochengin SA, Oliker VI. Computational algorithms for constructing reflectors. Comput Vis Sci 2003; 6: 15-21. DOI: 10.1007/s00791-003-0103-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.
[10] 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.
[11] 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.
[12] 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.
[13] Wang X-J. On design of a reflector antenna II. Calc Variations Partial Differ Equ 2004; 20: 329-341. DOI: 10.1007/s00526-003-0239-4.
[14] Doskolovich LL, Bykov DA, Mingazov AA, Bezus EA. Optimal mass transportation and linear assignment problems in the design of freeform refractive optical elements generating far-field irradiance distributions. Opt Express 2019; 27(9): 13083-13097. DOI: 10.1364/OE.27.013083.
[15] Yadav NK. Monge-Ampere problems with non-quadratic cost function: application to freeform optics. Eindhoven: Technische Universiteit Eindhoven; 2018.
[16] Bykov DA, Doskolovich LL, Mingazov AA, Bezus EA. Optimal mass transportation problem in the design of freeform optical elements generating far-field irradiance distributions for plane incident beam. Appl Opt 2019; 58(33): 9131-9140. DOI: 10.1364/AO.58.009131.
[17] 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.
[18] 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.
[19] Oliker V, Rubinstein J, Wolansky G. Supporting quadric method in optical design of freeform lenses for illumination control of a collimated light. Adv Appl Math 2015; 62: 160-183. DOI: 10.1016/j.aam.2014.09.009.
[20] Glimm T. A rigorous analysis using optimal transport theory for a two-reflector design problem with a point source. Inverse Probl 2010; 26(4): 045001. DOI: 10.1088/02665611/26/4/045001.
[21] Mingazov AA, Doskolovich LL, Bykov DA, Kazanskiy NL. The two reflector design problem for forming a flat wavefront from a point source as an optimal mass transfer problem. Computer Optics 2019; 43(6): 968-975. DOI: 10.18287/2412-6179-2019-43-6-968-975.
[22] 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.
[23] Bykov DA, Doskolovich LL, Mingazov AA, Bezus EA, Kazanskiy NL. Linear assignment problem in the design of freeform refractive optical elements generating prescribed irradiance distributions. Opt Express 2018; 26(21): 2781227825. DOI: 10.1364/OE.26.027812.
[24] Makarov B, Podkorytov A. Real analysis: Measures, integrals and applications. London: Springer-Verlag; 2013.
[25] 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.
[26] TracePro - software for design and analysis of illumination and optical systems. Source:
(https://www. lambdares.com/tracepro/).
Сведения об авторах
Мингазов Альберт Айдарович, в 2010 году окончил Самарский государственный университет (СамГУ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Математика». В 2012 году окончил магистратуру НИУ ВШЭ по направлению «Математика», в 2015 окончил аспирантуру Санкт-Петербургского отделения математического института имени В. А. Стеклова. Кандидат физико-математических наук (2015), научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН - филиала ФНИЦ «Кристаллография и фотоника» РАН. E-mail: mingazov88@smail.com .
Досколович Леонид Леонидович, в 1989 году с отличием окончил Куйбышевский авиационный институт (КуАИ, ныне - Самарский национальный исследовательский университет имени академика С.П. Королёва) по специальности «Прикладная математика». Доктор физико-математических наук (2001 год), профессор, главный научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН (ИС-ОИ РАН) — филиала ФНИЦ «Кристаллография и фотоника» РАН, профессор кафедры технической кибернетики Самарского университета и ведущий научный сотрудник научно-исследовательской лаборатории прорывных технологий дистанционного зондирования Земли Самарского университета. Специалист в области дифракционной оптики, лазерных информационных технологий, нанофотоники. E-mail: leonid@ipsiras.ru .
Быков Дмитрий Александрович, в 2009 году с отличием окончил Самарский государственный аэрокосмический университет имени академика С.П. Королёва (СГАУ, ныне — Самарский национальный исследовательский университет имени академика С. П. Королёва) по специальности «Прикладная математика и информатика». Доктор физико-математических наук (2017 г.), старший научный сотрудник лаборатории дифракционной оптики Института систем обработки изображений РАН (ИСОИ РАН) - филиала ФНИЦ «Кристаллография и фотоника» РАН и доцент Самарского университета. Области научных интересов: оптика резонансных дифракционных структур, магнитооптика, электромагнитная теория дифракции. E-mail: bykovd@smail.com .
Бызов Егор Владимирович, 1988 года рождения. В 2014 году с отличием окончил обучение в магистратуре Самарского государственного аэрокосмического университета имени академика С.П. Королёва (СГАУ, ныне Самарский национальный исследовательский университет имени академика С. П. Королева) по направлению «Прикладные математика и физика». Область научных интересов: методы расчетов формирующей неизобра-жающей оптики для светодиодов. E-mail: egor.byzov@smail.com .
ГРНТИ: 29.31.29
Поступила в редакцию 29 сентября 2021 г. Окончательный вариант - 26 ноября 2021 г.
Support quadric method in non-imaging optics problems that can be reformulated as a mass transfer problem
A.A. Mingazov1, L.L. Doskolovich1-2, D.A. Bykov1-2, E.V. Byzov1 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
The article deals with problems of generating desired illumination patterns, formulated in a special way. More precisely, we consider problems that can be reformulated as a Monge-Kantorovich mass transfer problem with some cost function. For all problems of this type, we uniformly formulate the support quadric method and show that it coincides with the gradient method for finding the maximum of a certain concave function.
Keywords: non-imaging optics, geometric optics, inverse problem, Monge-Kantorovich problem, support quadric method.
Citation: Mingazov AA, Doskolovich LL, Bykov DA, Byzov EV. Support quadric method in non-imaging optics problems that can be reformulated as a mass transfer problem. Computer Optics 2022; 46(3): 353-365. DOI: I0.18287/2412-6179-C0-I055.
Acknowledgements: The work was funded by the Ministry of Science and Higher Education of the Russian Federation under a government project of the FSRC "Crystallography and Photonics" RAS (numerical implementation of the calculation algorithm) and the Russian Science Foundation under grant #18-19-00326 (proof of the coincidence with the gradient method for the corresponding functional).
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 department 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@ipsiras.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 .
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 .
Received September 29, 2021. The final version - November 26, 2021.
