Хермен Г. Восстановление изображений по проекциям: Основы реконструктивной томографии

Подождите немного. Документ загружается.

178

ГЛАВА 9

т.е. при этом гарантируется, что величина [Лр}(г, ф) произвольно мала

при достаточно большом г.

На практике величину Sp можно вычислить лишь на конечном мно-

жестве точек. С учетом того, что следующим шагом будет нахождение

двумерного фурье-образа от &р

имело бы смысл вычислить значения @р

в эквидистантных точках на прямоугольной сетке (вспомним наше рассмот-

рение БПФ в двумерном случае), совмещенной с центрами ячеек и образую-

щей,

таким образом, кадр изображения. Оставшаяся часть сетки будет вы-

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

статочно большой, чтобы с полной гарантией можно было бы пренебречь

значениями &р вне этой области.

На второй стадии используем алгоритм БПФ для вычисления фурье-

образа ^2^/7 от &р. Как уже отмечалось выше, число точек, в которых

применяется алгоритм БПФ $2&р

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

ли значения ьвр. Пространственное разделение точек (вдоль координатных

осей),

в которых вычисляются значения &

@р

обратно пропорционально

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

жения и свободного поля вокруг него, которые были введены выше). Вы-

числение величин

б£р

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

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

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

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

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

сверточного алгоритма и фурье-алгоритма, может получить их самостоя-

тельно.

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

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

стого умножения [^2^р](/?, Ф) на \R\ необходимо умножить его на

II?"I х F

x/d

(\R\), где F

l/d

— функция «окна» с шириной полосы про-

странственных частот, обратно пропорциональной шагу между отсчетами.

Отметим важное свойство: F(0, Ф) = 0. Имея в виду соотношения (9.1),

(9.3) и (9.21), можно показать, что

И / *(г, 0)

<fr

<ty = 0, (9.23)

т.е. общая суммарная плотность изображения/* = 0, что весьма нежела-

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

общая плотность в оценке остается одной и той же. Сделаем два замеча-

ния по этому поводу.

Во-первых, соотношение (9.21) до некоторой степени двусмысленно.

Функция/* определена так, что она может всегда иметь (и обычно имеет)

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

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

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

ДРУГИЕ АЛГОРИТМЫ РЕКОНСТРУКЦИИ

179

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

реконструировать

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

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

ма свертываемые проекции р

тА

* &~

Ф имели нулевую суммарную плот-

ность [формула

(8.31)].

В приложении к фурье-алгоритму, описываемому

выражением (9.14), величина [$&yp](0

т) всегда умножается на нуль, и

значение суммарной плотности р

тА

теряет смысл. Тем не менее оба упомя-

нутых алгоритма дают реконструированные изображения, в которых сум-

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

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

фильтрации обратного проецирования, эта проблема зачастую оказывается

нерешаемой так же эффективно, как и в других алгоритмах реконструкции.

Поскольку на четвертой стадии используется алгоритм БПФ, вычисления

завершаются расчетом значений /* на прямоугольной сетке, совпадающей

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

го проецирования на первой стадии. Суммы значений /* на первой сетке

равны нулю. Если только она ненамного больше сетки, связанной с облас-

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

изображения будет существенно неточной. К счастью, этот эффект легко

можно исправить с помощью метода дополнительной нормировки (разд.

7.2). Заметим, что в этом случае мультипликативная нормировка непригод-

на.

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

ем сверточного алгоритма с обратной последовательностью операций.

ПРИМЕЧАНИЯ И ССЫЛКИ

Соотношения между преобразованиями Радона и Фурье были рассмотрены в

ряде математических статей (например, в [111] и библиографии к ней). Первое упо-

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

изображения дано в [13].

Эквивалентность соотношений (8.23) и (9.15) вытекает из теоремы свертки для

дискретного фурье-образа (например, [151).

Ряд исследователей, например Шепп и Логан

[142],

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

разований, в том числе и сверточный алгоритм, методом «фурье-реконструкции».

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

ния преобразования Фурье (хотя и то и другое в принципе возможно), данное назва-

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

сверточного алгоритма и фурье-алгоритма реконструкции существенны, и это требу-

ет для них различных названий. Хорошее представление обо всем том, что известно

о фурье-алгоритме, можно получить из знакомства с работой [117] и библиографии

к ней, особенно из работы [36]. В докладе [113] в некоторых экспериментальных ре-

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

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

[135].

В рабо-

те [155] обсуждалась возможность использования фурье-алгоритма для небольшого

числа ракурсов.

180

ГЛАВА 9

Алгоритм БПФ имеет чрезвычайно большое значение для обоих рассмотренных

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

священы целые книги, например [19].

Термин «rho-filtered layergram» был предложен в

[147].

Дальнейшее развитие ал-

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

[140],

в которой также содержа-

лось описание экспериментов с использованием данной методики и различного рода

функций «окна» и интерполяционных функций. Читатель может изучить этот воп-

рос по соответствующей литературе. Алгоритм получил обобщение на случай ре-

гистрации данных в расходящемся пучке в [36], а также в [23].

Менее существенные подробности рассмотренных выше алгоритмов реконструк-

ции читатель может найти в работе [79], из которой взяты рис. 9.4 и 9.5, а также в

работе [91].

Алгоритмы реконструкции сверточного типа

для веерных пучков

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

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

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

точки, и на каждом из них вычисляется лучевая сумма. Стандартные про-

екции являются данными того же типа.

Существует два основных подхода к разработке сверточных алгоритмов

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

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

ния по Радону, пригодного для использования для случая веерного пучка.

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

( f, 0)-пространстве для вычисления лучевых сумм на множестве парал-

лельных лучей по измеренным лучевым суммам на веерных лучах (этот

процесс интерполяции называется повторным разбиением), а затем приме-

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

В данном разделе обсуждаются оба эти подхода. Здесь мы также вер-

немся к вопросу о выборе функции «окна», на этом раз применительно к

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

10.1.

СВЕРТОЧНЫЙ АЛГОРИТМ ДЛЯ ВЕЕРНОГО ПУЧКА

Геометрия схемы регистрации исходных данных, которая будет рас-

смотрена, изображена на рис. 10.1. Источник рентгеновских лучей S посто-

янно находится на окружности радиусом & с центром в начале координат,

а детекторы располагают на полоске, имеющей форму дуги окружности с

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

множества расходящихся лучей с индексами а, /3, где величина /3 определяет

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

исходящей из источника луч. Подобное обозначение линий отличается от

ранее использованного обозначения f, в (в частности, на рис. 2.8). Безус-

ловно, каждому лучу с параметрами а, /3 соответствует луч с параметрами

(

0, причем 1 и 0 зависят от а и /3. Действительно, как легко видеть из

рис. 10.1, лучу с координатами а, /3 соответствует точка с координатами

Dsino, /3 + а в ( (, ^пространстве.

Как принято в этой книге, функция / двух полярных координат г и ф

обозначает искомую реконструируемую функцию. Напомним, что функция

182

ГЛАВА 10

Рис. 10.1. Геометрия системы регистрации исходных данных в веерном пучке.

Каждый луч пучка характеризуется двумя параметрами /3 и о. Пусть точка О совмещена с нача-

лом координат, а точка 5 — с источником реитгеиовских лучей, который всегда находится иа

окружности радиусом D с центром в точке О. При этом угол /3 + т/2 равен углу между отрез-

ком OS и осью В,

а — угол между рассматриваемым лучом и отрезком SO. Луч SO является

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

лить и параметры луча в и f. Пусть Р — точка пересечения рассматриваемого луча с перпендику-

ляром, опушенным из начала координат О иа луч. При этом в — угол между отрезком ОР и

осью В, a f — длина перпендикуляра ОР (Воспроизведено из работы [70] с согласия автора.)

Гу

ф) = о при г ^ Е. Предположим, что D > Е и, кроме того,

Ь = arcsin

(E/D),

где

0 < S = arcsin (E/D) < л/2. (10.1)

Для обоэначения интеграла от

функции

/ вдоль луча с а, /3 введем величину

Я(о,

р) =* [ #n(D sin

<т,

р + о). (Ю.2)

что следует из соотношения (10.1), если g(o, 0) = 0 при \о\ > 6.

Предположим, что у нас имеются все возможные проекции g(o, 0) в ин-

тервале значений О^0<27гиЫ^6. Можно показать, что формула об-

ращения Радона приводит к следующему соотношению: если обе произво-

дные функции g(a, 0) по переменным о и 0 существуют и обозначаются че-

рез g

}

(o, 0) и g

(o, 0) соответственно, то справедливо равенство

f(r.<*>)

А Г Г

~^—

<*'A

№ <

л Jo J-oo ° ~

где

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

183

G(c,

а) =

sin(<x

—

о) \W W J

sin

Г COS(/J — ф) Л Я

с = arete , — < с < ,

при

<5,

при

(10.4)

(Ю.5)

W = {[г cos(/? - ф)]

+ [D + г sin(/? - ф)У}

'\ W>0. (10.6)

Отметим, что величины a' \iW зависят от

/3,

гиф, но не от а. Геомет-

рический смысл о i\W состоит в следующем: если источник располагается

под углом /3, то лучи, которые идут по направлению г, ф, имеют парамет-

ры а', /3, а расстояние между источником и направлением, определяемым

г, ф, равно ИЛ

Отметим также, что формула (10.3), называется формулой

обращения

Радона для веерного пучка, предстает в форме сингулярного интеграла,

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

вание Гильберта функции G(a

/3, а') по переменной а. Следовательно, при-

веденный в разд. 8.1 анализ применим и в данном случае.

Предположим, что

G(a

/3, а'), рассматриваемая как функция перемен-

ной а, достаточно «хорошая» в точке а'. Тогда при произвольном выборе

множества [F

\A > 0}, удовлетворяющего условиям формулы (8.10), по-

лучаем семейство регуляризирующих функций \р

\А > 0), таких, что вы-

ражение (10.3) можно записать в форме

/(г,

ф) = lim

1 л2я ла

i?Jo

J-

(о'

- o)G(o

о') do dp

(10.7)

[формулы (8.9) и

(8.10)].

Преимущество перехода от записи (10.3) к (10.7)

состоит в том, что при любом значении А внутренний интеграл в выраже-

нии (10.7) легко вычисляется методом интегрирования по частям. Обознаь

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

функции

при

фиксированном А

через

т.е.

/>'

Ф)

= 4^2 J "Г

Рл(о'

*)G(ff,

Р.

«О

dp.

(10.8)

Выбирая некоторую функцию р

в качестве регуляризирующей, под-

ставляя значение G(o, /3, а') в соотношение (10.8) и интегрируя по частям

(при этом функция р

считается дифференцируемой [формула

(8.14)],

полу-

чим при и # 0

184

ГЛАВА 10

+ q

\о' - <r)cos о']д{о. ff) da

с!Ц

(10.9)

где

\u)= -ир

(и) sin

и, (10.10)

\u) = (р

(и) + M/^(M))/sin м, (10.11)

Выражения (10.10) и (10.11) дают предельные значения величин при и = 0.

(Часть громоздких выкладок при выводе была опущена.) Уравнение (10.9)

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

веерного пучка.

Таким образом, мы получили далеко не простой алгоритм реконструк-

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

При этом необходимо, во-первых, выбрать семейство функций «окна» F

и, во-вторых, установить величину А. Как только такой выбор будет сде-

лан, можно найти сворачивающие функции <7

(1)

и <7

(2)

согласно формулам

(10.10),

(10.11) и получить аппроксимацию для функции/(г, 6) по формуле

(10.9).

До сих пор для простоты анализа предполагалось, что проекции су-

ществуют при всех возможных значениях параметров а и /3, тогда как на

практике мы имеем лишь конечное число положений источника, а для каж-

дого из них — конечное число регистрируемых детекторами точек. Как по-

казано на рис. 5.2, проекции регистрируются в М-равноотстоящих по /8

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

в (2/V + 1)-равноотстоящих по углу точках с шагом X. Таким образом, ве-

личина g известна в точках с угловыми координатами лХ, тА в интервалах

значений

—

N^n^N,0^m^M

—

1, причем МА =

2ж.

Стандартные

проекции получаются именно этим способом (разд. 5.2). Даже несмотря на

то,

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

g(nX

тА), полученные в ре-

зультате эксперимента, мы используем те же обозначения g(n\, тА) для

этих оценок в оставшейся части главы.

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

ние по формуле (10.9) проводится в два этапа.

Вначале с помощью представления интеграла в виде интегральной сум-

мы по переменной а', умноженной на X, вычисляется внутренний интеграл

в выражении (10.9), в результате чего получаем соотношение вида

QcWKmtS) = Я £ cos(nA)q(n^ mA)q

((n - и)А)

п= -N

+ Acos(w7) £ д(пК тА)ч

\{п' - п)А). (10.12)

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

185

которое соответствует формуле (8.23). Заметим,

что

первая сумма пред-

ставляет собой дискретную свертку функции <?

(1)

исходной проекции, взве-

шенную

функцией косинуса,

вторая

—

просто дискретную свертку функ-

ции <7

(2)

с той же

проекцией.

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

/*(',<«

ъ I

0с(*'."'Д)- (10.13)

где

а' и W

определяются

из

соотношений (10.5)

(10.6)

при

замене

в них

тА

на

&. Вычисления

по

формуле (10.13) включают

себя интерполяцию

вычисленных значений

(o'

тА) по

значениям

(n'\

тА).

Назначение

операции интерполяции

уже

детально обсуждалось

разд.

8.5.

Заметим,

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

как

операцию «взве-

шенного обратного проецирования».

данной точке

координатами

г, ф и

при данном положении источника

тА

луч

параметрами

а', тА

соответст-

вует именно тому лучу, который проходит через Источник

тА в

точку

г, ф.

Вклад сворачиваемой лучевой суммы

(o',

тА) в

величину/*

точке

с ко-

ординатами

г, ф,

через которую проходит

луч,

обратно пропорционален

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

с г, ф и

источником

с тА.

Реализация данного алгоритма

на

ЭВМ такая

же, что и

при использова-

нии сверточного алгоритма

для

параллельного пучка,

что

следует

из

соот-

ношений (8.24)

(8.25), поэтому этот вопрос повторно рассматриваться

не

будет,

внимание читателя будет остановлено

на

выборе функции «окна».

10.2.

ВЫБОР СВОРАЧИВАЮЩЕЙ ФУНКЦИИ

В разд.

8.6 мы

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

чиваемой

интерполяционной функций

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

лельного пучка. Было показано,

что

линейная интерполяция дает лучшие

результаты,

чем

метод интерполяции

по

ближайшим значениям.

случае

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

ции

без

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

Задача настоящего раздела состоит

ответе

на

вопрос

том, каким

об-

разом выбор функции «окна» сказывается

на

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

изображения.

При

ответе

на

этот вопрос возможны

два

подхода,

имен-

но:

а) точно сформулировать критерии качества реконструированного изо-

бражения

найти функции «окна», оптимальные

точки зрения данных

критериев;

б) взять представительный набор функций «окна»

исследовать

их в со-

ответствии

некими фундаментальными критериями, такими,

как

эффек-

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

устойчивость

его к шу-

му,

также

соответствии

с их

«поведением»

при

работе

тест-

объектами, характерными

для той или

иной области применения метода.

186

ГЛАВА 10

Первый подход выглядит менее привлекательным с математической

точки зрения и не позволяет эффективно производить выбор функции «ок-

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

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

дующее требование: «Мне необходимо найти функцию «окна», которая да-

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

головы по ее рентгеновским проекциям?» Легче сформулировать общие

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

«оптимальная» функция «окна» может сильно зависеть от величин, харак-

теризующих эти параметры (примеры будут приведены ниже). После того

как критерий, согласно которому определяется оптимальность функции

«окна», зафиксирован, нахождение последней может оказаться весьма

сложным, поэтому при ее формулировании обычно делают упрощающие

предположения (такие, как возможность применения проекций по всем на-

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

нала).

Все это заставляет сомневаться в том, что выбранная функция «ок-

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

ходных данных. Однако даже с использованием упрощающих предположе-

ний часто трудно найти в замкнутой форме решение, определяюшее вид

функции «окна».

По причинам, перечисленным выше, здесь будет использован второй

подход к исследованию влияния вида функции «окна» на качество рекон-

струированного изображения. При выборе функций «окна» необходимо

принять во внимание следующие соображения.

В идеальном случае необходимо выбрать такую функцию «окна», кото-

рая равнялась бы единице в пределах своей полосы пропускания, а ширина

самой полосы была бы по возможности большой. Необходимость в ис-

пользовании функции «окна» иного вида состоит в уменьшении влияния по-

грешностей, обусловленных использованием конечного числа ракурсов, эф-

фектами дискретизации проекций и различного рода шумами. Как правило,

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

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

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

шума обычно выражены более существенно на высоких пространственных

частотах, чем на низких. «Хорошая» функция «окна» позволяет сгладить

указанные эффекты.

В соответствии с условиями, накладываемыми на вид функции F

, а

также со сделанными выводами о необходимости подавления ошибок без

излишне больших потерь информации функции «окна» должны представ-

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

(1 - 2U/A) ^ ^(1/) ^ 1 на интервале [О, А/2], Для исследования области

возможных значений функции рассмотрим функцию обобщенного «окна»

Хэмминга при трех значениях параметра а, а именно: 1,0, 0,8 и 0,54 (табл.

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

ют равной 1/а по соображениям, которые понятны из анализа разд. 8.6.

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

187

Для каждой из упомянутых выше функций «окна» были проведены по

три эксперимента: один — для исследования способности метода рекон-

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

чивости метода и, наконец, третий — для оценки эффективности метода

реконструкции поперечных сечений головы человека в рентгеновских лучах.

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

схема регистрации исходных данных, описанная в разд. 5.2, а реконструи-

рованная функция вычислялась в центрах 115-элементной сетки (разд. 4.1).

Шаг между отсчетами составлял 0,1504 см. Обсуждение этих трех экспери-

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

10.3.

ФУНКЦИЯ ИМПУЛЬСНОГО ОТКЛИКА

Поскольку процедура получения проекций и реконструкции изображений

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

ной операцией (разд. 6.3), то оптимальным алгоритмом реконструкции

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

исходного изображения, взятого в виде одиночной точки. К сожалению, из

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

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

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

при трансляции. Несмотря на то что форма распределения в реконструиро-

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

рассмотрим в качестве примера реконструкцию точечного объекта, поме-

щенного в начало координат. При этом лучевая сумма для центрального

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

Величину Р

(г, ф) мы используем для обозначения функции импульсного

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

Эта функция вычисляется при помощи подстановки исходных дискретных

проекций в формулы (10.12) и (10.13).

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

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

бщенного «окна» Хэмминга с а = 0,54. Осцилляции функции отклика вбли-

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

собе отображения. (Более наглядное представление о ней дают графики

рис. 10.3.)

В идеальном случае функция Р*(г, ф) должна равняться нулю повсюду

при г Ф 0, однако на практике дело обстоит иначе. Нежелательный эффект

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

жении одиночной точки, имеющей достаточно большую интенсивность

(рис. 10.2), могут маскировать реконструированное изображение другой

точки с меньшей интенсивностью. Это обстоятельство побуждает нас рас-

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