Акустический журнал, 2022, T. 68, № 4, стр. 440-448

1. ВВЕДЕНИЕ

В настоящее время в научных исследованиях и в клинической практике широко применяются такие методы неинвазивной неповреждающей диагностики, как магнитно-резонансная и компьютерная томография. Огромное значение этих методов было оценено мировым научным сообществом присуждением их авторам Нобелевских премий по физиологии и медицине. К сожалению, использование в этих методах потенциально опасных для здоровья человека жестких полей и излучений накладывает серьезные ограничения на их применение. Значительным потенциалом, с точки зрения безопасности и получения надежной диагностически значимой информации, обладают оптическая когерентная и оптическая диффузионная томография. Однако из-за сильного поглощения и рассеивания световых волн в биологических тканях чисто оптические методы способны эффективно визуализировать структуры лишь в приповерхностных областях. В этой связи одним из самых перспективных и безопасных направлений медицинской диагностики считается подход, основанный на оптоакустическом (ОА) эффекте.

Оптоакустическая визуализация основана на явлении, при котором поглощение оптического излучения неоднородностями биоткани вызывает их нагрев с последующим тепловым расширением. Следствием теплового расширения оптических поглотителей является генерация ультразвуковых колебаний, которые регистрируются антенной, расположенной на поверхности исследуемого объекта. Поскольку распространение и формирование ОА-сигналов зависит от теплофизических, акустических и оптических свойств среды, то взаимосвязь между этими факторами позволяет использовать зарегистрированные ультразвуковые сигналы для количественной оценки свойств среды, путем решения обратной задачи оптоакустики.

Помимо биомедицинской безопасности важным преимуществом ОА-визуализации является то, что этот подход объединяет возможности высококонтрастного оптического поглощения и глубокого ультразвукового проникновения. При этом качество изображения, восстановленного по зарегистрированным на поверхности ткани ультразвуковым сигналам, в значительной степени определяется алгоритмом оптоакустической реконструкции. Искажения и артефакты в реконструированном изображении могут появиться как из-за специфики метода реконструкции, так и вследствие неизбежных в реальных исследованиях препятствий и шумов различной природы. Обычно ОА-изображения содержат интенсивный нестационарный фон с артефактами, структурно подобными сигналам, отношение сигнал/фон как правило невелико, а цифровое изображение имеет невысокое качество, малое число уровней квантования, пятенный характер и нечеткие границы. Поэтому задача устранения искажений и артефактов при восстановлении оптоакустических изображений весьма актуальна для эффективного использования метода в клинической практике и научных исследованиях [1].

Цель настоящей работы заключается в разработке и тестировании in silico алгоритмов удаления искажений и артефактов из реконструированных двумерных и трехмерных оптоакустических изображений. Вначале дается описание прямой и обратной задач оптоакустики, а также описание предложенного итеративного алгоритма коррекции восстановленных ОА-изображений. В следующем разделе подробнее описан алгоритм, числовой эксперимент и результаты моделирования реконструкции ОА-изображений для двумерного и трехмерного случаев. В Заключении приведены выводы и перспективы разработанного алгоритма.

2. ИТЕРАТИВНЫЙ АЛГОРИТМ ОПТОАКУСТИЧЕСКОЙ РЕКОНСТРУКЦИИ

Оптоакустическая визуализация основана на эффекте термоупругости, когда при поглощении неоднородностями среды оптического излучения происходит превращение оптической энергии в тепловую. При умеренной плотности выделившейся энергии в области поглощения не происходит фазовых превращений и следствием теплового оптического поглощения является генерация звуковых волн.

Прямая задача ОА-томографии заключается в том, чтобы определить поле давления p(r, t) по известному распределению тепловых источников H(r, t) (здесь r – пространственная координата точки, t – время), возбужденных кратковременным световым (лазерным) импульсом (рис. 1).

В акустически однородной бесконечной среде искомая пространственно-временна́я зависимость p(r, t) определяется уравнением [2]:

где c – скорость звука, β – коэффициент изобарического расширения, c_p – теплоемкость при постоянном давлении, приходящемся на единицу массы. Если тепловой источник H(r, t) представить в виде произведения поглощенной энергии и временно́й функции подсветки H(r,t) = Q(r)I(t), то в случае короткого импульса H(r,t) = Q(r)δ(t), где δ(t) – дельта-функция Дирака.

Начальное акустическое давление, возникшее за счет поглощения импульсного излучения оптическими неоднородностями в момент времени t = 0, можно представить в виде p₀(r) = ΓQ(r), где Γ – безразмерный коэффициент Грюнайзена, характеризующий эффективность ОА-превращения поглощаемого света в звук.

В этом случае решение прямой задачи задается выражением [2]:

где V ' – объем, в котором распределены ОА-источники.

Восстановление распределения начального акустического давления p₀(r, t = 0) = p₀(r) по сигналам давления p_s(r, t), зарегистрированным на поверхности S объема V ', составляет суть обратной задачи оптоакустики [3].

Методы восстановления распределения источников в среде можно разделить на несколько категорий: алгоритм Фурье-реконструкции, алгоритм обращения во времени, метод обратных проекций [4, 5].

Метод обратных проекций [6] является классическим способом ОА-реконструкции. Он может реализовываться либо в пространственно-временно́й, либо в Фурье-области для нескольких конфигураций детектирования в плоской [6], цилиндрической [6] или в сферической геометрии [7]. При этом решения будут справедливы только для замкнутой идеальной поверхности (приемником является каждая точка поверхности). Кроме того, обычно считается, что целевой объект расположен в бесконечно-однородной среде без дисперсии с постоянными скоростью звука, коэффициентами поглощения и плотностью.

В реальных ситуациях эти условия не выполняются, что приводит к искажениям реконструируемых изображений. Специфика этих искажений различна при использовании разных методов реконструкции, разном расположении приемников и разной геометрии реконструируемого объекта.

Для компенсации искажений и преодоления сильной зависимости качества изображения от указанных факторов нами был разработан метод итеративного уменьшения искажений ОА-изображений, основанный на теореме Банаха о неподвижной точке.

Пусть

где x – неизвестный входной образ, а f (⋅) – отображение (оператор), переводящее входное изображение x в результат решения обратной задачи оптоакустики $y$. При этом оператор f (⋅) представляет себой результат последовательного применения операторов f₁(⋅) и f₂(⋅): f (⋅) = f₂(f₁(⋅)), где f₁(⋅) – оператор, определяющий решение прямой задачи оптоакустики (уравнение (2)), а f₂(⋅) – оператор, задающий решение обратной задачи.

По теореме Банаха о неподвижной точке известно, что если f – сжимающее отображение множества F ⊂ (M, ρ) самого в себя, (M, ρ) – полное метрическое пространство, а F – замкнутое множество, то тогда существует, и притом ровно одна, неподвижная точка x* ∈ F отображения f (если точка неподвижная, то f (x*) = x*).

При этом под сжимающим отображением понимают такое отображение f: F → F ⊂ (M, ρ), что ∃α ∈ [0,1): ρ(f (x), f (y)) ≤ αρ(x,y), ∀x, y ∈ F.

В соответствии с этим определением для поставленной задачи реконструкции F = R^{m × n} представляет себой цифровое изображение, состоящее из m × n пикселей (для двумерного случая), а ρ(x, y) = ||x–y|| – евклидово расстояние. Можно показать, что F ⊂ (M,ρ) является полным метрическим пространством, для которого справедлива теореме Банаха о неподвижной точке [8]. В этом случае неподвижная точка x* ∈ F может быть найдена методом последовательных приближений x_n = f (x_{n – 1}), где нулевое приближение x₀ ∈F – произвольная точка метрического пространства.

Таким образом, уточненное решение обратной задачи оптоакустики можно получить, сформировав последовательность изображений {I_q} такую, что I_q = ℜ(I_{q – 1}), где ℜ(⋅) – оператор, задаваемый отображениями f₁(⋅) и f₂(⋅) и переводящий начальное распределение давления p₀(r) в его нулевое приближение $p_{0}^{{(0)}}(r) = {{I}_{0}}$. Другими словами, в качестве нулевого приближения принимается решение обратной задачи оптоакустики p₀(r), а возникшие при этом искажения удаляются в ходе итеративного процесса.

В предлагаемом подходе проверялись два алгоритма последовательных приближений: алгоритм, основанный на теореме Пикара [9], и метод, предложенный в работе [10] (соответственно формулы (4) и (5)):

Для количественной оценки качества восстановления образа использовались такие критерии, как отношение сигнал/шум $20{{\log }_{{10}}}\left( {{{{{p}_{0}}} \mathord{\left/ {\vphantom {{{{p}_{0}}} {\left( {{{p}_{0}} - p_{0}^{{(q)}}} \right)}}} \right. \kern-0em} {\left( {{{p}_{0}} - p_{0}^{{(q)}}} \right)}}} \right)$, относительная ошибка $E = {{\left\| {{{p}_{0}} - p_{0}^{{(q)}}} \right\|} \mathord{\left/ {\vphantom {{\left\| {{{p}_{0}} - p_{0}^{{(q)}}} \right\|} {\left\| {{{p}_{0}}} \right\|}}} \right. \kern-0em} {\left\| {{{p}_{0}}} \right\|}}$ (здесь $p_{0}^{{(q)}}$ – q-я итерация алгоритма) и индекс структурного подобия SSIM, характеризующий близость изображений X и Y по яркости, контрасту и структуре. Индекс структурного подобия (SSIM от англ. sructure similarity) считается неофициальным стандартом для оценки качества изображений при наличии эталона. В общем случае значение SSIM(X,Y) определяется по формуле:

где x, y – координаты пикселя, $l(x,y) = \frac{{2{{\mu }_{x}}{{\mu }_{x}} + {{c}_{1}}}}{{\mu _{x}^{2} + \mu _{y}^{2} + {{c}_{1}}}}$ – функционал яркости, $c(x,y) = \frac{{2{{\sigma }_{x}}{{\sigma }_{y}} + {{c}_{2}}}}{{\sigma _{x}^{2} + \sigma _{y}^{2} + {{c}_{2}}}}$ – функционал контраста, $s(x,y) = \frac{{{{\sigma }_{{xy}}} + {{c}_{3}}}}{{{{\sigma }_{x}}{{\sigma }_{y}} + {{c}_{3}}}}$ – функционал структуры, а μ_x, μ_y, σ_x, σ_y, и σ_xy – локальные средние значения, стандартные отклонения и перекрестная ковариация соответственно для соответствующих изображений. Критерий SSIM принимает значения от –1 до 1. Значение 1 получается в том случае, если сравниваемые изображения одинаковы [11].

Заметим, что реализация предложенной итеративной схемы возможна в двух вариантах. Первый вариант реализуется формулами (4) и (5). Во втором варианте последовательность изображений {I_q} формируется итеративным процессом I_q = ℑ(I_q-1), где ℑ(⋅) – оператор, задаваемый теми же отображениями f₁(⋅) и f₂(⋅), но действующими в другом порядке, т.е. f (⋅) = f₁(f₂(⋅)). В этом случае в качестве входного изображения используется не восстановленное в ходе решения обратной задачи распределение давления $p_{0}^{{(0)}}(r)$, а сигналы давления p_s(r, t), зарегистрированные на поверхности S зондируемого объема V ', т.е. I₀ = p_s(r, t). При этом используются формулы аналогичные формулам (4) и (5), но для операторов ℑ₁(⋅) и ℑ₂(⋅) соответственно.

В следующем разделе будут приведены результаты тестирования обеих итерационных схем.

3. ТЕСТИРОВАНИЕ АЛГОРИТМА

Для числового моделирования задачи о распространении акустических волн использовался программный пакет k-Wave – набор инструментов для среды МАТЛАБ. Его использование позволяет моделировать системы с акустическими источниками и приемниками произвольных форм и размеров. При этом числовая модель основывается на переходе в k-пространство. Пространственные градиенты в этом пространстве вычисляются с использованием схемы быстрого преобразования Фурье. При вычислении временны́х градиентов используется скорректированная k-пространственная разностная схема [12].

Задавалась числовая модель, близкая по своим характеристикам к мягким биологическим тканям: среда однородна с плотностью ρ₀ = 1020 кг/м³ и скоростью звука c₀ = 1510 м/с. Задача решалась для двумерного и трехмерного случаев.

В качестве числовых фантомов (объектов для оптоакустической реконструкции) были выбраны диск и двумерная модель сосудистого дерева для 2D случая и трехмерная числовая модель аорты с аневризмой для 3D-пространства. Физический размер образца для двумерного случая составил 4.6 × 4.6 мм для двухмерного случая и 10.3 × 10.3 × 5.3 мм для трехмерного случая. В 2D случае датчики располагались линейно на верхней поверхности прямоугольного образца, в 3D-пространстве датчики были распределены на верхней плоскости параллелепипеда. Для реконструкции заданных объектов использовался метод обратных проекций с преобразованием Фурье.

На рис. 2а показаны ОА-источник p₀(x, y) в виде кругового диска, а на рис. 2б – результат его реконструкции $p_{0}^{{(0)}}$, реализуемый k-Wave алгоритмом. Форма реконструируемого изображения отображается достаточно точно. Однако изображение отягощено артефактами в виде дуг, прикасающихся к реконструированному образу, интенсивность восстановленного образа значительно меньше, чем у исходного образца, а его кромки размыты. Особенно это заметно на рис. 2г, где заданный ОА-источник p₀(x, y) представлен в изометрии (интенсивность источника отложена по оси аппликат). Результат реконструкции $p_{0}^{{(0)}}$ на рис. 2г представлен на переднем плане.

Два средних изображения на рис. 2г $p_{{{{\Re }_{2}}}}^{{(4)}}$ и $p_{{{{\Im }_{2}}}}^{{(4)}}$ представляют 4-ю итерацию схемы I_q = ℜ₂(I_{q – 1}) и 4-ю итерацию схемы I_q = ℑ₂(I_{q – 1}), соответственно.

Более детально результаты работы алгоритма можно рассмотреть на рис. 3. Здесь представлены линейные профили реконструированных изображений $p_{0}^{{(0)}}$ и $p_{{{{\Re }_{2}}}}^{{(4)}}$ и оригинального изображения p₀(x, y).

На рис. 3а изображено нулевое приближение $p_{0}^{{(0)}}$ кругового диска p₀. Стрелками обозначены положение и ориентация вертикального и горизонтального сечений реконструированного образа. Полученные в результате линейные профили исходного объекта p₀ (черная жирная линия) и его реконструкций изображены на рис. 3б и рис. 3в. Из рисунка видно, что в результате работы алгоритма кромки восстановленного и откорректированного изображения $p_{{{{\Re }_{2}}}}^{{(4)}}$ становятся более четкими, его интенсивность практически совпадает с оригиналом, а ошибки и артефакты реконструкции, имеющиеся в нулевом приближении $p_{0}^{{(0)}}$, удаляются, что приводит к существенному улучшению отношения сигнал/помеха. Аналогичные результаты получены и для более сложных 2D и 3D объектов.

На рис. 4 представлены результаты работы алгоритма для двумерного фантома сосудистого дерева. Здесь особо следует подчеркнуть, что предлагаемый алгоритм итеративной коррекции эффективно работает как для горизонтальных, так и для вертикальных линейных структур, которые особенно плохо реконструируются при ограниченном обзоре и линейном расположении датчиков только на верхней поверхности образца (см. например, вертикальные сегменты ветвей на рис. 4б).

Аналогичный результат получается и при 3D-моделировании, когда в качестве исходного объекта, который необходимо реконструировать, была использована числовая трехмерная модель аорты с аневризмой (рис. 5).

Здесь на рис. 5а представлено оригинальное изображение аорты, на рис. 5б – результат k-Wave реконструкции (в наших обозначениях – $p_{0}^{{(0)}}$), а на рис. 5в – откорректированная версия восстановленного изображения $p_{{{{\Re }_{2}}}}^{{(4)}}$. Более подробно отличия восстановленных изображений можно видеть на рис. 5г, где представлены результаты реконструкции в центральном x–z сечении смоделированного трехмерного сосудистого дерева (обозначения линий такие же, как на рис. 3). Как и в двумерном случае, алгоритм существенно улучшает качество восстановленного объекта и гораздо точнее воспроизводит его границы, кромки и интенсивность.

Количественные оценки качества реконструкции в терминах индекса структурного подобия SSIM для описанных восстановленных объектов сведены в таблице. Очевидно, что предложенный алгоритм значительно улучшает качество реконструкции ОА-изображений, как в двумерном, так и в трехмерном случае (полужирным шрифтом выделены лучшие показатели).

О скорости сходимости предложенных итеративных схем можно судить по результатам, представленным на рис. 6. Как отмечалось ранее, в качестве критериев эффективности коррекции изображений использовались такие характеристики, как индекс структурного подобия SSIM и относительная ошибка E. Из рисунков видно, что обе предложенные итеративные схемы – I_q = = ℜ(I_{q – 1}) и I_q = ℑ(I_{q – 1}) – дают устойчивое улучшение качества изображения как в терминах SSIM, так и E с ростом номера итераций q. В обоих случаях более быструю сходимость обеспечивают алгоритмы I_q = ℜ₂(I_{q – 1}) и I_q = ℑ₂(I_{q – 1}). На рисунке представлены зависимости для различных итерационных схем при реконструкции плоского кругового диска. Аналогичные закономерности наблюдались и при реконструкции других модельных ситуаций.

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

Основной целью работы была разработка и исследование числового алгоритма, предназначенного для коррекции артефактов и искажений, возникающих в результате восстановления ОА-образов. Ставилась задача разработки алгоритма, способного компенсировать особенности конкретного метода реконструкции. Предложенная итеративная схема корректировки ОА-изображений основана на теореме Банаха о неподвижной точке.

Для исследования эффективности предложенного алгоритма была построена числовая модель оптоакустического эксперимента, имитирующая биологическую среду с встроенным в нее объектом, подлежащим реконструкции. Рассматривались случаи приемной линейной (2D случай) или плоской (3D случай) акустических антенн, расположенных на поверхности исследуемых образцов.

Для решения обратной задачи – восстановления исходных оптоакустических источников – использовался программный пакет k-Wave Matlab toolbox. Пакет позволяет моделировать среду распространения звуковых волн с помощью таких параметров, как плотность и скорость звука. Решение обратной задачи оптоакустики, полученное с помощью алгоритмов k-Wave Matlab toolbox, корректировалось итеративным алгоритмом, основанным на теореме Банаха о неподвижной точке. Были предложены четыре итеративные схемы коррекции реконструированного изображения.

Качество реконструкции определялось путем как количественной, так и визуальной оценки полученных результатов. Для количественной оценки эффективности итеративного алгоритма улучшения качества реконструкции использовались индекс структурного сходства SSIM и относительная ошибка реконструкции E. Показано, что алгоритм демонстрирует существенное улучшение качества изображения по сравнению с исходной ОА-реконструкцией.

Полученные в данной работе результаты могут быть важны с точки зрения перспектив их дальнейшего практического применения в задачах биомедицинской визуализации. Окончательные выводы об эффективности предлагаемого в данной работе подхода можно будет сделать после его тестирования на массивах реальных экспериментальных данных.

Работа выполнена при поддержке Volkswagen Foundation project “Modeling, Analysis and Approximation Theory toward Applications in Tomography and Inverse Problems”.