CC BY
5
Алгоритм неупорядоченного распространения терагерцового волнового фронта для восстановления фазы при плотной многоплоскостной
регистрации данных
Е.Г. Циплакова1, А. Чопард 2'3, Н.С. Балбекин1, О.А. Смолянская1, Ж.-Б. Перро 4, Ж.-П. Гийе2, П. Моне2, Н.В. Петров1 1 Национальный исследовательский университет ИТМО, 197101, Россия, г. Санкт-Петербург, Кронверкский пр., д.49, лит. А;
2 Университет Бордо, 33405, Франция, г. Таланс, пл. Либерасьон, A31;
3 Компания Lytid, 92120, Франция, г. Монруж, ул. ла Фонтэн, д. 8;
4 Компания Optikan, 33000, Франция, г. Бордо, ул. Пермантад, д. 64
Аннотация
Итерационные алгоритмы восстановления фазы из дифракционных распределений в те-рагерцовом (ТГц) диапазоне частот являются многообещающим подходом вычислительной визуализации, способным обеспечить высокое пространственное разрешение восстанавливаемых изображений. Одним из часто применяемых алгоритмов стал SBMIR (англ. Single Beam Multiple Intensity Reconstruction), где в качестве исходных данных используются множественные распределения интенсивности дифракционного объектного поля. По сравнению с одноплоскостными методами, многоплоскостной подход позволяет ускорить сходимость, однако требует времязатратной регистрации данных на различном удалении приемника от объекта. Ранее нами был предложен метод записи в режиме однократного сканирования, который позволяет быстро получить исчерпывающий набор дифракционных распределений ТГц-объектного поля. В настоящей работе представлена апробация одного из современных алгоритмов на основе SBMIR и техники стохастического распространения волнового фронта (R-SBMIR, от англ. randomly ordered propagations) на экспериментальных данных ТГц-диапазона, снятых при помощи метода однократного сканирования. В отличие от ряда классических алгоритмов, которые могут потребовать серии численных экспериментов для определения оптимальных распределений интенсивности из большого набора данных, неупорядоченное распространение оценочного волнового фронта гарантирует восстановление изображений с наилучшим контрастом и разрешением без предварительной настройки параметров. Показано, что комплексное внедрение техники однократного сканирования с последующей обработкой данных при помощи R-SBMIR имеет потенциал для автоматизации метода восстановления фазы в ТГц-диапазоне.
Ключевые слова: восстановление фазы, терагерцовое излучение, итерационный алгоритм, метод стохастической оптимизации, фазовая визуализация, квантово-каскадный лазер.
Цитирование: Циплакова, Е.Г. Алгоритм неупорядоченного распространения терагер-цового волнового фронта для восстановления фазы при плотной многоплоскостной регистрации данных / Е.Г. Циплакова, А. Чопард, Н.С. Балбекин, О. А. Смолянская, Ж.-Б. Перро, Ж.П. Гийе, П. Моне, Н.В. Петров // Компьютерная оптика. - 2023. - Т. 47, № 6. - С. 901912. - DOI: 10.18287/2412-6179-CO-1253.
Citation: Tsiplakova EG, Chopard A, Balbekin NS, Smolyanskaya OA, Perraud GB, Guillet JP, Mounaix P, Petrov NV. An algorithm of unordered wavefront propagation in terahertz phase retrieval with dense multiplane data acquisition. Computer Optics 2023; 47(6): 901-912. DOI: 10.18287/2412-6179-CO-1253.
Введение
Одним из важных приложений излучения терагерцового (ТГц) диапазона в силу его уникальных свойств стали задачи визуализации [1]. ТГц-излучение способно проникать сквозь многие диэлектрики (полимеры, ткань, дерево и др.), которые непрозрачны в оптическом диапазоне спектра. Вместе с тем субмиллиметровое излучение не оказывает ионизирующего воздействия, что делает его перспективным для применения в задачах неразрушающего кон-
троля качества продукции, визуализации внутренних частей объектов [2, 3] и исследования биологических образцов [4, 5]. Также колебательно-вращательные переходы многих органических и неорганических молекул относятся к ТГц-диапазону, поэтому, наряду с геометрией объектов, появляется возможность идентифицировать их характерный химический состав [6]. При этом длины волн ТГц-излучения позволяют получать изображения с достаточным для решения ряда прикладных задач пространственным разрешением субмиллиметрового масштаба. Сово-
купность этих характеристик привела к росту числа исследований, связанных с применением ТГц-излучения в области бесконтактной визуализации объектов.
В настоящее время много работ в этой области направлено на развитие методов восстановления фазового распределения поля, необходимого для полного анализа трехмерных свойств объектов. Активно развивается техника ТГц-количественной фазовой визуализации с источниками как импульсного [7 - 9], так и непрерывного излучения [10 - 27]. Методы количественной фазовой визуализации с использованием источников непрерывного действия предполагают более простую и компактную конфигурацию приемопередающей системы. Среди них можно выделить цифровую голографию и интерферометрию [10 - 16], восстановление волнового фронта методом Гартмана [17], шерографию [18], птихографию [19, 20] и другие методы восстановления фазы [21 - 27].
Практическая реализация методов фазовой визуализации в ТГц-диапазоне была впервые продемонстрирована в работе [11], исследующей схему внеосе-вой голографии с пропускающим объектом. В качестве источника использовался диод Ганна (0,1 ТГц, мощность 35 мВт), а в качестве детектора - диод с барьером Шоттки в режиме растрового сканирования. Позже были продемонстрированы работы, где в качестве детектора использовались матричные микроболометрические приемники, позволяющие проводить ТГц-фазовую визуализацию в реальном времени [12, 13]. В настоящее время, благодаря усовершенствованию ТГц-оборудования, увеличению мощности источников и чувствительности детекторов, число исследований в этой области продолжает расти [14 -16]. Одним из примеров использования на практике высокомощного ТГц-источника для получения амплитудно-фазовых изображений является работа [16], обозревающая применение лазера на свободных электронах мощностью от 10 до 200 В в задачах цифровой и компьютерной голографии.
Другим подходом для решения фазовой задачи являются итерационные алгоритмы восстановления фазы [19 - 54], которые берут свое начало в работе Герчберга-Сакстона [28] и находят широкое применение: для вычисления структуры дифракционных оптических элементов [29, 30], в электронной микроскопии [31, 32], астрономии [33], для коррекции аберраций [34] и задач биологии [35, 36]. К преимуществам итерационных методов относятся использование однолучевой схемы регистрации, что делает их более энергоэффективными и устойчивыми к вибрациям по сравнению с методами внеосевой голографии; простота юстировки освещающего объект пучка; отсутствие мнимого изображения по сравнению с методами осевой голографии [27, 38]. Первые алгоритмы поиска фазы использовали в качестве исходных данных единственное заданное монохроматиче-
ское распределение интенсивности в некоторой плоскости, перпендикулярной оптической оси. Решение обратной задачи дифракции в этом случае может требовать большого числа итераций, а также иметь проблему сходимости к локальному минимуму. В связи с этим впоследствии были предложены многокадровые алгоритмы восстановления фазы, решающие проблему стагнации классических алгоритмов. В качестве входных данных в многокадровых методах используется набор распределений, зарегистрированный при вариации какого-либо параметра [39]: расстояния до объекта (многоплоскостные методы) [40 - 42], длины волны освещающего объект пучка [42 - 45], фазовой модуляции [46 - 48], степени оптической нелинейности [49 - 51].
В данной работе мы рассмотрим применение многоплоскостных алгоритмов, предполагающих использование нескольких распределений интенсивности, зарегистрированных в зоне дифракции Френеля на различном удалении от объекта. Многоплоскостной подход позволяет предотвратить проблему стагнации, характерную для алгоритмов, использующих единственное распределение интенсивности в качестве данных для восстановления поля в объектной плоскости. Одним из алгоритмов, использующих многоплоскостной подход, является SBMIR (от англ. Single Beam Multiple Intensity Reconstruction) [40]. Сбор данных для задачи восстановления фазы при помощи алгоритма SBMIR, как правило, предполагает последовательное перемещение детектора с некоторым продольным шагом в пространстве дифракционного поля объекта, а запись каждого отдельного распределения осуществляется в неподвижном положении детектора в течение некоторого времени для последующего усреднения данных с целью повышения отношения сигнал / шум. Будем сокращенно называть данный режим регистрации данных SM (от англ. Stop Motion). Выбор оптимальных положений детектора для записи дифракционных распределений в режиме SM представляет определенную сложность и часто должен производиться эмпирически, путем проведения серий экспериментальных измерений.
В работе [52] нами был предложен другой метод регистрации данных - экспресс-техника однократного сканирования (OTG от англ. On-The-Go). Использование эффективной конфигурации источника в виде квантово-каскадного лазера (ККЛ) и высокочувствительного матричного приемника позволило быстро получить исчерпывающий объем данных в ТГц-диапазоне частот с малым продольным шагом, порядка сотен мкм. Техника OTG дала возможность проводить подбор оптимального числа и расположения плоскостей регистрации данных на этапе численной обработки, исключая тем самым необходимость в многократном проведении физических измерений.
Несмотря на то, что техника OTG существенно упрощает процесс измерений, все еще остается опре-
деленная сложность подбора подходящих распределений интенсивности из их большого числа в плотном наборе данных для реализации алгоритма SBMIR. Для нахождения наилучших параметров восстановления волнового фронта часто требуется проведение серии численных расчетов.
Потенциальным решением данной проблемы может стать алгоритм на основе SBMIR, использующий распространение волнового фронта в стохастическом порядке [53, 54]. В нем дистанция распространения между плоскостями регистрации является варьируемой величиной. Преимуществами алгоритма стохастического распространения волнового фронта является то, что они обладают повышенной стойкостью к стагнации, а также отсутствие необходимости в тщательном подборе распределений интенсивности и расстояний между плоскостями регистрации. Наибольшую эффективность алгоритмы стохастического распространения показывают при работе с большим числом распределений интенсивности, поскольку в таком случае общее количество возможных путей распространения поля сильно возрастает, что приводит к ускоренной сходимости. Вместе с этим предложенная нами ранее техника OTG способна обеспечить подходящие наборы данных, которые содержат большое число плотно записанных распределений.
В рамках данной работы представлены результаты восстановления фазы из плотно зарегистрированных в объеме дифракционного поля ТГц-распределений интенсивности при помощи оптимизированного SBMIR-алгоритма, использующего стохастическое распространение волнового фронта R-SBMIR (от англ. randomly ordered propagations) [53]. Показано, что совмещение предложенной ранее экспресс-техники OTG с последующей численной обработкой данных при помощи алгоритма неупорядоченного распространения волнового фронта позволяет нивелировать необходимость подбора параметров регистрации и открывает возможность автоматизации метода итеративного восстановления фазы.
Используемые методы и подходы
К настоящему моменту разработано большое число систем генерации и детектирования ТГц-излучения [55]. Одним из достижений в области ТГц-технологий можно считать появление квантово-каскадных лазеров (ККЛ), способных обеспечить высокую выходную мощность, более 100 мВт, в режиме непрерывной генерации [56, 57]. При этом ККЛ способны генерировать излучение в широком диапазоне частот (1 - 5 ТГц). Использование ККЛ в задачах фазовой визуализации уже было исследовано в ряде работ [58 - 61]. В том числе мы исследовали применение ККЛ в качестве источника непрерывного излучения в задаче восстановления фазы из интенсивности [52]. Применение излучения ККЛ совместно с высокочувствительным матричным приемником позволило зарегистрировать распределения интенсивности с
высоким соотношением сигнал / шум, обеспечивая высокое качество восстановленных амплитудных и фазовых распределений.
В данной работе в качестве входных данных мы использовали дифракционные распределения, снятые при помощи аналогичной конфигурации источника и приемника, как в работе [52]. В качестве источника использовался ККЛ (частота излучения V = 2,5 ТГц, длина волны X = 119,91 мкм), а в качестве детектора -матричный приемник от 12Б на основе массива микроболометров, сопряженного с антеннами (288 х 384, размер пикселя Ах = 50 мкм). Поскольку разрешение достигает порядка длины волны, высокочастотное излучение ККЛ позволяет уверенно исследовать объекты с меньшим размером неоднородности, порядка 150 - 200 мкм, по сравнению с излучением используемого нами ранее диода Ганна (частота = 0,287 ТГц, длина волны = 1044,57 мкм) [62, 63]. Другим преимуществом используемой комбинации источника и приемника является возможность быстрой регистрации распределений интенсивности в режиме ОТв. В работе [52] мы приводим сравнение техник многоплоскостной регистрации данных, 8М и ОТв, и результаты последующей обработки при помощи известного алгоритма восстановления фазы 8ВМГЯ.
На рис. 1 представлены графики зависимости числа кадров Ыям и Ыото, снятых при помощи соответствующих техник сбора данных, от времени t и координаты приемника г, отсчитываемой от первой плоскости регистрации. Графики обоих методов 8М и ОТв построены в промежутке перемещения приемника от 0 до 10 мм. Для режима регистрации данных 8М введем также количество кадров, используемых для усреднения в каждой плоскости, N (отмечены на рис. 16 для рассматриваемых плоскостей регистрации). В связи с неавтоматизированным процессом измерений значение N варьировалось в пределах от 15 до 50 в каждой плоскости. Таким образом, общее время регистрации в 8М режиме зависит от Щ, а итоговое распределение интенсивности в плоскости рассчитывается как:
JNSM (У)="
Z N^ ( y)
Nj
При последовательном перемещении детектора с продольным шагом Аг = 2 мм и записи в стационарном состоянии, в выбранном промежутке расстояний г от 0 до 10 мм было получено Щш=6 распределений интенсивности. С учетом съемки кадров для усреднения для этого потребовалось около 22 секунд (рис. 1а).
Предложенная в работе [52] экспресс-техника ОТв реализует сбор данных при непрерывном перемещении детектора с постоянной скоростью (10 мм/с), а запись осуществляется с частотой 25 кадр/с. Таким образом, в промежутке 10 мм перемещения детектора за время t = 1 с было зарегистрировано Щота=25 дифракционных распределений (рис. 16).
a) z, мм 10 20 с 6)
Рис. 1. Сравнение техник сбора данных. Зависимость числа распределений интенсивности, полученных в режиме (a) SM (stop motion) с регистрацией Njраспределений в каждой плоскости Nsm для последующего усреднения; (б) OTG (on-the-go),
от времени t и координаты приемника z
При этом, как уже было отмечено в работе [52], ввиду использования высокочувствительной камеры последнего поколения на основе массива микроболометров, сопряженного с антеннами, распределения интенсивности, снятые в режиме 8М в одну экспозицию, были пригодны для использования в качестве исходных данных в 8ВМЖ. Тем не менее, как было показано в [52], скорость сбора данных в режиме ОТв все равно превосходит скорость сбора данных в режиме 8М, даже без регистрации кадров для усреднения.
В настоящей работе исследовалась схема регистрации данных на пропускание [26]. В качестве объекта использовалась пластиковая карта с нанесенными на нее рельефными изображениями (рис. 2). В данной работе мы используем данные, позволяющие сформировать два фрагмента рельефных изображений данного объекта, а именно: изображения идентификационного кода пластика (далее - объект 1) и стрелки (объект 2) (выделены красным на рис. 2).
Рис. 2. Фотография исследуемого объекта
Сбор данных для восстановления амплитудно-фазовых изображений этих объектов проводился в режиме ОТв на удалении х1 = 11,9 мм и 11,4 мм от объектов соответственно. При скорости перемещения детектора 10 мм/с были зарегистрированы наборы данных, содержащие порядка Ыото = 500 (далее индекс ОТв будет опускаться, предполагая использование наборов данных, снятых методом ОТв по умолчанию) распределений интенсивности дифракционного поля с про-
дольным шагом /\2=400 мкм. В работе [52] уже были получены фрагменты изображений данного тестового объекта, но применение алгоритма Я-8ВМГЯ для обработки данных ранее не рассматривалось.
Основой итерационных алгоритмов восстановления фазы, использующих многоплоскостной подход, является процедура распространения оценочной функции волнового фронта между плоскостями регистрации и плоскостью объекта. Моделирование распространения волнового фронта между плоскостями регистрации данных в зоне дифракции Френеля производится при помощи метода углового спектра [64]. В дискретном представлении уравнение для расчета поля через представление в виде углового спектра плоских волн описывается следующим образом:
и (х, у, х ) = Т-1 [ Т [и 0 (у')] • Н ((, Л, А)],
где и (X, у') - распределение комплексной амплитуды поля в исходной плоскости; (X, у') - координаты в исходной плоскости; (х, у) - координаты в конечной плоскости; Ах - расстояние распространения; Т и Т-1 -операторы быстрого прямого и обратного Фурье-преобразования.
Передаточная функция Н (/х, /у, Ах), описывающая процессы распространения, представляется в следующем виде:
Н (/х, /у, X ) =
ехр^^Ах-(Х/х)) +(Х/у)) ^, /2 + Х ;
где fx,fy - координаты в частотной плоскости; 1 -длина волны излучения.
Для удобства введем оператор углового спектра Fas, описывающий распространение волнового фронта между плоскостями:
U(x,y,z = z0 +Az) = Fas (Uo (xyz0),Az).
Ранее в работе [52] восстановление волнового фронта проводилось при помощи алгоритма SBMIR. В данной работе мы сравниваем эффективность восстановления волнового фронта из плотно зарегистрированных наборов распределений с помощью алгоритмов SBMIR и R-SBMIR. Ниже мы рассмотрим особенности данных алгоритмов более подробно.
SBMIR реализует следующую последовательность шагов: 1. На первом шаге итерации формируется начальное приближение поля в первой плоскости регистрации:
Vi (х,y,zi) = ^Ii (х,y,zi) exp(ф (x,y,zi)) =
= |«i (x, y, zi )| exp (ф (x, y, zi)),
где щ (x,y, zi) является истинным комплексным распределением объектного поля в первой плоскости, I (x, y, zi) - зарегистрированное в первой плоскости распределение интенсивности, ф1 (x, y, zi) - начальное фазовое приближение, представляющее собой константное значение или набор случайных значений.
2. Комплекснозначное распределение Vi (x,y, zi) численно распространяется во вторую плоскость регистрации данных на расстояние Az, формируя распределение в данной плоскости:
V2 (x, y, z2 = zi + Az ) = Fas (Vi (x, y, zi), Az).
3. Далее вычисляется новая оценочная функция поля путем замены амплитудной части сформированного комплекснозначного распределения на квадратный корень из интенсивности, снятой во второй плоскости регистрации:
V2
= 4h ^ VITexp (гф2).
4. Шаги (2) и (3) повторяются последовательно для каждой плоскости регистрации вплоть до последней с индексом N, где в результате получаем комплексное распределение Vn = 4INexp ^n ).
5. Далее распространение оценочного поля и последующая замена модуля комплексной амплитуды в каждой плоскости регистрации на соответствующий корень из интенсивности производится в обратном направлении:
Vn-i (x, y, zn-i ) = Fas ( Vn (x, y, Zn ), - Az);
Vn-i = VI
; = J I
ex
p (n-i ).
Vi (x,y,zi) = Fas (u2 (x2,y2,z2), -Az);
v! =JTi Д = Vliexp(ф).
Vi
6. Можем оценить комплексное распределение поля
в плоскости объекта:
vj (x,y,0) = Fas (vi (x,y,zi), - zi ) = |v0| exp(гф0).
Здесь верхний индекс обозначает номер итерации. Таким образом, данная последовательность шагов представляет собой одну итерацию алгоритма SBMIR.
Особенностью алгоритма R-SBMIR является распространение оценочной функции поля между плоскостями регистрации в стохастическом порядке, так что расстояние распространения не принимает фиксированного значения и не равняется продольному шагу регистрации камерой Az, как в SBMIR (рис. 3).
Объектная
Рис. 3. Последовательность распространения оценочной функции поля между плоскостями регистрации в ходе итерационного процесса (а) ББМШ и (6) К-ББМШ
Рассмотрим подробнее итерационный процесс алгоритма К-8ВМ1Я.
1. Перед началом итеративного процесса формируются 2 вектора N и Щь, которые задают неупорядоченную последовательность использования плоскостей регистрации по ходу прямого и обратного распространения оценочного волнового фронта соответственно. Из второго вектора N исключается элемент, равный последнему в Щ После этого векторы объединяются в один: ЩШп<1 = N + Щь, который представляет собой последовательность использования плоскостей регистрации в ходе одной итерации Я-8ВМ1Я.
2. Вектору Щгапа ставится в соответствие вектор координат плоскостей регистрации в сформированном порядке:
гтпа = г + (уЩгаМ ~ !■) ' Аг .
3. Далее формируется вектор путей неупорядоченного распространения оценочного волнового фронта между плоскостями регистрации:
2
N-i
Az _ ¡ \z1 2
rand l rand rand
где Az'mnd _ z'^m - z'mnd. Верхний индекс обозначает порядковый номер элемента в векторе.
4. Используя оператор распространения Fas и вектор Azmnd, производится распространение оценочного волнового фронта между соседними парами плоскостей INrand (x, y). При этом на каждом шаге итерации так же, как в SBMIR, производится замена амплитудной части оценочной функции на соответствующее экспериментально измеренное распределение.
5. На последнем шаге итерации комплексное распределение поля фокусируется в плоскость объекта посредством распространения на расстояние -zrNand .
6. Для последующих итераций выполняются описанные шаги 1- 5, но дополнительно после формирования вектора zrand оценочная функция поля должна быть распространена на расстояние z1rand . Оценка сходимости алгоритмов производилась с
помощью вычисления нормированного среднеквад-ратического отклонения (НСКО) [65]:
E _
У
Z(l v(x,У)|- u(x,y) |)
Z |u (x, y )|
где V (х, у) и и (х, у) - комплексная амплитуда оценочного и истинного распределения поля в объектной плоскости соответственно.
Результаты и обсуждение
Ранее мы показали, что использование распределений интенсивности с неоптимальным шагом регистрации Дг может приводить к стагнации алгоритма 8ВМ1Я, вносить дополнительные неоднородности и уменьшать контраст восстановленных амплитудных и фазовых изображений [21, 52]. В свою очередь, данная проблема явно фигурировала при работе с данными плотного набора распределений, снятых методом ОТв [52]. Таким образом, хотя предложенная техника экспресс-регистрации позволяет нивелировать необходимость точной настройки шага детектирования на этапе постановки эксперимента, выбор оптимальных распределений из большого набора данных для внедрения в алгоритм 8ВМ1Я все еще представляется трудоемкой задачей, требующей отдельного внимания. В настоящей работе мы исследуем применение техники неупорядоченного распространения волнового фронта для решения проблемы времязатратного подбора подходящих распределений из большого объема данных, снятых методом ОТв, на этапе численной обработки при помощи 8ВМГЯ.
Распределения дифракционного поля удобно ха-растеризовать при помощи чисел Френеля:
F * _
D2 4Xz
где Б - характерный размер объекта, 1 - длина волны излучения, х - расстояние от объекта до плоскости регистрации дифракционного распределения. Для излучения ККЛ и исследуемых объектов характеризующие числа Френеля сформированного поля вычислялись с использованием наибольшего линейного размера объекта: 10 мм - для объекта 1 и 5 мм - для объекта 2. На рис. 4а, б приведены продольные осевые срезы картины дифракции для рассматриваемых объектов. Из-за различий в размере объектов продольные распределения дифракционного поля относятся к разным зонам дифракции Френеля. Это, в свою очередь, обуславливает необходимость подстройки продольного шага Аг между используемыми распределениями в каждом отдельном случае для корректного восстановления волнового фронта при помощи классического алгоритма 8ВМ1Я.
Мы провели серию восстановлений волнового фронта при помощи алгоритмов 8ВМ1Я и Я-8ВМ1Я из наборов данных, снятых методом ОТв. В ходе численной обработки экспериментальных данных при помощи итерационных алгоритмов было установлено, что использование распределений интенсивности, зарегистрированных на больших расстояниях от объекта, негативно сказывается на пространственном разрешении, поскольку они обладают низкой числовой апертурой. Поэтому восстановление волнового фронта производилось с использованием только близко зарегистрированных к объектам распределений интенсивности. Предельное расстояние регистрации распределений интенсивности определяется длиной волны используемого источника, геометрией объекта и размером сенсора приемника. Кроме того, в ходе численной обработки было установлено, что использование в Я-8ВМ1Я распределений интенсивности, снятых в пределах протяженного интервала расстояний от объекта, способствует ухудшению качества восстанавливаемых амплитудных и фазовых изображений. Для восстановления распределений с минимальным уровнем шума при помощи Я-8ВМ1Я необходимо было установить ограничение либо на расстояние распространения в неупорядоченной последовательности плоскостей регистрации, либо на число используемых распределений интенсивности. Для исследуемых объектов 1 и 2 мы использовали распределения интенсивности, снятые на расстоянии не более 30 мм. Число распределений интенсивности варьировалось в пределах от 5 до 20 в зависимости от используемого продольного шага Дг между плоскостями регистрации.
Вначале в качестве исходных данных мы использовали 20 распределений с минимальным шагом Д7 = 400 мкм. В этом случае восстановленные при помощи 8ВМЖ фазовые и амплитудные распределения обоих объектов содержат большое число неоднородностей (рис. 4е, ж ,д, и). Из-за недостаточного интервала распространения заменяемая амплитудная часть в каждой
новой плоскости в незначительной степени отличалась от предыдущего распределения, не позволяя сформировать изменения в новом фазовом приближении, что привело к стагнации алгоритма 8ВМЖ. На графике НСКО (рис. 4л) объекта 1 синей кривой, демонстрирующей возрастающую тенденцию, показан случай восстановления с рассматриваемыми параметрами Аг = 400 мкм, N = 20. В итоге при обработке данных с
помощью 8ВМЖ потребовалось провести несколько численных экспериментов с разным количеством задействованных плоскостей и значением Дг, прежде чем был достигнут наилучший результат. При этом не все рассмотренные комбинации шага Дг и числа плоскостей N смогли обеспечить сходимость классического алгоритма 8ВМ1Я. Достаточным значением НСКО было условно принято значение 0,18.
а)
Рис. 4. (а, 6) Продольное сечение амплитуды зарегистрированного дифракционного поля по центру (у = 0мм) при использовании
в качестве объекта изображения идентификационного кода пластика и стрелки соответственно; (в - е) амплитудные и (ж - к) фазовые распределения, восстановленные с помощью алгоритма (в, д, ж, и) ББМШ при Аг = 400 мкм
и (г, е, з, к) К-ББМШ; (л) НСКО в зависимости от числа итераций при использовании в качестве объекта изображения идентификационного кода; (м) фазовые сечения части изо6ражения, выделенные вертикальной 6елой линией на изо6ражении (к)
ЧГ1
I
!
Л^ Номер итерации м) у, мм
Рис. 4 (продолжение). (а, б) Продольное сечение амплитуды зарегистрированного дифракционного поля по центру (у = 0мм) при использовании в качестве объекта изображения идентификационного кода пластика и стрелки соответственно; (в - е) амплитудные и (ж - к) фазовые распределения, восстановленные с помощью алгоритма (в, д, ж, и) 8БЫШ при Аг = 400 мкм и (г, е, з, к) К-8БЫШ; (л) НСКО в зависимости от числа итераций при использовании в качестве объекта изображения идентификационного кода; (м) фазовые сечения части изображения, выделенные вертикальной белой линией на изображении (к)
Вследствие того, что в плоскости объекта мы не имели точно измеренной интенсивности, сравнение проводилось с распределением, усредненным по нескольким наилучшим восстановленным амплитудным изображениям. Параметры, которые обеспечили сходимость 8ВМ1Я, определены методом перебора: Аг = 2,4 мм, число плоскостей N = 5 (рис. 4л, желтая кривая).
Далее было исследовано восстановление волнового фронта объектов при помощи алгоритма Я-8ВМ1Я, используя так же, как и в первом расчете, 20 ближайших к объекту плоскостей с минимальным шагом Аг = 400 мкм. Оптимальное значение НСКО было сразу достигнуто, при этом за меньшее число итераций. Восстановленные распределения поля с использованием алгоритма Я-8ВМ1Я из неотсортированного набора первых 20 картин интенсивности, снятых методом ОТв, представлены на рис. 4г, з, е, к. В данном случае распределения четко определяют границы и форму объекта, если сравнивать с изображениями, полученными методом 8ВМ1Я из аналогичного набора данных.
Рассматривая восстановление волнового фронта, прошедшего объект 2, используемые значения Аг и N не смогли обеспечить сходимость алгоритма 8ВМ1Я. Фазовое запаздывание, вносимое рельефом объекта 2, отмечено на рис. 4м для разных параметров 8ВМ1Я в сравнении с Я-8ВМ1Я. На рис. 4к линией отмечено положение сравниваемых фазовых сечений. Числа Френеля, характеризующие дифракционные распре-
деления поля объекта 2 на заданном расстоянии, имели малые значения, в пределах от 1,7 до 4,6, в связи с чем спекл-картины записанных распределений интенсивности в малой степени отличались друг от друга (рис. 4б) и восстановление волнового фронта при помощи 8ВМ1Я представлялось затруднительным. Алгоритм Я-8ВМ1Я позволил получить наилучший фазовый контраст восстановленного изображения объекта благодаря стохастическому характеру распространения поля между большим числом распределений (рис. 4к). Среди большого количества пар плоскостей, между которыми производится распространение поля в ходе итерационного процесса Я-8ВМ1Я, присутствуют комбинации с достаточным различием в амплитудных распределениях, что способствует улучшению нового фазового приближения. При этом с увеличением числа распределений интенсивности, снятых в пределах расстояний, которые бы обеспечили требуемое содержание высоких пространственных частот, эффективность алгоритма Я-8ВМ1Я повышается, поскольку общее число шагов и способов распространения оценочного волнового фронта возрастает. Поэтому метод ОТв хорошо применим в качестве техники сбора данных для последующей автоматизированной обработки с помощью Я-8ВМ1Я. Также для повышения разрешения восстанавливаемых амплитудных и фазовых распределений в случае использования дифракционных картин, характеризуемых малыми числами Френеля, и детектора с малой площадью сенсора, возможно применение
техники экстраполяции данных за пределы области регистрации [66 - 68].
Нужно отметить, что фазовое запаздывание, вносимое объектом 2 (рельефной стрелкой) в верхней ее части, так и не было до конца восстановлено, поскольку зондирующий пучок ККЛ содержал фазовую сингулярность в своем распределении в области расположения объекта [52]. Тем не менее, данное искажение в восстановленных изображениях не связано с результатом применения алгоритма R-SBMIR.
Вывод
В работе исследовано применение техники стохастического распространения волнового фронта для восстановления фазы ТГц-объектной волны при плотной многоплоскостной регистрации данных. Использование алгоритма неупорядоченного распространения оценочного волнового фронта R-SBMIR между плоскостями регистрации позволяет значительно упростить процесс обработки данных, исключая необходимость в проведении серии численных экспериментов для подбора оптимальных параметров сбора данных, а в случае использования техники OTG - восстановления изображений: продольного шага между распределениями интенсивности и их числа. Вместе с тем большой объем экспериментальных данных, полученный с помощью техники OTG, повышает эффективность работы стохастического алгоритма благодаря увеличенному числу путей распространения поля. Таким образом, техника OTG вместе с последующей численной обработкой данных при помощи высокоскоростного алгоритма неупорядоченного распространения могут выступать взаимодополняющими техниками для автоматизации метода восстановления фазы из дифракционных распределений.
Благодарности
Работа выполнена при поддержке международного гранта РФФИ № 21-52-15035 НЦНИ_а.
References
[1] Mittleman DM. Twenty years of terahertz imaging. Opt Express 2018; 26(8): 9417-9431. DOI: 10.1364/OE.26.009417.
[2] Amenabar I, Lopez F, Mendikute A. In introductory review to THz non-destructive testing of composite mater. J Infrared Millim Terahertz Waves 2013; 34(2): 152-169. DOI: 10.1007/s10762-012-9949-z.
[3] Zeitler JA, Gladden LF. In-vitro tomography and nondestructive imaging at depth of pharmaceutical solid dosage forms. Eur J Pharm Biopharm 2009; 71(1): 2-22. DOI: 10.1016/j.ejpb.2008.08.012.
[4] Son JH, ed. Terahertz biomedical science and technology. Boca Raton: CRC Press; 2014. ISBN: 978-1-4665-7045-0.
[5] Wan M, Healy JJ, Sheridan JT. Terahertz phase imaging and biomedical applications. Opt Laser Technol 2020; 122: 105859. DOI:10.1016/j.optlastec.2019.105859.
[6] Mohr T, et al. Two-dimensional tomographic terahertz imaging by homodyne self-mixing. Opt Express 2015; 23(21): 27221-27229. DOI: 10.1364/OE.23.027221.
[7] Petrov NV, et al. Application of terahertz pulse timedomain holography for phase imaging. IEEE Trans Terahertz Sci Technol 2016; 6(3): 464-472. DOI: 10.1109/TTHZ.2016.2530938.
[8] Karpowicz N, et al. Comparison between pulsed terahertz time-domain imaging and continuous wave terahertz imaging. Semicond Sci Technol 2005; 20(7): S293. DOI: 10.1088/0268-1242/20/7/021.
[9] Johnson JL, Dorney TD, Mittleman DM. Interferometric imaging with terahertz pulses. IEEE J Sel Top Quantum Electron 2001; 7(4): 592-599. DOI: 10.1109/2944.974230.
[10] Li Z, et al. Terahertz synthetic aperture in-line holography with intensity correction and sparsity autofocusing reconstruction. Photonics Res 2019; 7(12): 1391-1399. DOI: 10.1364/prj.7.001391.
[11] Mahon RJ, Murphy JA, Lanigan W. Digital holography at millimetre wavelengths. Opt Commun 2006; 260(2): 469473.
[12] Ding SH, et al. Continuous-wave terahertz digital holography by use of a pyroelectric array camera. Opt Lett 2011; 36(11): 1993-1995.
[13] Zolliker P, Hack E. THz holography in reflection using a high resolution microbolometer array. Opt Express 2015; 23(9): 10957-10967.
[14] Locatelli M, et al. Real-time terahertz digital holography with a quantum cascade laser. Sci Rep 2015; 5(1): 13566.
[15] Zhang Y, et al. Lensless Fourier-transform terahertz digital holography for real-time full-field phase imaging. Photonics Res 2022; 10(2): 323-331. DOI: 10.1364/PRJ.435769.
[16] Choporova Y, Knyazev B, Pavelyev V. Holography with high-power CW coherent terahertz source: optical components, imaging, and applications. Light: Advanced Manufacturing 2022; 3(3): 525-541. DOI: 10.37188/lam.2022.031.
[17] Richter H, et al. Terahertz wavefront measurement with a Hartmann sensor. Appl Phys Lett 2012; 101(3): 031103. DOI: 10.1063/1.4737164.
[18] Agour M, et al. Terahertz referenceless wavefront sensing by means of computational shear-interferometry. Opt Express 2022; 30(5): 7068-7081. DOI: 10.1364/oe.450708.
[19] Valzania L, et al. Terahertz ptychography. Opt Lett 2018; 43(3): 543-546.
[20] Rong L, et al. Continuous-wave terahertz reflective ptychography by oblique illumination. Opt Lett 2020; 45(16): 4412-4415. DOI: 10.1364/OL.400506.
[21] Petrov NV, et al. The features of optimization of a phase retrieval technique in THz frequency range. Speckle 2012: V Int Conf on Speckle Metrology 2012; 8413: 387-391. DOI: 10.1117/12.978688.
[22] Junkin G. Planar near-field phase retrieval using GPUs for accurate THz far-field prediction. IEEE Trans Antennas Propag 2012; 61(4): 1763-1776. DOI: 10.1109/TAP.2012.2220324.
[23] Gao X, Li C, Fang GY. The realization of terahertz image reconstruction with high resolution based on the amplitude of the echoed wave by using the phase retrieval algorithm. Chinese Phys Lett 2013; 30(6): 068401. DOI: 10.1088/0256-307X/30/6/068401.
[24] Ren Y, et al. Development of terahertz two-dimensional phase gratings for multiple beam generation based on a high-accuracy phase retrieval algorithm. Opt Express 2021; 29(12): 17951-17961. DOI: 10.1364/OE.425838.
[25] Rong L, et al. Transport of intensity equation-based terahertz lensless full-field phase imaging. Opt Lett 2021; 46(23): 5846-5849. DOI: 10.1364/ol.442625.
[26] Petrov NV, et al. Terahertz multiple-plane phase retrieval. Digital Holography and Three-Dimensional Imaging 2020: HF4G.8. DOI: 10.1364/DH.2020.HF4G.8.
[27] Jin X, et al. Iterative denoising phase retrieval method for twin-image elimination in continuous-wave terahertz in-line digital holography. Opt Lasers Eng 2022; 152: 106986.
[28] Gerchberg RW. A practical algorithm for the determination of phase from image and diffraction pictures. Optik 1972; 35(2): 237-246.
[29] Kotlyar VV, Khonina SN, Soifer VA. Iterative calculation of diffractive optical elements focusing into a three-dimensional domain and onto the surface of the body of rotation. J Mod Opt 1996; 43(7): 1509-1524. DOI: 10.1080/09500349608232822.
[30] Pavelyev VS, Soifer VA, Duparre M, Kowarschik R, Ludge B, Kley B. Iterative calculation, manufacture and investigation of DOE forming unimodal complex distribution. Opt Lasers Eng 1998; 29(4-5): 269-279. DOI: 10.1016/S0143-8166(97)00115-2.
[31] Miao J, Charalambous P, Kirz J, Sayre D. Extending the methodology of X-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens. Nature 1999; 400: 342-344.
[32] Zuo J, Vartanyants I, Gao M, et al. Atomic resolution imaging of a carbon nanotube from diffraction intensities. Science 2003; 300: 1419.
[33] Fienup C, Dainty J. Phase retrieval and image reconstruction for astronomy. Ch 7. In Book: Stark H, ed. Image Recovery: Theory and Application. Academic Press; 1987: 231-275.
[34] Fienup JR, et al. Hubble Space Telescope characterized by using phase-retrieval algorithms. Appl Opt 1993; 32(10): 1747-1767.
[35] Shevkunov I, et al. Super-resolution microscopy for biological specimens: lensless phase retrieval in noisy conditions. Biomed Opt Express 2018; 9(11): 5511-5523.
[36] Shechtman Y, et al. Phase retrieval with application to optical imaging: a contemporary overview. IEEE Signal Process Mag 2015; 32(3): 87-109.
[37] Boikov IV, Zelina YaV. Approximate methods of solving amplitude-phase problem for continuous signals [In Russian]. Measurement Techniques 2021; 5: 37-46. DOI: 10.32446/0368-1025it.2021-5-37-46.
[38] Rong L, Pan F, Xiao W, Li Y, Wang F. Twin image elimination from two in-line holograms via phase retrieval. Chinese Opt Lett 2012; 10(6): 060902.
[39] Nalegaev SS, Petrov NV, Bespalov VG. Special features of iteration methods for phase problem in optics [In Russian]. Nauchno-Tekhnicheskii Vestnik Informatsionnykh Tekhnologii, Mekhaniki i Optiki 2012; 6(82): 30-35.
[40] Pedrini G, Osten W, Zhang Y. Wave-front reconstruction from a sequence of interferograms recorded at different planes. Opt Lett 2005; 30(8): 833-835. DOI: 10.1364/ol.30.000833.
[41] Allen LJ, Oxley MP. Phase retrieval from series of images obtained by defocus variation. Opt Commun 2001; 199(1-4): 65-75. DOI: 10.1016/S0030-4018(01)01556-5.
[42] Petrov NV, Bespalov VG, Volkov MV. Wavefront reconstruction with a reference-free digital CCD-registration of multispectral speckle-patterns [In Russian]. Nanosyst Phys Chem Math 2011; 2(1): 82-90.
[43] Bao P, et al. Phase retrieval using multiple illumination wavelengths. Opt Lett 2008; 33(4): 309-311. DOI: 10.1364/OL.33.000309.
[44] Petrov NV, Bespalov VG, Gorodetsky AA. Phase retrieval method for multiple wavelength speckle patterns. Proc SPIE 2010; 7387: 73871T. DOI: 10.1117/12.871433.
[45] Chen N, Yeom J, Lee B. Optimized phase retrieval algorithm with multiple illuminations. In Book: Osten W, ed. Fringe 2013: 7th international workshop on advanced optical imaging and metrology. - Berlin, Heidelberg: Springer; 2014: 337-340. DOI: 10.1007/978-3-642-36359-7_60.
[46] Katkovnik V, et al. Computational wavelength resolution for in-line lensless holography: phase-coded diffraction patterns and wavefront group-sparsity. Proc SPIE 2017; 10335: 1033509. DOI: 10.1117/12.2269327.
[47] Shevkunov I, et al. Super-resolution microscopy for biological specimens: lensless phase retrieval in noisy conditions. Biomed Opt Express 2018; 9(11): 5511-5523. DOI: 10.1364/BOE.9.005511.
[48] Katkovnik V, et al. Computational super-resolution phase retrieval from multiple phase-coded diffraction patterns: simulation study and experiments. Optica 2017; 4(7): 786794. DOI: 10.1364/OPTICA.4.000786.
[49] Ivanov VY, Sivokon VP, Vorontsov MA. Phase retrieval from a set of intensity measurements: theory and experiment. J Opt Soc Am A 1992; 9(9): 1515-1524. DOI: 10.1364/JOSAA.9.001515.
[50] Lu CH, et al. Phase retrieval using nonlinear diversity. Appl Opt 2013; 52(10): D92-D96. DOI: 10.1364/AO.52.000D92.
[51] Lu JT, Lu CH, Fleischer JW. Enhanced phase retrieval using nonlinear dynamics. Opt Express 2016; 24(22): 2509125102.
[52] Chopard A, et al. Single-scan multiplane phase retrieval with a radiation of terahertz quantum cascade laser. Appl Phys B 2022; 128(3): 63. DOI: 10.1007/s00340-022-07787-x.
[53] Xing C, Qi F, Guo S. Improved terahertz phase imaging with single-beam multiple-plane reconstruction. 46th Int Conf on Infrared, Millimeter and Terahertz Waves (IRMMW-THz) 2021: 1-2.
[54] Xu C, et al. Enhancing multi-distance phase retrieval via unequal interval measurements. Photonics 2021; 8(2): 48. DOI: 10.3390/PHOTONICS8020048.
[55] Valusis G, et al. Roadmap of terahertz imaging 2021. Sensors 2021; 21(12): 4092. DOI: 10.3390/s21124092.
[56] Liang G, Liu T, Wang QJ. Recent developments of terahertz quantum cascade lasers. IEEE J Sel Top Quantum Electron 2016; 23(4): 1-18. DOI: 10.1109/JSTQE.2016.2625982.
[57] Wang X, et al. High-power terahertz quantum cascade lasers with~ 0.23 W in continuous wave mode. AIP Advances 2016; 6(7): 075210. DOI: 10.1063/1.4959195.
[58] Hack E, Zolliker P. Terahertz holography for imaging amplitude and phase objects. Opt Express 2014; 22(13): 16079-16086. DOI: 10.1364/OE.22.016079.
[59] Huang H, et al. Continuous-wave terahertz multi-plane inline digital holography. Opt Lasers Eng 2017; 94: 76-81. DOI: 10.1016/J.OPTLASENG.2017.03.005.
[60] Locatelli M, et al. Real-time terahertz digital holography with a quantum cascade laser. Sci Rep 2015; 5(1): 13566. DOI: 10.1038/srep13566.
[61] Deng Q, et al. High-resolution terahertz inline digital holography based on quantum cascade laser. Opt Eng 2017; 56(11): 113102. DOI: 10.1117/1.OE.56.11.113102.
[62] Petrov NV, et al. Terahertz diffractive reflection phase imaging. 2020 45th Int Conf on Infrared, Millimeter, and Terahertz Waves (IRMMW-THz) 2020: 1-1. DOI: 10.1109/IRMMW-THz46771.2020.9370854.
[63] Petrov NV, et al. Terahertz phase retrieval imaging in reflection. Opt Lett 2020; 45(15): 4168-4171. DOI: 10.1364/OL.397935.
[64] Mandel L, Wolf E. Optical coherence and quantum optics. Cambridge University Press; 1995.
[65] Fienup JR. Invariant error metrics for image reconstruction. Appl Opt 1997; 36(32): 8352-8357. DOI: I0.1364/A0.36.008352.
[66] Huang Z, Kuang C, Xu L, Cao L. Multiplane digital holography based on extrapolation iterations. Opt Commun 2021; 481: 126526.
[67] Latychevskaia T, Fink HW. Resolution enhancement in digital holography by self-extrapolation of holograms. Opt Express 2013; 21(6): 7726-7733.
[68] Balbekin NS, Kulya MS, Belashov AV, Gorodetsky A, Petrov NV. Increasing the resolution of the reconstructed image in terahertz pulse time-domain holography. Sci Rep 2019; 9(1): 180.
Сведения об авторах
Циплакова Елизавета Георгиевна, 2000 года рождения, окончила бакалавриат Национального исследовательского университета ИТМО по специальности «Фотоника и оптоинформатика» в 2022 году. В настоящее время обучается на первом курсе магистратуры Университета ИТМО по специальности «Техническая физика». Область научных интересов: алгоритмы восстановления фазы, терагерцовая оптика, цифровая голография. E-mail: e.tsiplakova@itmo.ru .
Чопард Адриен, 1993 года рождения, получил диплом инженера по оптике и фотонике в Высшей школе Института оптики (Institut d'Optique Graduate School (IOGS-Saclay)), Палезо, Франция, в 201б году и степень магистра в области нанотехнологий в Королевском технологическом институте Швеции (Royal Institute of Technology of Sweden), Стокгольм, Швеция. В 2021 году защитил диссертацию в области физики, направленную на разработку прикладных решений в области ТГц-технологий, при сотрудничестве между лабораторией IMS (UMR CNRS 5218) Университета Бордо, Таланс, Франция и компанией Lytid, Париж, Франция; в настоящее время отвечает за разработку прикладных решений в компании Lytid. E-mail: adrien.chopard@u-bordeaux.fr .
Балбекин Николай Сергеевич, 1989 года рождения, к.ф.-м.н., заведующий учебными лабораториями Университета ИТМО, заместитель руководителя лаборатории цифровой и изобразительной голографии. Автор 75 публикаций. Область научных интересов: аналоговая голография, цифровая голография, терагерцовая спектроскопия. E-mail: nbalbekin@niuitmo.ru .
Смолянская Ольга Алексеевна, к.ф.-м.н., доцент университета ИТМО, является заведующей Лаборатории изучения культурного наследия в Университете ИТМО. Область научных интересов: терагерцовая оптика и спектроскопия в искусстве, медицине и биологии; изучение культурного наследия; фемтосекунд-ная и нелинейная оптика, нанобиофотоника. Является членом SPIE, IEEE и Оптического общества им. Д.С. Рождественского. Преподавание курсов: биофотоника, биофизика и медицинская оптическая томография. Является заместителем руководителя отдела Life Science & Health Университета ИТМО, членом Pho-tonics21 - The European Technology Platform и членом научной школы Российской Федерации «Физика и астрономия». Лауреат индивидуальных грантов Санкт-Петербурга; премии Джорджа Сороса. E-mail: smolyanskaya@itmo. ru .
Перро Жан-Батист, 1987 года рождения, получил степень магистра в области авиационного технического обслуживания в институте Evering в 2011 году и степень магистра физики в Университете Бордо, Франция, в 2013 году. Получил степень доктора философии по направлению «Лазерная материя и науковедение» в 2018 году в лаборатории IMS университета Бордо в сотрудничестве с CEA Tech; диссертация на тему: «Терагерцевая визуализация со скоростью видеосигнала для промышленного неразрушающего контроля (full field terahertz imaging at video rate for industrial non-destructive testing)». В настоящее время является генеральным директором и соучредителем компании Optikan - стартапа, созданного в 2021 году, миссией которого является разработка, производство и коммерциализация решений неразрушающего контроля в ТГц-диапазоне для промышленности. E-mail: _iean-baptiste.perraud@optikan.com .
Гийе Жан-Поль, 1984 года рождения, получил степень магистра и степень доктора философии в области терагерцовой визуализации в 2007 и 2010 годах соответственно, лаборатория IES, Монпелье, Франция. Тема первой исследовательской работы - терагерцовая микроскопия ближнего поля с использованием компонентов на основе электроники. До 2013 года занимался исследованиями в области неразрушающего контроля с использованием терагерцовых систем импульсного и непрерывного действия. В 2013 году стал адъюнкт-профессором лаборатории IMS (UMR CNRS 5218) Университета Бордо, где в настоящее время занимается терагерцовой визуализацией и его приложениями. E-mail: _ iean-paul.guillet@u-bordeaux.fr .
Моне Патрик, получил степень инженера по материаловедению в 1988 году и степень доктора философии по квантовым устройствам в 1992 году в Университете Лилля (Ecole Universitaire D'Ingénieur de Lille (EUDIL)). После этого присоединился к CNRS в Институте микроэлектроники и нанотехнологий (IEMN UMR 8520). В 2002 году переехал в Бордо, где занялся терагерцовой спектроскопией и визуализацией. Сейчас является старшим научным сотрудником IMS UMR CNRS 5218, Университет Бордо. В настоящее время интересуется методами равновесной и терагерцовой спектроскопии с временным разрешением; участвует в проектах 2D- и SD-терагерцовой визуализации, 3D компьютерной томографии для исследования объектов культурного наследия, аддитивных технологий. Также областью научных интересов являются диэлектрические метаматери-алы и их применение в микроволновом и терагерцовом диапазонах частот; биомедицинские приложения тера-герцовой спектроскопии и визуализации. Является соавтором более 210 статей. E-mail: patrick.mounaix@u-bordeaux.fr .
Петров Николай Владимирович, 1983 года рождения, д.ф.-м.н. C 2016 года руководит лабораторией цифровой и изобразительной голографии в Университете ИТМО. С 2021 года ведущий профессор. Является соавтором более 300 публикаций, из которых более 150 в системе цитирования Scopus. Является старшим членом международного сообщества OPTICA. Область научных интересов: прикладная голография, фемтосекундная оптика, ТГц-имаджинг. E-mail: n.petrov@niuitmo.ru .
ГРНТИ: 29.31.29
Поступила в редакцию 21 ноября 2022 г. Окончательный вариант - 11 июля 2023 г.
An algorithm of unordered wavefront propagation in terahertz phase retrieval with dense multiplane data acquisition
E.G. Tsiplakova 1, A. Chopard2'3, N.S. Balbekin 1, O.A. Smolyanskaya 1, J.-B. Perraud4, J.-P. Guillet2, P. Mounaix2, N.V. Petrov1
1ITMO University, 197101, Saint Petersburg, Russian Federation, Kronverksky Pr. 49;
2IMS Laboratory, University of Bordeaux, UMR CNRS 5218, 33405, Talence, France, 351 Cours de la Libération Bâtiment A31;
3 Lytid SAS, 92120, Montrouge, France, rue la Fontaine 8;
4 Optikan SAS 33000, Bordeaux, France, rue Permentade 64
Abstract
Iterative phase retrieval algorithms from multiple diffraction patterns in the terahertz (THz) frequency range are a promising tool of computational imaging capable of providing high spatial resolution of reconstructed phase images. One of the commonly used algorithms is SBMIR, which employs multiple intensity distributions of the diffraction object wavefield as input data. Compared with single-frame methods, the multi-plane approach allows for a faster convergence, but requires time-consuming data acquisition from a receiver positioned at a variety of distances from the object. Previously, we proposed a method for THz data acquisition in a single scan mode, which allows one to quickly obtain an exhaustive set of diffraction distributions. In this paper we evaluate an up-to-date phase retrieval algorithm based on the SBMIR/R-SBMIR method (which utilizes stochastic wavefront propagation) on the experimental data captured by a single-scan technique. Unlike a number of conventional phase retrieval algorithms, which may require a series of numerical experiments for determining optimal intensity distributions from a large dataset, disordered propagation of the estimation wavefront guarantees the high-contrast and high-resolution image reconstruction without pre-setting the parameters. It is shown that the package use of the single-scan technique with the subsequent data processing using the R-SBMIR algorithm has application potential for automation of the multi-plane phase retrieval in the THz range.
Keywords: phase retrieval, terahertz radiation, iterative algorithm, stochastic optimization, phase imaging, quantum cascade laser.
Citation: Tsiplakova EG, Chopard A, Balbekin NS, Smolyanskaya OA, Perraud GB, Guillet JP, Mounaix P, Petrov NV. An algorithm of unordered wavefront propagation in terahertz phase retrieval with dense multiplane data acquisition. Computer Optics 2023; 47(6): 901-912. DOI: 10.18287/2412-6179-CO-1253.
Acknowledgements: This work was supported by the Russian Foundation for Basic Research and Centre National de la Recherche Scientifique under grant No. 21-52-15035/21.
Authors' information
Elizaveta Georgievna Tsiplakova (b. 2000) graduated from ITMO University in 2022 with a bachelor's degree in Photonics and Optoinformatics. Currently she is the first-year student of a master's degree at ITMO University in the specialty Technical Physics. Research interests: phase retrieval algorithms, terahertz optics, digital holography. E-mail: e.tsiplakova@itmo.ru .
Adrien Chopard was born in France, in 1993. He received the Diplome d'Ingénieur in Optics and Photonics from the Institut d'Optique Graduate School (IOGS-Saclay), Palaiseau, France, in 2016 and the M.Sc. degree in Engineering specialization in Nanotechnology from the Royal Institute of Technology of Sweden, (KTH), Stockholm, Sweden. In 2021, He finalized its industrial Ph.D. degree in Physics on applicative developments using terahertz technologies, in collaboration between the IMS Laboratory (UMR CNRS 5218), University of Bordeaux, Talence, France and the firm Lytid, Paris, France, for which he is currently in charge of the development of applicative solutions. E-mail: adrien.chopard@u-bordeaux.fr .
Nikolay Sergeevich Balbekin (b. 1989), PhD (Phys.-Math.), head of academic laboratories at ITMO University, deputy head of the Digital and Display Holography Laboratory. He has co-authored 75 publications. Research interests: analog holography, digital holography, terahertz spectroscopy. E-mail: nbalbekin@niuitmo.ru .
Olga Alekseevna Smolyanskaya is an associate professor at ITMO University. She is also the head of Laboratory of Heritage Science at ITMO University. Her research interests include terahertz optics and spectroscopy in art, heritage science, medicine and biology, femtosecond and nonlinear optics, tissue optical clearing, and nanobi-ophotonics. She is a member of SPIE, IEEE, and D.S. Rozhdestvensky Optical Society. Teaching courses: Bio-photonics, Biophysics, and Medical Optical Tomography. Smolyanskaya O.A. is a vice-chief of Life Science & Health at ITMO University and a member of Photonics21 - The European Technology Platform, and a member of scientific school of Russian Federation «Physics and astronomy». Winner of Saint-Petersburg individual grants; George Soros award. E-mail: smolyanskaya@itmo.ru .
Jean Baptiste Perraud (b. 1987) received the M.S. degree in Aeronautical Maintenance Engineering from Evering Institute in 2011, and the M.S. degree in Physics from the University of Bordeaux, France, in 2013. He obtained his PhD "Laser Matter and Nanoscience" in 2018 at IMS laboratory of Bordeaux University in cooperation with CEA Tech. His work demonstrated full field terahertz imaging at video rate for industrial nondestructive testing (NDT). He is now CEO and co-founder of Optikan, a start-up created in 2021 that develops, manufactures, and commercializes terahertz non-destructive testing solutions for industry. E-mail: jean-baptiste.perraud@optikan.com .
Jean-Paul Guillet was born in Lyon, France, in 1984. He received the M.Sc. and Ph.D. degrees in Terahertz Imaging from the IES Lab, Montpellier, France, in 2007 and 2010, respectively. The topic of this first research work was terahertz near-field microscopy using electronics-based components. During his postdoctoral studies, he focused on non-destructive testing application using pulsed and continuous wave terahertz systems, until 2013. In 2013, he was an Associate Professor with the IMS (UMR CNRS 5218), University of Bordeaux, where he is currently focusing on terahertz imaging and applications. His research is focused on imaging and tomography, using femtosecond pulses, FMCW approach, and a terahertz camera, including the development of terahertz subsystems. E-mail: jean-paul.guillet@u-bordeaux.fr .
Patrick Mounaix received the degree of Engineer of Material Science in 1988 from the Ecole Universitaire D'Ingénieur de Lille (EUDIL) and PhD on quantum devices in 1992 from the University of Lille. Then he joined CNRS at the Institute for Microelectronic and Nanotechnology (IEMN UMR 8520). He moved to Bordeaux in 2002 where he developed Terahertz Spectroscopy and Imaging Activities. Now he is senior researcher at IMS UMR CNRS 5218, Bordeaux University and he's currently interested in equilibrium and time-resolved terahertz spectroscopy techniques for solid and liquid materials. He's also involved in 2D & 3D terahertz imaging projects chiefly on 3D Computed Tomography for Art science, additive manufacturing and NDE industrial applications. He works also on all dielectric metamaterials and their applications at microwave and terahertz frequency range. Recently, he also works on biomedical applications of terahertz spectroscopy and imaging. Dr Mounaix is coau-thored more than 210 papers. E-mail: patrick.mounaix@u-bordeaux.fr .
Nikolay Vladimirovich Petrov (b. 1983), Dr. Sci. (Phys.-Math.). Since 2016 he is a head of the Digital and Display Holography laboratory in ITMO University. Leading Professor since 2021. His research interests are currently focused on holography, femtosecond optics, and terahertz imaging. He has co-authored more than 150 publications in the Scopus database. He is a senior member of OPTICA. E-mail: n.petrov@niuitmo.ru .
Received November 21, 2022. The final version - July 11, 2023.
