ЖЭТФ, 2022, том 162, вып. 6 (12), стр. 811-822

ФОРМИРОВАНИЕ ИЗОБРАЖЕНИЙ В ФАНТОМНОЙ

ВОЛОКОННОЙ ЭНДОСКОПИИ МЕТОДОМ РЕДУКЦИИ

ИЗМЕРЕНИЙ

Д. А. Балакин^*, Д. П. Агапов^**, П. П. Гостев^***, С. А. Магницкий^****, Д. Н. Фроловцев^†,

А. С. Чиркин^‡

Московский государственный университет им. М. В. Ломоносова, физический факультет

119991, Москва, Россия

Поступила в редакцию 15 июня 2022 г.,

после переработки 3 августа 2022 г.

Принята к публикации 7 августа 2022 г.

Для формирования фантомных изображений (ФИ) в волоконно-оптических системах применен метод

редукции измерений. ФИ восстанавливались путем численного моделирования полного цикла получения

ФИ в базовой схеме волоконного многомодового эндоскопа, начиная с распространения случайно моду-

лированного света по многомодовому волокну и заканчивая расчетом ФИ с помощью численных методов

однопиксельной визуализации. Для получения высококачественных восстановленных ФИ привлекалась

априорная информация об исследуемом объекте. Показано, что при числе используемых шаблонов осве-

щения объекта, существенно меньшем числа пикселей восстанавливаемого изображения, погрешность

восстановления объекта может быть уменьшена, несмотря на ограничение числа мод, распространя-

ющихся по оптическому волокну и формирующих некогерентное излучение. Для сравнения для фор-

мирования ФИ были использованы несколько вариантов метода сжатых измерений и традиционный

корреляционный метод. Продемонстрировано, что предложенный вариант метода редукции измерения в

применении к многомодовым волокнам показывает сравнимые, а часто лучшие результаты, чем тради-

ционный корреляционный метод и метод сжатых измерений.

DOI: 10.31857/S004445102212001X

основе получаемой косвенной информации об объ-

EDN: LBWCFT

екте с использованием фантомного принципа синте-

за изображений. Это требует не только разработки

специальных физических методов управления свой-

1. ВВЕДЕНИЕ

ствами света, но и разработки специальных матема-

тических алгоритмов.

В последние несколько лет наблюдается повы-

шенный интерес к новому направлению оптики

волоконной фантомной оптике (ВФО), которая на-

Перспективным направлением может стать при-

ходится на стыке волоконной, квантовой, статисти-

менение методов ВФО для получения изображе-

ческой оптики, оптики фантомных изображений и

ний в эндоскопических системах, где наблюдение за

интеллектуальных систем формирования и обработ-

объектом происходит через оптическое волокно или

ки изображений. Ключевым моментом в ВФО яв-

жгут волокон. Интерес прежде всего связан с воз-

ляется приготовление пространственно структури-

можностью поместить основные оптические элемен-

рованного света и формирование изображений на

ты вне зонда, что позволяет уменьшить его разме-

ры и тем самым улучшить качество инвазивной ди-

^* E-mail: balakin_d_a@physics.msu.ru

агностики. Отметим также, что принцип фантом-

^** E-mail: dimaagapov@mail.ru

ных изображений (ФИ) позволяет исключить влия-

^*** E-mail: fongostev@gmail.com

ние флуктуаций и различных неоднородностей сре-

^**** E-mail: sergeymagnitskiy@gmail.com

ды на трассе распространения излучения на восста-

^† E-mail: dfrolovtsev@gmail.com

^‡ E-mail: aschirkin@physics.msu.ru

новленное изображение.

811

Д. А. Балакин, Д. П. Агапов, П. П. Гостев и др.

ЖЭТФ, том 162, вып. 6 (12), 2022

Изначально, в отличие от формирования обыч-

функции.

ных изображений, при детектировании фантомного

Типичная модель измерения при обработке изоб-

изображения [1-3] для получения информации об

ражений имеет вид

объекте использовалась корреляция фотонов двух

пучков, один из которых (в объектном канале) осве-

ξ = Af + ν,

(1)

щает объект, после чего регистрируется простран-

где вектор ξ

∈ X результат измерения (по-

ственно неразрешающим (однопиксельным) датчи-

следовательность показаний однопиксельного дат-

ком, а второй (в восстанавливающем канале) реги-

чика), неизвестный вектор f

∈ F характеризует

стрируется многопиксельной камерой. В этом ва-

степень прозрачности объекта исследования, опе-

рианте ФИ, который будем в дальнейшем назы-

ратор A моделирует измерительный преобразова-

вать традиционным, для формирования изобра-

тель (в рассматриваемой схеме оператор A опреде-

жений предлагалось вычислять корреляционную

лен развернутыми в строки шаблонами освещения),

функцию второго порядка сигнала однопиксельного

ν ∈ X погрешность измерений, F,X конечно-

датчика в объектном канале и сигнала, соответству-

мерные евклидовы пространства. В этих обозначе-

ющего каждому отдельному пикселю в восстанав-

ниях традиционный алгоритм восстановления фан-

ливающем канале. Позднее появился вычислитель-

томных изображений описывается формулой

ный вариант ФИ, в котором объект освещается слу-

чайными шаблонами [4], формируемыми простран-

(

)(

)

∑

ственным модулятором света (см. [5] и цитируемую

ξ_j -

ξ_k

a_j -

a_l

там литературу), или цифровым микрозеркальным

dim X

j=1

k=1

l=1

устройством, а регистрация фотонов в восстанав-

(2)

ливающем канале заменяется ее численным моде-

где a_j j-я строка матрицы оператора A, т. е. раз-

лированием, что позволяет обойтись без многопик-

вернутый в строку j-й шаблон освещения, ξ_j пока-

сельной камеры. Позднее такой вариант ФИ и про-

зания однопиксельного датчика при освещении объ-

изводные от него алгоритмы начали называть од-

екта j-м шаблоном.

нопиксельной визуализацией. Если формирование

В современных системах ФИ с применением оп-

шаблонов освещения воспроизводимо, вместо чис-

тических волокон используется несколько основных

ленного моделирования можно предварительно за-

подходов. Так, в работе [7] обсуждается идея исполь-

регистрировать освещение в плоскости объекта, вре-

зования жгута одиночных волокон. В [8,9] для об-

менно заменив его камерой (на этапе формирования

лучения объекта применяется одиночное многомо-

изображения объекта многопиксельная камера по-

довое волокно, на входном торце которого распола-

прежнему не используется) [6]. Такой алгоритм пло-

гается жгут одиночных волокон. Еще один тип си-

хо подходит для систем реального времени, пред-

стем в ВФО включает в себя многопиксельные моду-

назначенных для генерации изображений объекта с

ляторы света и одиночные многомодовые волокна.

частотой 10 Гц и более, так как для построения од-

Этот подход экспериментально реализован в рабо-

ного изображения объекта требуется накопить боль-

тах [6,10,11]. В таких системах формирование изоб-

шое количество данных для расчета корреляцион-

ражения может быть осуществлено с помощью раз-

ной функции.

личных математических алгоритмов (см., например,

На рис.

приведена схема формирования

[12]). Результаты последних исследований показы-

фантомных изображений при облучении объекта

вают, что наиболее перспективными для создания

светом, прошедшим через оптическое волок-

волоконно-оптических фантомных эндоскопов явля-

но. Лазерное излучение проходит через фазо-

ются системы, использующие одиночные многомо-

вый/амплитудный пространственный модулятор

довые волокна. Такие системы позволяют создать

света, в дальней зоне которого профиль интенсив-

миниатюрные эндоскопы, обладающие минималь-

ности пучка представляет из себя спекл-картину.

ной степенью инвазивности, т. е. позволяют прово-

С помощью микрообъектива излучение заводится

дить вмешательства с минимальной травматизаци-

в многомодовое волокно. После распространения

ей окружающих тканей.

в волокне свет облучает объект исследования.

В принципе, ФИ могут быть синтезированы с

Излучение, прошедшее через объект, собирается

применением различных структур света, освещаю-

однопиксельным датчиком, сигнал с которого по-

щих объект (шаблонов освещения). На сегодняшний

ступает на персональный компьютер (ПК на рис. 1)

день предпочтение отдается статистически незави-

для вычисления пространственной корреляционной

симым шаблонам со случайной структурой. В дан-

812

ЖЭТФ, том 162, вып. 6 (12), 2022

Формирование изображений в фантомной волоконной эндоскопии

как новые строки матрицы оператора A оказывают-

ся линейно зависимыми от старых. Число мод опре-

деляет верхний предел ранга оператора A. В связи

с этим заметим, что хотя волоконные фантомные

изображения (ВФИ) могут быть получены в схемах

с различными конфигурациями волокон [8-10], наи-

более перспективными являются одиночные много-

модовые волокна. Это связано с тем, что соотноше-

Рис. 1. Схема для реализации фантомной визуализации с

ние размера волокна и количества пространствен-

помощью многомодового волокна

ных мод в многомодовых волокнах лучше, чем в

жгутах одномодовых и многомодовых волокон.

ной работе мы ограничимся рассмотрением форми-

Для повышения качества восстановления векто-

рования ФИ только таким светом.

ра f требуется привлечь дополнительную информа-

На рис.

представлена базовая схема фор-

цию о нем. Эта информация, как правило, ограни-

мирования фантомных изображений при облуче-

чивает класс объектов, изображения которых мо-

нии объекта светом со случайной пространствен-

гут быть точно восстановлены. Иными словами, по-

ной структурой, прошедшим через многомодовое

грешность обработки уменьшается при восстановле-

оптическое волокно. В такой схеме с теоретиче-

нии изображений ¾более правдоподобных¿ объектов

ской точки зрения учитываются практически все ас-

в ущерб восстановлению изображений ¾менее прав-

пекты формирования ФИ в волоконно-оптических

доподобных¿. Например,

схемах, основой которых являются многомодовые

• кусочно-постоянный вид оптических характе-

волокна. Лазерное излучение проходит через фа-

ристик объекта, если он состоит из оптически одно-

зовый/амплитудный пространственный модулятор

родных материалов;

света, в дальней зоне которого профиль интенсив-

• близость оптических характеристик объекта в

ности пучка представляет из себя спекл-картину.

близких его точках.

С помощью микрообъектива излучение заводится

В данной работе рассмотрен ряд математических

в многомодовое волокно. После распространения в

методов, позволяющих сформировать ФИ на осно-

волокне свет облучает объект исследования. Излу-

ве данных, получаемых в схеме, изображенной на

чение, прошедшее через объект, собирается одно-

рис. 1, а именно, традиционный метод, несколько

пиксельным датчиком, сигнал с которого поступа-

вариантов метода сжатых измерений и предложен-

ет на ПК. Для получения ФИ используются специ-

ный в работе вариант метода редукции измерения,

альные математические алгоритмы, входными дан-

в котором используется дополнительная информа-

ными которых являются пространственное распре-

ция о равенстве нулю многих компонент восстанав-

деление интенсивности света в плоскости объекта и

ливаемого изображения в собственном базисе моде-

интегральная интенсивность прошедшего через объ-

ли интерпретации измерения. Показано, что исполь-

ект света, собранная однопиксельным детектором.

зование дополнительной информации об объекте ис-

Заметим, что исторически первый традиционный

следования позволяет значительно уменьшить чис-

алгоритм, основанный на вычислении ковариацион-

ло шаблонов освещения, которыми требуется осве-

ной функции, и все последующие, включая разраба-

тить объект для формирования фантомного изоб-

тываемые в настоящее время, алгоритмы являются

ражения, и, следовательно, ускорить регистрацию.

статистическими в том смысле, что требуют доволь-

но значительного числа измерений, полученных при

различных пространственных профилях освещения

2. МЕТОД СЖАТЫХ ИЗМЕРЕНИЙ

объекта (шаблонов освещения).

2.1. Общая идея метода сжатых измерений

В ВФО проблема конечного числа шаблонов

освещения обостряется из-за того, что в оптическом

Как правило, для уменьшения числа шабло-

волокне от модулятора до объекта распространяет-

нов освещения, которыми требуется осветить объ-

ся ограниченное число мод. Поэтому даже увеличе-

ект для формирования фантомного изображения,

ние числа шаблонов освещения выше определенного

используется метод сжатых измерений [13-18]. В

предела позволяет лишь более точно оценить состав-

нем оценка вектора f определена решением задачи

ляющие изображения, информация о которых уже

минимизации по g (здесь g вектор, характери-

была получена из предшествующих измерений, так

зующий степень прозрачности некоторого объекта,

813

Д. А. Балакин, Д. П. Агапов, П. П. Гостев и др.

ЖЭТФ, том 162, вып. 6 (12), 2022

выступающего как пробный объект в смысле срав-

• норма L¹ в заданном базисе,

нения реального результата измерений и результата

Ω(g) = ∥T g∥₁

(4)

измерений, который был бы получен при измерении

пробного объекта) функционала

(основной используемый вариант), см. [6, 24];

Δ(Ag, ξ) + αΩ(g) ∼ min,

(3)

• отношение квадратов норм L¹ и L²,

g∈F

Ω(g) = ∥T g∥²¹/∥T g∥², см. [22].

где функционал Δ характеризует отклонение ожи-

даемого (при объекте, моделируемом вектором g)

2.3. Характеристики регулярности,

результата измерения Ag и фактически полученно-

основанные на конечных разностях

го ξ, а регуляризующий функционал Ω характеризу-

Существуют также варианты метода сжатых из-

ет регулярность объекта (как правило гладкость

мерений, в которых в формуле (3) используются

его изображения как функции координат).

другие характеристики регулярности изображения:

Основные используемые варианты выбора функ-

• полная кривизна [18, 25]

ционала Δ:

• функционал наименьших квадратов [6, 17-21]

∑

(



)

Δ(Ag, ξ) = ∥Ag - ξ∥²;

Ω(g) =

(D^2xg)_i

+

(D^2yg)_i



;

• логарифм функции правдоподобия [22]. При

i=1

пуассоновской статистике фотоотсчетов

• квадратичная кривизна [25]

∑

(



)

Δ(Ag, ξ) =

(ξ_j ln((Ag)_j + ǫ) -

²

Ω(g) =

(D^2xg)_i

² +

(D^2yg)_i

;

(5)

j=1

i=1

- ((Ag)_j + ǫ) - ln ξ_j !) ,

• сглаженная (при β = 0) изотропная (в смысле

где ǫ среднее число фотоотсчетов, вызванных тем-

чувствительности к ориентации координатных осей

новыми шумами, ! обозначает факториал.

в плоскости изображения) полная вариация [26]

Выбор функционала Δ, как правило, опреде-

ляется доступной информацией о модели измере-

∑

√

Ω(g) =

β² + |D^1xg|²ⁱ + |D^1yg|²ⁱ;

(6)

ния (известна ли статистика фотоотсчетов и распре-

i=1

деление темновых шумов, определяющие функцию

правдоподобия).

• анизотропная полная вариация [18]

∑

(

)

2.2. Использование разреженности в базисе

Ω(g) =

|D^1xg|_i + |D^1yg|_i

заданного преобразования

(7)

i=1

Коэффициенты пропускания реальных объек-

= ∥D^1xg∥₁ + ∥D^1yg∥₁;

тов в близких точках часто имеют схожие вели-

• альтернативный вариант анизотропной полной

чины. Это утверждение может быть математиче-

вариации

ски сформулировано как разреженность вектора f

Ω(g) = ∥D^1xg∥ + ∥D^1yg∥.

(8)

в некотором (выбираемом априори, в том числе

на основе информации о характерном виде объек-

Здесь (D^kx,yg)_i

конечная разность k-го порядка

та) базисе, заданном преобразованием T [6, 18-22].

изображения g в направлении оси x (y) в позиции

В качестве этого базиса, как правило, применя-

i-го пикселя.

ется базис дискретного косинусного преобразова-

ния (DCT) [22], см. сравнение ряда преобразова-

2.4. Свойства метода сжатых измерений.

ний (тождественного преобразования, дискретного

Выбор параметра регуляризации

вейвлет-преобразования, DCT) в работе [23]. Основ-

ное требование к преобразованию как можно бо-

Если восстанавливаемое изображение содержит

лее быстрое убывание упорядоченных по невозрас-

компоненты, не влияющие на полученное изображе-

танию модулей координат вектора f. На этом спосо-

ние и восстанавливаемые за счет информации о раз-

бе основаны следующие варианты выбора функци-

реженности, для обработки применяется вариант

онала Ω в формуле (3):

регуляризованного метода наименьших квадратов с

814

ЖЭТФ, том 162, вып. 6 (12), 2022

Формирование изображений в фантомной волоконной эндоскопии

поиском минимума регуляризующего функционала

В дополнение к таким широко известным спосо-

при требовании воспроизведения полученных изоб-

бам выбора значения параметра регуляризации α,

ражений [27], формально получающийся как предел

как принцип невязки [29] (к которому близка реко-

решения задачи (3) при α → +0. Для этого случая

мендация [18] использовать значение α, при котором

доказана [14] оптимальность выбора нормы L¹ в вы-

для полученной оценки

f отличие прогнозируемо-

бранном базисе как регуляризующего функционала

го результата от наблюденного ∥Af - ξ∥² ≈ dimX),

при отсутствии шума и определенных ограничени-

метод L-кривой [30, 31], метод обобщенной кроссва-

ях, накладываемых на модель измерения. Требует-

лидации

[32], в

[25] предложен способ, в кото-

ся, чтобы матрица оператора AT обладала ограни-

ром несколько раз моделируется получение изоб-

ченной изометрией, а именно, если для любого нату-

ражения объекта, описываемого гауссово размытой

рального s ≤ dim F существует такая константа изо-

оценкой максимального правдоподобия, вычисляет-

метрии δ_s, что (1 - δ_s)∥z∥² ≤ ∥Az∥² ≤ (1 + δ_s)∥z∥²

ся среднее уменьшение значения логарифма функ-

для любого вектора z ∈ F, имеющего не более s

ции правдоподобия при переходе от размытой оцен-

ненулевых компонент, и δ_s + δ_2s + δ_3s ≤ 1, то точ-

ки к результату моделирования и ищется значение

но восстановим любой вектор, имеющий не более s

α, при котором результат обработки настолько же

ненулевых компонент. Смысл этого условия приме-

уменьшает значение логарифма функции правдопо-

нительно к схеме на рис. 1 объекты, распреде-

добия по сравнению с оценкой максимального прав-

ления прозрачностей которых представимы с помо-

доподобия.

щью малого числа векторов из выбранного базиса,

В силу его схожести с рассматриваемым далее

при измерении дают достаточно сильно различаю-

вариантом метода редукции измерения также от-

щиеся результаты.

метим метод восстановления изображения при по-

мощи усеченного сингулярного разложения [8, 33],

Кроме указанной причины, выбор функционала

обычно не относимый к методам сжатых измере-

Ω часто определяется вычислительными соображе-

ний, но формально получающийся по формуле (3)

ниями. Приведенные выше функционалы, в кото-

при Ω(g) = max{0 ≤ k ≤ rk A|(I - P_k)g = 0},

рых используется не норма L¹ изображения в неко-

где P_k проектор на линейную оболочку первых k

тором базисе, а конечные разности изображения,

упорядоченных по невозрастанию соответствующих

позволяют решать задачу (3) локальными мето-

им сингулярных чисел правых сингулярных векто-

дами (в частности, стохастической оптимизацией),

ров оператора A, rk ранг оператора или матри-

поскольку изменение пикселя восстанавливаемого

цы. В нем восстановленное изображение формиру-

изображения изменяет лишь ограниченное число

ется действием на ξ оператора, псевдообратного к A

слагаемых в регуляризующем функционале, соот-

после усечения до нуля k наименьших сингулярных

ветствующих этому пикселю и его соседям [18,25].

чисел оператора A.

Использование сглаженной изотропной полной ва-

В заключение раздела отметим возможность

риации (6) вместо несглаженной упрощает вычисле-

комбинирования различных регуляризующих функ-

ния, поскольку сглаженная изотропная полная ва-

ционалов, как это сделано, например, в работе [34],

риация дифференцируема, а несглаженная нет.

где использована линейная комбинация нормы L²,

Выбор в качестве функционала Δ логарифма

нормы L¹ и полной вариации. В этом случае требу-

функции правдоподобия приводит к формальной

ются более сложные (из-за необходимости миними-

эквивалентности соответствующего варианта мето-

зации или нахождения корня функции нескольких

да сжатых измерений и метода апостериорного мак-

параметров, а не одного) способы выбора значений

симума, использованного, например, в [28], в ко-

регуляризующих параметров, обобщающие принцип

тором оценивание выполняется максимизацией вы-

невязки [35,36], метод L-кривой [37], метод обобщен-

численного по правилу Байеса апостериорного рас-

ной кроссвалидации [38], или основанные на исполь-

пределения прозрачности, где априорное распреде-

зовании тестовых измерений [39]. Также применение

ление случайной прозрачности объекта исследова-

этих методов затрудняется тем, что, как правило, их

ния, например, пропорционально экспоненте отри-

свойства проанализированы для регуляризующего

цательной полной вариации изображения [28]. Со-

функционала суммы слагаемых вида ∥Qg∥² для

ответствующее заданному регуляризующему функ-

заданного набора операторов Q, в то время как в ме-

ционалу априорное распределение при этом может

тоде сжатых измерений обычно используются более

не быть распределением вероятностей ввиду ненор-

сложные регуляризующие функционалы, указанные

мируемости.

выше (исключение функционал (8)).

815

Д. А. Балакин, Д. П. Агапов, П. П. Гостев и др.

ЖЭТФ, том 162, вып. 6 (12), 2022

(а)

(б)

(в)

(г)

(а)

(б)

(в)

(г)

(д)

(е)

(ж)

(з)

(д)

(е)

(ж)

(з)

(и)

Рис.

2. Компьютерное моделирование восстановления

Рис.

3. Компьютерное моделирование восстановле-

фантомного изображения объекта исследования (а),

ния фантомного изображения объекта исследования

109.375 × 109.375 мкм, 128 × 128 пикселей, по результа-

(а),109.375×109.375 мкм, 128×128 пикселей, по результа-

там освещения объекта 1024 бинарными (на входе опто-

там освещения объекта 1024 бинарными (на входе оптово-

волокна) шаблонами освещения: (б) традиционным мето-

локна) шаблонами освещения: (б) традиционным методом

дом (2), методами сжатых измерений при выборе в ка-

(2), методами сжатых измерений при выборе в качестве ре-

честве регуляризующего функционала нормы L¹ в базисе

гуляризующего функционала нормы L¹ в базисе (в) DCT

(в) DCT и (г) преобразования Хаара, см. формулу (4) в

и (г) преобразования Хаара, см. формулу (4) в разд. 2.2,

разд. 2.2, (д) квадратичной кривизны, см. формулу (5) в

(д) квадратичной кривизны, см. формулу (5) в разд. 2.3,

разд. 2.3, (е) анизотропной полной вариации (7) и (ж) ее

(е) анизотропной полной вариации (7) и (ж) ее альтерна-

альтернативного определения (8), см. разд. 2.3, (з) линей-

тивного определения (8), см. разд. 2.3, (з) линейной ре-

ной редукции измерения, см. разд. 3.1, и (и) нового алго-

дукции измерения, см. разд. 3.1, и (и) нового алгоритма

ритма редукции измерения, см. разд. 3.2

редукции измерения, см. разд. 3.2

3. МАТЕМАТИЧЕСКИЙ МЕТОД РЕДУКЦИИ

В рассматриваемом случае, когда процесс изме-

ИЗМЕРЕНИЯ

рения не возмущает объект и восстанавливается

именно изображение объекта, U

= I. Примером

Другой вариант восстановления фантомного

ситуации, когда это не так, служит восстановле-

изображения дает математический метод редукции

ние изображения отдельного слоя трехмерного объ-

измерений

[40], в котором ищется преобразова-

екта. Если погрешность оценивания h(R, U), обу-

ние из некоторого класса (например, линейных

словленная линейным преобразованием результата

преобразований), минимизирующее погрешность

измерения оператором R, определена как средне-

использования результата преобразования как

квадратичная максимальная по всем значениям f

оценки изображения объекта исследования.

h(R, U) = sup E∥Rξ - Uf∥², то [40, §5.1] минимизи-

f∈F

рующий ее оператор линейной редукции, если вы-

3.1. Линейная несмещенная редукция

полнено условие

измерения

U (I - A^-A) = 0,

(9)

Пусть в модели измерения (1) f априори про-

извольный вектор, погрешность ν случайный век-

равен

R_∗ = U(Σ^-1/2νA)^-Σ^-1/2ν,

(10)

тор, Eν = 0 и известен невырожденный корреляци-

онный оператор шума Σ_ν : ∀x ∈ X Σ_ν x = Eν(x, ν), а

где^- обозначает псевдообращение. Линейной редук-

измерение на идеальном измерительном преобразо-

ции соответствует минимальная погрешность

вателе, результат которого интересует исследовате-

(

)

ля, моделируется линейным оператором U : F → U.

h(R_∗, U) = tr

U (A^∗Σ^-1νA)^-U^∗

(11)

816

ЖЭТФ, том 162, вып. 6 (12), 2022

Формирование изображений в фантомной волоконной эндоскопии

(а)

(б)

(в)

(г)

(а)

(б)

(в)

(г)

(д)

(е)

(ж)

(з)

(д)

(е)

(ж)

(з)

(и)

Рис.

4. Компьютерное моделирование восстановления

Рис.

5. Компьютерное моделирование восстановления

фантомного изображения объекта исследования (а),

109.375 × 109.375 мкм, 128 × 128 пикселей, по результа-

там освещения объекта 1024 бинарными (на входе оптово-

локна) шаблонами освещения: (б) традиционным методом

(2), методами сжатых измерений при выборе в качестве ре-

гуляризующего функционала нормы L¹ в базисе (в) DCT

и (г) преобразования Хаара, см. формулу (4) в разд. 2.2,

(д) квадратичной кривизны, см. формулу (5) в разд. 2.3,

(е) анизотропной полной вариации (7) и (ж) ее альтерна-

тивного определения (8), см. разд. 2.3, (з) линейной ре-

дукции измерения, см. разд. 3.1, и (и) нового алгоритма

редукции измерения, см. разд. 3.2

Если же условие (9) не выполнено, то априори до-

как собственный базис этого оператора. Собствен-

пускает оценивание с конечной погрешностью (рав-

ный базис обладает экстремальным свойством, со-

ной (11)) лишь составляющая UA^-Af. Условие (9)

гласно которому с. к. погрешность

означает, что изменения объекта, не отражающие-

∑

ся на результатах измерений, не должны влиять на

h(Π_sR_∗, Π_sU) =

σ^2j =

восстанавливаемое изображение.

j=1

= inf{h(ΠR_∗, ΠU)| rk Π ≥ s},

(12)

3.2. Использование дополнительной

информации об объекте исследования

s = 1,...,dimU,

Пусть задача редукции измерения разрешима

где Π_s ортогональный проектор на линейную обо-

при заданном U и требуется, чтобы среднеквадра-

лочку L(e₁, . . . , e_s) векторов e₁, . . . , e_s, Π ортого-

тичная (с. к.) погрешность h(R_∗, U) оценивания Uf

нальный проектор. Согласно (12), с. к. погрешность,

не превосходила заданную величину. Для ответа на

сопутствующая оцениванию ортогональной проек-

этот вопрос в [41, гл. 2, §3] введено понятие соб-

ции сигнала идеального измерительного преобразо-

ственного базиса модели интерпретации измерения.

вателя на линейное подпространство размерности

Например, в рассмотренном выше случае её орто-

s или более, минимальна, если это подпростран-

нормированный собственный базис e_j, j = 1, 2, . . . ,

ство есть линейная оболочка первых s базисных век-

определяется путем решения задачи [40, гл. 8, вве-

торов, упорядоченных по неубыванию соответству-

дение и §8.1] на собственные значения оператора

ющих собственных значений. Таким образом, соб-

U (A^∗Σ^-1νA)^-1U^∗e_j

= σ^2je_j, σ21 ≤ σ22 ≤ ..., т.е.

ственный базис оптимально (в указанном смысле)

817

Д. А. Балакин, Д. П. Агапов, П. П. Гостев и др.

ЖЭТФ, том 162, вып. 6 (12), 2022

(а)

(б)

(в)

(г)

(а)

(б)

(в)

(г)

(д)

(е)

(ж)

(з)

(д)

(е)

(ж)

(з)

(и)

Рис.

6. Компьютерное моделирования восстановления

Рис.

7. Компьютерное моделирования восстановления

фантомного изображения объекта исследования (а),

109.375 × 109.375 мкм, 128 × 128 пикселей, по результатам

освещения объекта 1024 шаблонами освещения, соответ-

ствующими случайной модуляции фазы в каждом пикселе:

(б) традиционным методом (2), методами сжатых измере-

ний при выборе в качестве регуляризующего функционала

нормы L¹ в базисе (в) DCT и (г) преобразования Хаара,

см. формулу (4) в разд. 2.2, (д) квадратичной кривизны,

см. формулу (5) в разд. 2.3, (е) анизотропной полной ва-

риации (7) и (ж) ее альтернативного определения (8), см.

разд. 2.3, (з) линейной редукции измерения, см. разд. 3.1,

и (и) нового алгоритма редукции измерения, см. разд. 3.2

отделяет наиболее зашумленные компоненты фор-

ний характеризуется ее математическим ожиданием

мируемой оценки от наименее зашумленных.

и ковариационным оператором. Построенный кри-

Ранее в [42] предложен алгоритм редукции из-

терий существенно сложнее (13), поскольку требует

мерений при априорной информации, аналогичной

найти максимальное n ≤ dim U и бинарную матрицу

используемой в методе сжатых измерений, о бли-

C ранга n, удовлетворяющую условию

зости оптических характеристик объекта в близких

его точках и формализации этой информации разре-

((CT ΣR∗ξT^∗C^∗)^-1CT û, CT û) < nτ²,

женностью изображения как вектора в априори за-

данном базисе. Как и при линейной редукции, пред-

где ΣR∗ξ = U(A^∗Σ^ν1A)^-U^∗ ковариационный опе-

полагалось выполнение условия (9). Алгоритм осно-

ратор оценки R_∗ξ. Наконец, в [44] предложено ис-

ван на проверке статистических гипотез о равенстве

пользование в качестве базиса, в котором проверя-

компонент восстанавливаемого изображения нулю

ются гипотезы о компонентах, собственного базиса

при общем виде используемого для проверки гипо-

модели интерпретации измерения в силу его экстре-

тез критерия. Также в [42] показано его примене-

мального свойства (12) и показано, что оптималь-

ние для критерия вида (13), см. ниже, но в задан-

ный критерий, охарактеризованный в теореме 1 в

ном исследователем базисе, при этом внедиагональ-

[43], в этом базисе приобретает наиболее простой

ные компоненты матрицы корреляций формируемо-

вид (13) благодаря диагональности матрицы кова-

го изображения не учитывались. Позднее в теореме

риаций восстанавливаемого изображения в этом ба-

1 в [43] определен оптимальный критерий, учитыва-

зисе.

ющий корреляционные связи компонент формируе-

С учетом сделанных замечаний алгоритм редук-

мой оценки, в ситуации, когда погрешность измере-

ции приобретает следующий вид.

818

ЖЭТФ, том 162, вып. 6 (12), 2022

Формирование изображений в фантомной волоконной эндоскопии

(а)

(б)

(в)

(г)

(а)

(б)

(в)

(г)

(д)

(е)

(ж)

(з)

(д)

(е)

(ж)

(з)

(и)

Рис.

8. Компьютерное моделирования восстановления

Рис.

9. Компьютерное моделирования восстановления

фантомного изображения объекта исследования (а),

109.375 × 109.375 мкм, 128 × 128 пикселей, по результатам

освещения объекта 1024 шаблонами освещения, соответ-

ствующими случайной модуляции фазы в каждом пикселе:

(б) традиционным методом (2), методами сжатых измере-

ний при выборе в качестве регуляризующего функционала

нормы L¹ в базисе (в) DCT и (г) преобразования Хаара,

см. формулу (4) в разд. 2.2, (д) квадратичной кривизны,

см. формулу (5) в разд. 2.3, (е) анизотропной полной ва-

риации (7) и (ж) ее альтернативного определения (8), см.

разд. 2.3, (з) линейной редукции измерения, см. разд. 3.1,

и (и) нового алгоритма редукции измерения, см. разд. 3.2

1. Вычисление оценки û = R_∗ξ (10) линейной

Если выполнено условие (9), шаги 1-4 реша-

несмещённой редукции.

ют задачу восстановления изображения. Но инте-

2. Переход к собственному базису модели интер-

рес представляет противоположный случай, так как

претации измерения, в котором искомое значение

именно он реализуется в случае dim X < dim U (чис-

TUf интересующей исследователя характеристики

ло шаблонов освещения меньше числа пикселей).

объекта исследования считается разреженным, т.е.

Тогда по тем же причинам, что и при линейной ре-

содержит много нулевых компонент по сравнению с

дукции, после шага 4 восстанавливается лишь со-

его размерностью. Обозначим соответствующее пре-

ставляющая UA^-Af, и для восстановления состав-

образование T, а его результат T û.

ляющей U(I -A^-A)f требуется, как и в методе сжа-

3. Формирование оценки T û_thr по правилу

тых измерений, привлечь дополнительную инфор-

= 0, если гипотеза ¾(T û)_i = 0¿ принимает-

мацию. Воспользуемся в качестве такой информа-

ся, а именно, если выполнено условие

ции используемым в методе сжатых измерений регу-

ляризующим функционалом Ω, считая, что восста-

(T û)²ⁱ < τ(T ΣR∗ξT^∗)ii,

(13)

навливаемое изображение тем более правдоподобно,

чем меньше значение этого функционала. Таким об-

= (T û)_i, i

= 1, . . ., dimU, где

разом, завершающие шаги алгоритма имеют следу-

ΣR∗ξ = U(A^∗Σ^ν1A)^-U^∗ ковариационный опера-

ющий вид.

тор оценки R_∗ξ, τ ≥ 0 параметр критерия (13).

5. Вычисление наиболее соответствующей до-

Возвращение к исходному базису:

полнительной информации версии составляющей

= T^-1T û_thr.

U (I - A^-A)f (любым) решением задачи минимиза-

819

Д. А. Балакин, Д. П. Агапов, П. П. Гостев и др.

ЖЭТФ, том 162, вып. 6 (12), 2022

ции

модуляция осуществляется с помощью жидкокри-

Ω(û_thr + U(I - A^-A)g) ∼ min .

(14)

сталлического модулятора света.

g∈F

Изменение пространственного профиля интен-

Обозначим результат действия оператора

сивности при распространении света в многомодо-

U (I - A^-A) на ее решение v_thr.

вом волокне рассчитывалось с помощью библиоте-

6. Восстановление изображения как суммы

ки pyMMF (https://github.com/vongostev/pyMMF).

û_thr + v_thr.

При моделировании использованы следующие зна-

В непосредственно рассматриваемом случае

чения параметров: диаметр оптического волокна

U = I задача (14) упрощается до

31.25 мкм, длина волокна 50 см, отсутствие изги-

ба, аппроксимация пространственного распределе-

Ω(û_thr + v) ∼ min

(15)

v∈F:Av=0

ния показателя преломления параболическим про-

филем с n₁ = 1.4613 и a = 31.25 мкм, численная

Величина τ ≥ 0 есть параметр алгоритма. Он

апертура 0.275, длина волны 632.8 нм.

отражает приемлемый для исследователя компро-

На рис. 2-5 показаны результаты применения ме-

мисс между подавлением шума (чем больше τ, тем

тодов восстановления изображения, рассмотренных

больше подавляется шум) и искажением изображе-

в предыдущих разделах, при случайной погрешно-

ний, компоненты которых близки к нулю. Для выбо-

сти измерений с дисперсией 0.1 и бинарной ампли-

ра значения τ исследователь может смоделировать

тудной модуляции. Использованные в вариантах ме-

регистрацию тестового изображения, содержащего

тода сжатых измерений значения параметра регуля-

требуемые интересующие исследователя ¾тонкие¿

ризации выбраны методом L-кривой (см. разд. 2.4).

детали, и использовать максимальное значение па-

В случае метода редукции измерения при исполь-

раметра, при котором эти детали сохраняются, или,

зовании дополнительной информации об объекте в

аналогично [39], минимизировать среднюю погреш-

качестве функционала Ω в (14) и (15) использована

ность восстановления по набору тестовых данных.

анизотропная полная вариация (7), поскольку имен-

Шаг 3 можно также интерпретировать как замену

но этот функционал показал лучшие результаты по

исходного оператора U на такой, ядру которого при-

сравнению с другими вариантами метода сжатых

надлежат наиболее поражённые шумом компоненты

измерений. На рис. 3, 4 видно, как ограниченность

оценки.

диаметра оптоволокна приводит к ухудшению ка-

чества восстановленного изображения в целом и в

особенности у его краев.

4. СРАВНЕНИЕ МЕТОДОВ

Видно, что предлагаемый вариант метода редук-

Сравнение описанных выше алгоритмов было

ции измерения показывает сравнимые и часто луч-

проведено при моделировании схемы, изображен-

шие результаты, чем варианты метода сжатых изме-

ной на рис. 1. Возможны различные способы выбо-

рений. Примечательно, что несмотря на теоретиче-

ра шаблонов освещения, от детерминированных (на-

ское преимущество варианта полной вариации (7),

пример, шаблонов Адамара) до (псевдо)случайных.

вариант (8) на рис. 3-5 показал близкие результа-

В настоящей работе использованы два варианта по-

ты к результатам варианта (7), а в ряде случаев

следних.

результаты использования функционала (8) оказа-

В первом варианте использована псевдослучай-

лись лучше (7). Представляется, что это связано

ная бинарная амплитудная модуляция, т.е. на вход

с большим числом перепадов яркости на объектах,

оптического волокна подается псевдослучайное

изображенных на рис. 3-5 по сравнению с рис. 2, в

изображение, яркость каждого пикселя которого с

связи с чем ¾более мягкая¿ регуляризация функци-

одинаковой вероятностью принимает либо нулевое,

оналом (8) приводит к меньшей погрешности благо-

либо единичное значение. Такая модуляция может

даря меньшему сглаживанию этих перепадов.

быть реализована с помощью микрозеркального

На рис. 6-9 показаны результаты компьютерного

модулятора. Во втором варианте использована

моделирования при получении шаблонов освещения

псевдослучайная фазовая модуляция, т.е. на вход

случайной фазовой модуляцией. Сравнение с рис. 2-

оптического волокна подается псевдослучайное

5 показывает, что сделанные выводы верны и в этом

изображение, каждый пиксель которого имеет

случае.

единичную яркость, но вносится псевдослучайная

фазовая задержка от 0 до 2π. Обычно фазовая

820

ЖЭТФ, том 162, вып. 6 (12), 2022

Формирование изображений в фантомной волоконной эндоскопии

5. ЗАКЛЮЧЕНИЕ

ЛИТЕРАТУРА

А. В. Белинский, Д. Н. Клышко, ЖЭТФ 105, 487

В работе предложен и применен новый метод ре-

(1994) [A. V. Belinskii and D. N. Klyshko, JETP 78,

дукции измерений для формирования изображений

259 (1994)].

с помощью одиночных многомодовых волокон. По-

T. B. Pittman, Y. H. Shih, D. V. Strekalov, and

казано, что конечное число мод в волокне не явля-

A. V. Sergienko, Phys. Rev. A 52, R3429 (1995).

ется принципиальным препятствием для получения

изображений в ВФО. Проведенный в работе анализ

А. Гатти, Э. Брамбилла, М. Баке и др., в Кван-

товое изображение, под ред. М. И. Колобова,

ряда вариантов метода сжатых измерений и вариан-

А. С. Чиркина, Физматлит, Москва (2009), с. 96

та метода редукции измерения для формирования

[A. Gatti, E. Brambilla, M. Bache et al., in Quantum

изображений в ВФО показал осуществимость этой

Imaging, ed. by M. I. Kolobov, Springer (2007), p.

задачи при числе шаблонов освещения, существенно

79].

меньшем числа пикселей восстанавливаемого изоб-

ражения, несмотря на ограниченное число мод, спо-

J. H. Shapiro, Phys. Rev. A 78, 061802 (2008).

собных распространяться по оптическому волокну

Д. П. Агапов, И. А. Беловолов, П. П. Гостев

от модулятора до объекта. Выбор алгоритма вос-

и др., ЖЭТФ 162, 215 (2022) [D. P. Agapov,

становления ВФО оказывает значительное влияние

I. A. Belovolov, P. P. Gostev et al., JETP 135, 188

на погрешность результата, в связи с чем приведен

(2022)].

пример ситуации, когда в методе сжатых измерений

L. V. Amitonova and J. F. de Boer, Opt. Lett. 43,

использование в качестве регуляризующего функ-

5427 (2018).

ционала полной вариации с нормой L² приводит к

лучшим результатам, чем использование функцио-

C. Liu, J. Chen, J. Liu et al., Opt. Exp. 26, 10048

нала полной вариации с нормой L¹. Продемонстри-

(2018).

ровано, что предлагаемый вариант метода редук-

T. Fukui, Y. Kohno, R. Tang et al., J. Lightwave

ции измерения показывает сравнимые результаты с

Technol. 39, 839 (2021).

вариантами метода сжатых измерений. Кроме то-

T. Fukui, Y. Nakano, and T. Tanemura, J. Opt. Soc.

го, применение метода редукции приводит к менее

Amer. B 38, 379 (2021).

ярко выраженным артефактам в области, где опти-

ческое волокно слабо передает освещение. Ср., на-

10.

A. M. Caravaca-Aguirre and R. Piestun, Opt. Exp.

пример, рис. 5е и 5и, где метод сжатых измерений

25, 1656 (2017).

восстанавливает удаленные от центра участки как

11.

L. M. Xiang, Y. Li, J. Gao et al., Opt. Exp. 28, 13662

две области постоянной яркости, а метод редукции

(2020).

не восстанавливает резкую границу областей. При

этом для получения рис. 5и для формулирования

12.

D. Yang, M. Hao, G. Wu et al., Opt. Lasers Eng.

дополнительной информации использован такой же

149, 106827 (2022).

регуляризующий функционал, как и в варианте ме-

13.

D. L. Donoho, IEEE Trans. Inf. Theory 52, 1289

тода сжатых измерений, использованном для полу-

(2006).

чения рис. 5е. Необходимость использования такой

14.

E. J. Candes and T. Tao, IEEE Trans. Inf. Theory

дополнительной информации, отличающая предло-

52, 5406 (2006).

женный вариант метода редукции измерений от его

непосредственных предшественников, как указано в

15.

E. J. Candes and M. B. Wakin, IEEE Signal

разд. 3.2, обусловлена тем, что без нее в результа-

Processing Magazine 25, 21 (2008).

те редукции была бы восстановлена лишь составля-

16.

M. F. Duarte, M. A. Davenport, D. Takhar et al.,

ющая изображения, непосредственно влияющая на

IEEE Signal Processing Magazine 25, 83 (2008).

результат измерений.

17.

S. Han, H. Yu, X. Shen et al., Appl. Sci. 8, 1379

Полученные результаты будут полезны при раз-

(2018).

работке эндоскопических систем реального времени

на основе технологии ВФИ.

18.

G. M. Gibson, S. D. Johnson, and M. J. Padgett, Opt.

Exp. 28, 28190 (2020).

Финансирование. Исследование выполнено за

счет гранта Российского научного фонда (проект

19.

P. Zerom, K. W. C. Chan, J. C. Howell et al., Phys.

№ 21-12-00155).

Rev. A 84, 061804 (2011).

821

Д. А. Балакин, Д. П. Агапов, П. П. Гостев и др.

ЖЭТФ, том 162, вып. 6 (12), 2022

20.

W. Gong and S. Han, Phys. Lett. A 376, 1519 (2012).

32.

G. H. Golub, M. Heath, and G. Wahba,

Technometrics 21, 215 (1979).

21.

W. Gong and S. Han, Sci. Rep. 5, 9280 (2015).

33.

P. C. Hansen, BIT 27, 534 (1987).

22.

P. A. Morris, R. S. Aspden, J. E. C. Bell et al., Nat.

34.

M. Marquez, P. Meza, H. Arguello et al., Opt. Exp.

Commun. 6, 5913 (2015).

27, 17795 (2019).

23.

J. Du, W. Gong, and S. Han, Opt. Lett. 37, 1067

35.

S. Lu and S. V. Pereverzev, Numerische Math. 118,

(2012).

1 (2010).

24.

X. Shi, X. Huang, S. Nan et al., Laser Phys. Lett. 15,

36.

Z. Wang, J. Comp. Appl. Math. 236, 1815 (2012).

045204 (2018).

37.

M. Belge, M. E. Kilmer, and E. L. Miller, Inverse

25.

L. Mertens, M. Sonnleitner, J. Leach et al., Sci. Rep.

Problems 18, 1161 (2002).

7, 42164 (2017).

38.

C. Brezinski, M. Redivo-Zaglia, G. Rodriguez et al.,

26.

T. Jiang, W. Tan, X. Huang et al., J. Opt. 23, 075201

Numerische Math. 94, 203 (2003).

(2021).

39.

J. Chung and M. I. Espanol, Inverse Problems 33,

27.

O. Katz, Y. Bromberg, Y. Silberberg et al., Appl.

074004 (2017).

Phys. Lett. 95, 131110 (2009).

40.

Ю. П. Пытьев, Методы математического мо-

делирования измерительно-вычислительных си-

28.

S. H. Chan and Y. M. Lu, in IEEE Global Conference

стем, Физматлит, Москва (2012).

on Signal and Information Processing (GlobalSIP)

Atlanta, GA (2014), p. 312.

41.

Ю. П. Пытьев, Математические методы интер-

претации эксперимента, Высшая школа, Москва

29.

А. В. Гончарский, А. С. Леонов, А. Г. Ягола, Ж.

(1989).

вычислит. матем. и матем. физ. 13, 294 (1973)

[A. V. Goncharskii, A. S. Leonov, and A. G. Yagola,

42.

D. A. Balakin, A. V. Belinsky, and A. S. Chirkin,

USSR Comp. Math. Math. Phys. 13, 25 (1973)].

Quantum Inf. Process. 18(3) (2019).

30.

P. C. Hansen and D. P. O’Leary, SIAM J. Sci. Comp.

43.

D. A. Balakin and A. V. Belinsky, Quant. Inf.

14, 1487 (1993).

Process. 19, 316 (2020).

31.

C. L. Lawson and R. J. Hanson, Solving Least Squares

44.

D. A. Balakin and Yu. P. Pyt’ev, Pattern Recognition

Problems SIAM, Philadelphia, PA (1995).

and Image Analysis 31, 601 (2021).

822