Кристаллография, 2019, T. 64, № 1, стр. 161-166

ВВЕДЕНИЕ

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

Традиционно используемые алгоритмы реконструкции предполагают, что зондирование ведется монохроматическим излучением. То есть для зондирования используется только одна выделенная каким-либо способом длина волны рождаемого источником рентгеновского излучения спектра. Использование полихроматического излучения (используется весь рожденный источником спектр или какая-то его часть) усложняет интерпретацию результатов реконструкции, если реконструкция проводится классическими (не ориентированными под задачу полихроматического зондирования) методами [1]. Связано это с тем, что восстанавливаемой функцией в классической задаче реконструкции для метода рентгеновской томографии является функция распределения коэффициента линейного ослабления объекта, которая однозначно описывает морфологию зондируемого объекта только в случае монохроматического зондирования. То есть объект сканируется монохроматическим излучением и восстанавливается распределение коэффициента для используемой энергии. При этих условиях прямая задача томографии может быть сведена к преобразованию Радона. Тогда обратная задача или задача реконструкции заключается в восстановлении функции по ее известному преобразованию Радона.

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

МАТЕМАТИЧЕСКАЯ МОДЕЛЬ ФОРМИРОВАНИЯ ТОМОГРАФИЧЕСКОЙ ПРОЕКЦИИ ПРИ ПОЛИХРОМАТИЧЕСКОМ ЗОНДИРОВАНИИ

Задача решается в следующей постановке. Пусть двухкомпонентный объект зондируется полихроматическим излучением, спектр которого описывается функцией ${{I}_{0}}(E)$. Для данного спектра ослабляющие свойства одного компонента отличаются от свойств другого. Например, зуб с металлическим включением, допустим, коронкой, зондируется рентгеновским спектром, сформированным молибденовым анодом. Тогда материал зуба – это один компонент, а металлическая коронка – второй. Компоненты объекта пространственно разделены (в одном пикселе одновременно не может быть двух компонентов) и их поглощательные свойства затабулированы [4]. В модели используется приближение параллельной схемы измерения (рис. 1). Оно правомерно, если расстояние источник–образец достаточно велико. Рассмотрим 2D-задачу, т.е. будем восстанавливать распределение компонентов в одном из сечений объекта. Система “источник–детектор” неподвижна. Детектор позиционно-чувствительный. Объект укреплен на гониометре. Величина φ описывает текущий проекционный угол, r – определяет ячейку детектора, для которой будет записано выражение, связывающее величину регистрируемого сигнала и распределение фракций компонентов (рис. 1).

Согласно закону Бугера–Ламберта–Бера наблюдается экспоненциальное ослабление зондирующего излучения при его прохождении через объект:

Здесь $I(\varphi ,r)$ – регистрируемое излучение, пройденное через объект, повернутый на угол φ, ячейкой детектора, положение которой определено $r$. $L(r,\varphi )$ описывает траекторию регистрируемых ячейкой лучей. $\mu (l,E)$ – коэффициент линейного ослабления объекта для энергии $E$ для текущей точки $l$ на луче, заданном $L(r,\varphi )$. Данный коэффициент есть линейная комбинация вкладов ${{c}_{1}}(l)$ и ${{c}_{2}}(l)$ от двух входящих в состав объекта компонент с коэффициентами ${{\mu }_{1}}(E)$ и ${{\mu }_{2}}(E)$ соответственно: Решение системы алгебраических уравнений, для которых правой частью являются регистрируемые сигналы (1) относительно ${{c}_{1}}$ и ${{c}_{2}}$, позволит описать морфологию изучаемого объекта.

ОБ ОБРАТНОЙ ЗАДАЧЕ В ДВУХ ПОСТАНОВКАХ

Рассмотрим алгебраический подход к решению обратной задачи, т.е. будем решать алгебраическую систему уравнений [1]. Воспользуемся векторной записью. Введем матрицу весов $W$, размером ${{n}_{\varphi }}n \times {{n}^{2}}$, которая описывает геометрию эксперимента и содержит значения 0 и 1. 1 – если рассматриваемый луч пересекает пиксел, 0 – не пересекает. Здесь ${{n}_{\varphi }}$ – число проекционных углов, $n$ – число ячеек детектора соответственно, $n \times n$ – размер восстанавливаемого сечения. Векторы ${{c}_{1}}$ и ${{c}_{2}}$ размером $n \times n$ есть восстанавливаемые изображения, а компонентами вектора I размером ${{n}_{\varphi }}n$ являются измеренные значения. Рассмотрим только дискретные значения энергий используемого для зондирования спектра. Перепишем математическую модель формирования проекций, применив введенные обозначения:

где $j$ – индекс дискретного значения энергии в спектре. Решим задачу реконструкции как оптимизационную задачу большой размерности. Обозначим через $Q({{c}_{1}},{{c}_{2}})$ нормированную на величину спектра невязку между результатом измерения ${{I}^{{{\text{exper}}}}}$ и рассчитанными значениями проекций ${{I}^{k}}$ на k-ой итерации: Оптимизационную задачу можно сформулировать следующим образом, минимизируя на наборе ограничений функционал Поскольку одновременно в пикселе не могут находиться оба компонента (ожидается нулевое значение одного из вкладов), регуляризирующий член есть сумма квадратов перемноженных компонент векторов ${{c}_{1}}$ и ${{c}_{2}}$. Параметр $\beta $ регулирует вклад слагаемых. Это задача квадратичного программирования, но из-за своей высокой размерности стандартными программными реализациями решателей QP она не может быть решена. Поэтому был выбран подход, в котором ограничения учитываются с помощью логарифмических барьеров: Вдали от границы логарифмические штрафы не влияют на оптимизируемый функционал, а при подходе к границе очень быстро растут, внося соответствующие поправки в шаг итерации [5]. Параметр $t$ регулирует степень близости решения к границам. В ходе итерационной минимизации постепенно увеличивается значение $t$, что позволяет уменьшить влияние логарифмических членов и приблизиться к их границам, а значит, и к истинному значению минимума. Метод градиентного спуска использован для решения оптимизационной задачи относительно ${{c}_{1}}$ и ${{c}_{2}}$. Метод вычисления градиента функции $Q({{c}_{1}},{{c}_{2}})$ представлен в [6].

В серийных томографах для решения задачи реконструкции применяется метод свертки и обратной проекции (Filtered Back Projection (FBP)) [1]. Метод выбирается производителями за его быстродействие. Результат работы метода есть пространственное распределение некоторого усредненного коэффициента линейного ослабления зондируемого объекта. Рассмотрим приближения, без которых метод не применим, и обсудим последствия их влияния на результат реконструкции. Обратимся к выражению (1) и перепишем его, введя оператор преобразования Радона:

где ${{\mu }_{E}}(\varphi ,r)$ – точки, лежащие на траектории луча $(\varphi ,r)$. Дальнейшая работа ведется с результатом логарифмирования $P(\varphi ,r)$. Согласно [7], это Авторы [8] предлагают проводить нормировку И далее все авторы сходятся к векторному решению называя ${{R}^{{ - 1}}}$ оператором свертки и обратной проекции. Отметим, восстанавливаемое изображение – это распределение “усредненного” по спектру коэффициента ослабления материалом зондируемого объекта. Для всех химических элементов зависимость коэффициента ослабления от энергии нелинейна. Положение и величина скачка поглощения на энергетической оси являются уникальным кодом для элемента. В связи с этим интерпретация полученных результатов остается отдельной трудоемкой задачей, часто требующей привлечения знаний из области, откуда поступил исследуемый объект. То есть интерпретацию результата трудно автоматизировать. Последнее является желательным во многих областях применения метода компьютерной томографии. Есть и другая проблема – возникновение артефакта типа “эффект чаши” [9] на восстановленном изображении. Его появление обусловлено тем, что излучение мягкого и жесткого диапазонов (в спектре есть и то и другое) ослабляется по-разному, т.е. происходит так называемое ужесточение пучка (Beam Hardening (BH)). Компенсация эффекта ужесточения не может быть реализована в ядре метода FBP. Поэтому разными группами проводятся исследования, целью которых является создание эффективных методов BH-коррекции на этапе предобработки данных, посылаемых на реконструкцию. В частности, недавно был предложен алгоритм автоматического поиска параметров BH-коррекции, базирующийся на использовании Радоновского инварианта и не требующий эталонных измерений [10].

В следующем разделе описан предлагаемый подход к интерпретации результатов компьютерной томографии при полихроматическом зондировании двухкомпонентных объектов, позволяющий в будущем автоматизировать интерпретацию результата. Карты распределения компонентов анализируются совместно с результатом, полученным методом FDK [11].

ОБСУЖДЕНИЕ РЕЗУЛЬТАТОВ

Приведенные результаты модельного эксперимента получены для двухкомпонентного фантома. Пространственные распределения компонентов ${{c}_{1}}$ и ${{c}_{2}}$ приведены на рис. 2. Зависимости коэффициентов ослабления от энергии для каждого из компонентов и спектр, использованные при расчете синограммы, представлены на рис. 3. На рис. 4 приведена синограмма, рассчитанная с использованием выражения (1). Реконструкция проводилась двумя методами: FBP (рис. 5а) и алгебраическим (рис. 5б, 5в). При решении оптимизационной задачи использовались следующие значения параметров: β = 0.05, t = 1. Процедура выбора значений состоит в нахождении таких величин, чтобы выполнялись сразу два условия: значение градиентов по всем трем слагаемым выражения (6) соизмеримы, т.е. ни одно из слагаемых в суммарном шаге не доминирует; суммарный шаг не выводит за границы целевого множества. В экспериментах параметры подбирались экспертом вручную, в дальнейшем их подбор планируется осуществлять автоматически. Если за 20 итераций значение функционала (6) изменялось меньше чем на 0.001, то величина параметра $t$ удваивалась. Начальное значение функционала – 2011.815, значение на момент остановки – 19.232. Использованное для реконструкции число итераций равно 1800. Результаты реконструкции ${{c}_{1}}$ и ${{c}_{2}}$ алгебраическим методом представлены на рис. 5б и 5в соответственно. Отметим, что три изображения, приведенные на рис. 5, дополняют друг друга. Результат реконструкции методом FDK предъявляет весь объект целиком, не давая никаких указаний к интерпретации. В то время как рис. 5б точно указывает местоположение первого компонента, а распределение ${{c}_{2}}$ (рис. 5в) показывает в дополнение к местоположению второго место отсутствия первого. Однако вывод о том, что гало вокруг места отсутствия первого есть артефакт, можно сформулировать лишь после сравнения результатов на рис. 5в и 5а. В дальнейшем планируется продолжить исследования с целью уменьшения величины артефактов и создания метода совместной интерпретации получаемых разными методами результатов.

ЗАКЛЮЧЕНИЕ

Описан новый подход к решению задачи томографической реконструкции при сканировании объектов полихроматическим излучением. Поскольку результат реконструкции при использовании традиционных алгоритмов в таких условиях представляет собой распределение некоторого усредненного коэффициента ослабления рентгеновского излучения, полученное распределение не всегда удается легко связать с морфологией объекта. Предложено решать задачу реконструкции в постановке, когда восстанавливаются распределения компонентов объекта. Задача реконструкции решается алгебраическим методом. Функционал, содержащий невязку и регуляризирующий член, минимизируется на наборе ограничений. Обсуждаются результаты модельного эксперимента с двухкомпонентным объектом, в котором продемонстрирована потенциальная возможность разделения компонент. Дополнение ансамбля интерпретируемых изображений результатом реконструкции методом FBP позволяет выделить на восстановленных распределениях компонентов сомнительные области. Одним из направлений дальнейшей работы будет создание алгоритма совместного анализа результатов реконструкции двумя методами. Во время работы над текстом статьи была опубликована работа [11], в которой приведен обзор последних аппаратных и программных подходов к решению задачи, возникающей в томографии в условиях полихроматического зондирования. Это подчеркивает необходимость создания отдельных инструментов при работе в таких условиях и актуальность представленных в данной статье работ.

Авторы выражают благодарность Н.П. Власовой за техническое содействие в подготовке материалов (ФИЦ ИУ ИСА РАН).

Исследования, посвященные способу модификации алгебраического метода томографии для работы с мультикомпонентными данными, и численная реализация использованных алгоритмов проводились при финансовой поддержке Российского научного фонда (проект № 14-50-00150). Реализация модельных расчетов, включая оптимизацию используемых распределений, выполнялась при поддержке Федерального агентства научных организаций (соглашение № 007-ГЗ/Ч3363/26).