Сенсорные системы, 2020, T. 34, № 1, стр. 64-71

ВВЕДЕНИЕ

Вычислительная компьютерная томография, широко используемая в медицине и промышленности, позволяет исследовать внутреннюю структуру объекта без нарушения его целостности. В основе этого метода лежит численное решение обратной задачи поглощения по набору измеренных рентгеновских проекций. Одним из распространенных способов реконструкции томографических изображений является метод свертки и обратной проекции (Filtered back projection, FBP) (Kak, Slaney, 1988). Несмотря на то что этот метод менее устойчив к шумам, чем алгебраические итерационные методы, он позволяет восстановить изображение за существенно меньшее время ($O({{N}^{3}})$, где N – линейный размер изображения). Тем не менее актуальной задачей компьютерной томографии является проведение диагностики в режиме реального времени, в этом случае скорости работы FBP оказывается недостаточно.

Существует несколько подходов, позволяющих увеличить скорость реконструкции методом свертки и обратной проекции. Один из таких подходов основан на применении теоремы о центральном сечении (Fourmont, 2003; Potts, Steidl, 2000; O’Connor, Fessler, 2006). Большинство алгоритмов данной группы использует быстрое преобразование Фурье на неоднородной сетке с последующим переходом от полярных координат к декартовым в фурье-пространстве при помощи интерполяции. При этом качество реконструкции существенно зависит от выбора метода интерполяции. Подробный анализ возникающих артефактов, искажающих восстановленное изображение, приведен в работе (Potts, Steidl, 2002).

В работе (Andersson, 2005) описан метод ускорения операции обратного проецирования с использованием перехода к логполярным координатам. Показано, что в этом случае обратное проецирование может быть представлено в виде свертки с конечным ядром и выполнено при помощи быстрого преобразования Фурье.

Еще один подход к ускорению метода FBP предложен в работах (Basu, Bresler, 2000; Xiao et al., 2002). Авторы предлагают восстанавливать не все изображение целиком, а используя свойства преобразования Радона, вычислять синограммы, соответствующие четырем квадрантам изображения, и восстанавливать их по отдельности. В этом случае, согласно теореме Найквиста-Шеннона-Котельникова, для восстановления квадранта с линейным размером в 2 раза меньше исходного изображения, можно использовать в 2 раза меньше проекций без потери качества реконструкции. Разбиение изображения на квадранты можно продолжать последовательно до тех пор, пока размер независимо восстанавливаемых участков не достигнет величины в 1 пиксел. В этом случае вычислительная сложность данного алгоритма составляет $O({{N}^{2}}\log N)$. Недостаток этого метода состоит в большом количестве промежуточных интерполяций, что может приводить к накоплению ошибки (Basu, Bresler, 2001).

Все описанные методы требуют $O({{N}^{2}}\log N)$ операций умножения. В данной работе предлагается модификация алгоритма свертки и обратной проекции, позволяющая восстановить изображение за $O({{N}^{2}}\log N)$ операций сложения и $O({{N}^{2}})$ операций умножения.

МЕТОД СВЕРТКИ И ОБРАТНОЙ ПРОЕКЦИИ

В основе математического аппарата компьютерной томографии лежит интегральное преобразование Радона, которое ставит в соответствие некоторой непрерывной финитной функции $f(x,y)$, определенной на ${{\mathbb{R}}^{2}}$, множество ее линейных интегралов вдоль всех возможных прямых $L$:

В пространстве ${{\mathbb{R}}^{2}}$ любая прямая линия L может быть однозначно задана двумя параметрами. В компьютерной томографии наиболее распространено определение прямой L углом θ наклона нормали $\vec {n} = (\cos \theta ,\sin \theta )$ к этой прямой и расстоянием от начала координат r (рис. 1, а). Такая параметризация соответствует параллельной схеме томографических измерений: объект вращается вокруг своей оси в центре системы координат, а линейный детектор зафиксирован напротив источника излучения, или, что эквивалентно, объект неподвижно зафиксирован, а линейный детектор вращается одновременно с источником. Также данная система координат актуальна для среднего слоя в конусной схеме измерений, когда объект просвечивается не параллельным, а конусным пучком рентгеновского излучения.

Проекцией функции $f(x,y)$ называется множество всех точек в ее радон-образе, соответствующее определенному углу θ: ${{p}_{\theta }}(r) \equiv {{\left. {R[f](\alpha ,r)} \right|}_{{\alpha = \theta }}}.$ Важным свойством преобразования Радона является его обратимость. Это свойство позволяет восстановить исходную функцию $f(x,y)$ по набору измеренных проекций. Рассмотрим метод свертки и обратной проекции, широко применяемый для восстановления томографических изображений. Суть метода состоит в последовательном применении двух операций. На первом шаге выполняется свертка проекций с некоторой фильтрующей функцией $h(r)$:

Затем выполняется операция обратного проецирования:

где $r = x\cos \theta + y\sin \theta $. При выражение (3) соответствует точной аналитической формуле обратного преобразования Радона, и $\tilde {f}(x,y) = f(x,y)$. Данный фильтр обычно называют в честь авторов фильтром рам-лак (Ramachandran, Lakshminarayanan, 1971) или рамп-фильтром.

В реальных измерительных системах получить непрерывный набор проекций не представляется возможным, поэтому все интегральные операторы заменяются соответствующим суммированием, а непрерывная свертка преобразуется в одномерный линейный фильтр. В компьютерной томографии число проекционных углов P выбирается таким образом, чтобы оно было соизмеримо с размером изображения N. В этом случае операция свертки может быть выполнена за $O({{N}^{3}})$ операций в пространстве изображения или за $O({{N}^{2}}\log N)$ операций в фурье-пространстве с использованием алгоритма быстрого преобразования Фурье (БПФ). Операция обратного проецирования требует $O({{N}^{3}})$ операций. В следующих разделах будут описаны методы, позволяющие ускорить дискретные операции свертки и обратного проецирования.

УСКОРЕНИЕ СВЕРТКИ ПРИ ПОМОЩИ РЕКУРСИВНОГО ФИЛЬТРА

При разработке быстрых алгоритмов, выполнимых за $O({{N}^{2}}\log N)$ операций, время, требуемое для вычисления свертки, становится соизмеримым со временем, требуемым для обратного проецирования. В связи с этим становится актуальным вопрос о возможности ускорения операции свертки.

Рассмотрим идеальный рамп-фильтр, импульсный отклик которого задается выражением (4). В общем случае, он имеет сингулярность в точке $r = 0$ (Wei et al., 2005). Однако на практике спектр измеренных проекций не является бесконечным ($\left| \omega \right| < W$, где W – максимальная частота, присутствующая в спектре), что позволяет избежать сингулярности. Для дискретного сигнала $W = 1{\text{/}}2\Delta $, где Δ – шаг дискретизации. Без потери общности при выборе соответствующих координат можно полагать, что $\Delta = 1$, тогда импульсный дискретный отклик фильтра (4) приобретает вид

где $sinc(r) = \sin (r){\text{/}}r$. Последнее выражение можно упростить:

Дискретная свертка проекции ${{p}_{\theta }}$ с ядром (6) может быть записана в виде операции фильтрации с конечной импульсной характеристикой (КИХ-фильтр):

где $L = 2{{L}_{0}} + 1$ – длина ядра фильтра. Несмотря на то что значение функции $h(n)$ квадратично затухает с ростом n, уменьшение длины ядра фильтра ведет к существенным искажениям восстановленного изображения. Для достижения минимальной ошибки необходимо использовать фильтр длиной, соизмеримой с размером изображения N (рис. 2, а). Таким образом, вычисление свертки N проекций с фильтром (7) требует $O({{N}^{3}})$ операций. При этом важно заметить, что дальнейшее увеличение длины фильтра не ведет к изменению результата свертки, так как значения, лежащие на расстоянии больше N, не участвуют в вычислении.

Вопрос снижения вычислительной сложности операции свертки является актуальным во многих задачах обработки сигналов. Так, в работах (Deriche, 1987; Young, van Vliet, 1995) приведены вычислительно эффективные способы аппроксимации фильтра Гаусса, применяемого для размытия изображений, его первой и второй производных при помощи фильтров с бесконечной импульсной характеристикой (также называемых рекурсивными или БИХ-фильтрами).

Разностное уравнение, описывающее дискретный БИХ-фильтр, имеет вид

где M – порядок фильтра, ${{a}_{k}}$ и ${{b}_{k}}$ – коэффициенты, характеризующие фильтр. Преимуществом БИХ-фильтра является то, что при $M \ll N$ для обработки одной проекции требуется $O(N)$ операций. Таким образом, для свертки всей синограммы требуется $O({{N}^{2}})$ операций.

Важно заметить, что импульсный отклик фильтра (8) является однонаправленным, тогда как импульсный отклик рамп-фильтра является симметричным ($r(n) = r( - n)$). Симметричный рекурсивный фильтр может быть представлен в виде суммы (Deriche, 1993) или произведения (Young, van Vilet, 1995) казуальной и антиказуальной составляющей. Оба подхода широко применяются для аппроксимации гауссового фильтра в задачах обработки изображений. В рамках данной работы мы ограничимся применением подхода Дерише.

Представим импульсную переходную функцию фильтра (6) в виде суммы казуальной и антиказульных составляющих:

Тогда

Представим результирующий функцию ${{\tilde {p}}_{\theta }}$ в виде суммы выходов двух рекурсивных фильтров, аппроксимирующих ${{h}^{ + }}$ и ${{h}^{ - }}$:

В силу симметрии $a_{k}^{ + } = a_{k}^{ - }$ и $b_{k}^{ + } = b_{k}^{ - }$, поэтому достаточно найти коэффициенты только для казуального фильтра. Для нахождения коэффициентов будем минимизировать среднеквадратичную ошибку между откликом КИХ-фильтра (7) и БИХ-фильтра (10) на смещенную импульсную функцию вида

где $n = - {{L}_{0}}...{{L}_{0}}$. Функция такого вида позволяет более точно аппроксимировать поведения на границах по сравнению с обычной импульсной функцией.

Эта задача может быть решена любым оптимизационным алгоритмом, например, методом сопряженных градиентов Пауэлла (Powell, 1964) или симплекс-методом (Gao, Han, 2012). Найденные коэффициенты для фильтров разного порядка приведены в табл. 1. Рисунок 2, б показывает зависимость ошибки восстановления изображения от порядка фильтра. Можно заметить, что уже при использовании рекурсивных фильтров четвертого порядка ошибка становится сопоставима с ошибкой, полученной при восстановлении изображения при помощи свертки с КИХ-фильтром.

В отличие от фильтра Гаусса, рамп-фильтр не имеет шкалируемых параметров. Для обработки изображений известного размера достаточно найти коэффициенты 1 раз и не пересчитывать их для каждого изображения. На практике зачастую используют другие типы фильтров, менее чувствительные к высокочастотным шумам. Тем не менее все применяемые фильтры являются симметричными, для их представления в виде БИХ-фильтра можно использовать процедуру, аналогичную описанной выше.

УСКОРЕНИЕ ОБРАТНОГО ПРОЕЦИРОВАНИЯ ПО СХЕМЕ БРЕЙДИ $(s,t)$-ПАРАМЕТРИЗАЦИЯ

В классической реализации алгоритма FBP, операция обратного проецирования является самой затратной с вычислительной точки зрения и определяет общую асимптотическую сложность метода $O({{N}^{3}})$. В данном разделе предлагается быстрый способ вычисления дискретного обратного проецирования при помощи перехода к другой системе координат.

Введем $(s,t)$-параметризацию прямой таким образом, чтобы некоторая точка $({{x}_{0}},{{y}_{0}})$ на плоскости изображения $(x,y)$ определяла прямую на плоскости $(s,t)$. Набор проекций в пространстве параметров $(s,t)$ иногда называют линограммой (Edholm, Herman, 1987). Такая параметризация является естественной, например, для позитронной эмиссионной томографии (Hamill et al., 2003).

Пусть функция $f(x,y)$ задана в области $0 \leqslant x \leqslant N,$ $0 \leqslant y \leqslant N$. Разобьем множество всех прямых на четыре класса:

преимущественно-вертикальные прямые с уклоном вправо $L_{v}^{ + }$ $(3\pi {\text{/}}4 \leqslant \theta < \pi )$,

преимущественно-вертикальные прямые с уклоном влево $L_{v}^{ - }$ $(0 \leqslant \theta < \pi {\text{/}}4)$,

преимущественно-горизонтальные прямые с уклоном вверх $L_{h}^{ + }$ $(\pi {\text{/}}2 \leqslant \theta < 3\pi {\text{/}}4)$

и преимущественно-горизонтальные прямые с уклоном вниз $L_{h}^{ - }$ $(\pi {\text{/}}4 \leqslant \theta < \pi {\text{/}}2)$.

Тогда параметры s и t для $L_{h}^{ \pm }$ задают координаты двух точек прямой, лежащих на вертикальных границах области определения функции: $(0,s)$ и $(N,s + t)$; а для $L_{v}^{ \pm }$ – координаты двух точек, лежащих на горизонтальных границах изображения: $(s,0)$ и $(s + t,N)$ ( рис. 1, б). Параметр $t$ может принимать значения от –N до 0 для $L_{h}^{ - }$ и $L_{v}^{ - }$ и от 0 до N для $L_{h}^{ + }$ и $L_{v}^{ + }$. Между пространствами параметров $(s,t)$ и $(r,\theta )$ существует взаимно однозначная связь:

где $p = 1$ для $L_{h}^{ \pm }$ и $p = - 1$ для $L_{v}^{ \pm }$. Эта связь позволяет перейти к $(s,t)$ координатам в случае, если исходные данные заданы в виде синограммы. Выполнение интерполяции с ядром конечного размера требует $O({{N}^{2}})$ операций.

Можно заметить, что $(s,t)$-параметризация обладает морфологической симметрией: для каждой из областей определения этой параметризации точка в пространстве изображения соответствует прямой в $(s,t)$-пространстве, а точка в $(s,t)$-пространстве соответствует прямой в пространстве изображений. Такая симметрия позволяет установить связь между оператором прямого проецирования (преобразования Радона) $R[f](s,t)$ и соответствующим ему оператором обратного проецирования $B[p](x,y)$. Обозначим $P_{h}^{ + }(s,t)$, $P_{h}^{ - }(s,t)$, $P_{v}^{ + }(s,t)$ и $P_{v}^{ - }(s,t)$ линограммы, соответствующие введенным семействам прямых. Заметим, что операторы прямого проецирования для преимущественно-горизонтальных прямых ${\text{R}}_{h}^{ \pm }[f]$ могут быть получены из соответствующих операторов для преимущественно-вертикальных прямых ${\text{R}}_{v}^{ \pm }[f]$ путем предварительного транспонирования изображения (поворота против часовой стрелки и отражения относительно оси x):

Перепишем выражения для операторов прямого (1) и обратного (3) проецирования в соответствии с введенной $(s,t)$-параметризацией:

Сравнивая полученные выражения, заметим, что оператор ${\text{B}}_{v}^{ \pm }$, задающий отображение из пространства $(s,t)$ в пространство $(x,y)$, может быть выражен через оператор прямого проецирования следующим образом:

Таким образом, при использовании $(s,t)$-параметризации операция обратного проецирования эквивалентна операции прямого проецирования со сменой знака одного из параметров. Аналогичным образом, с учетом (12),оператор ${\text{B}}_{h}^{ \pm }$ записывается в виде

ОБРАТНОЕ ПРОЕЦИРОВАНИЕ ПО СХЕМЕ БРЕЙДИ

В дискретном пространстве преобразование Радона от функции вдоль заданной прямой может быть аппроксимировано суммированием значений функции в точках, принадлежащих дискретному представлению этой прямой. Такое приближение позволяет использовать ускоренные методы вычисления оператора прямого проецирования. М. Брейди заметил, что дискретные представления прямых с близкими углами наклона обладают большим количеством общих точек (Brady, 1998). В этом случае для нахождения суммы вдоль каждой из этих прямых нет необходимости считать повторяющийся участок дважды. Брейди предложил последовательно находить частичные суммы для сегментов длины ${{2}^{i}}$, $i = 1...{{\log }_{2}}(N + 1)$. Рекурсивная реализация описанного алгоритма для двумерного и трехмерного случаев приведена в работе (Ершов и др., 2017). Важной особенностью описанного алгоритма является то, что он требует $O({{N}^{2}}\log N)$ операций, причем все эти операции являются сложением, а не умножением. В работе (Прун и др., 2013) было предложено использовать схему Брейди для ускорения вычислений обратных проекций для алгебраического метода восстановления томографических изображений.

Согласно выражению (14), операция обратного проецирования может быть представлена в виде операции прямого проецирования со сменой знака одного из параметров. Это позволяет использовать схему Брейди для ускорения алгоритма свертки и обратной проекции, позволяющего восстановления исходного изображения по набору проекций, заданных в пространстве $(s,t)$. На первом этапе выполняется сверка исходных проекций $P_{h}^{ \pm }$ и $P_{v}^{ \pm }$ с заданным фильтром. Если проекции заданы в $(r,\theta )$-пространстве, выполняется переход к $(s,t)$-координатам (11). Затем для каждого набора проекций выполняется дискретное обратное проецирование по схеме Брейди в соответствии с выражениями (14) и (15):

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

На рис. 3 представлен результат работы описанного алгоритма на модельном фантоме Шеппа-Логана размером 256 × 256 пикселов. Для смоделированного изображения при помощи алгоритма Брейди была найдена полная дискретная линограмма. Изображение восстанавливалось описанным методом из набора проекций в $(s,t)$-координатах.

ЗАКЛЮЧЕНИЕ

В данной работе предложена модификация алгоритма свертки и обратной проекции, позволяющая ускорить обработку изображений. Предложена аппроксимация рамп-фильтра парой рекурсивных фильтров. Показано, что для восстановления изображения в хорошем качестве достаточно фильтров четвертого порядка. Показано, что при определенном выборе параметров, задающих прямые, оператор обратного проецирования может быть выражен через оператор прямого проецирования (преобразования Радона). Это позволяет использовать схему Брейди для быстрого выполнения операции дискретного обратного проецирования. Вычислительная сложность предложенного алгоритма составляет $O({{N}^{2}})$ операций умножения, требуемых для свертки, и $O({{N}^{2}}\log N)$ операций сложения, требуемых для обратного проецирования.

Работа выполнена при частичной финансовой поддержке РФФИ (проекты 18-29-26020, 18-29-26033).