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

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

-0,60 -0,40 -0,20 0,00

Частота

0,20 0.40 0,60

-0,60 -0,40 -0,70 0,00

Частота

0,20 0,40 0,60

Рис. 8.1 (продолжение).

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

139

lim [ф *

](V)

= [Jf

£|(ю).

(8.9)

/4-оо

Численную оценку гильберт-образа

[<&Ф]{и)

теперь можно произвести

путем вычисления свертки [ф *

](v)

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

точно большим, с тем чтобы левая и правая часть соотношения (8.9) были

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

об-

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

по величине, дифференцируема и т.д., то значение свертки [ф *

](v)

можно найти достаточно просто. В частности, если ф — производная неко-

торой другой функции, то значение [ф *

](v)

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

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

мула

(8.6)].

Предыдущий анализ основан на том, что семейство регуляризуюших

функций существует. Как можно показать ниже и строго доказать в

Рис. 8.2. Графики функций р

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

на рис. 8.1.

Ось и отмечена на графиках как ось расстояний, ось р(ы) — как ось величин. Отметим, что функ-

ция р антисимметрична, т.е. р(-и) = -р(и) (формула

(8.10)].

3 .

*0 ,00

2,00 3,00 4,00

Расстояние

'0,00 1,00

г ,oo э,оо 4,оо

Расстояние

5,00

6,00

Рис. 8.2 (прололжение)

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

141

разд.

16.1, существует большое разнообразие упомянутых семейств. Это

утверждение выражает следующая

теорема

регуляризации:

Пусть для каждого вещественного значения А существует вещественная

интегрируемая функция F

, удовлетворяющая при и ^ 0 условиям:

а) 0 ^ F

(U) < 1, F

(U) = 0, если U > А/2;

б) F

(U) — монотонно невозрастающая функция U;

в) lim F,(U) = 1.

А - <» "

Пусть

/>» = -2 Г F

(V)sm(2nUu)dU, (8.10)

Тогда множество

IP^'^

> 0) является семейством регуляризующих

функций.

По соображениям, которые станут понятными ниже, функцию

обыч-

но называют

функцией

«окна» с

шириной

полосы

пропускания

А.

Нетрудно получить семейство функций «окна», удовлетворяющих усло-

виям «а» — «в». В табл. 8.1 приведены некоторые из них и даны их наибо-

лее часто употребительные названия. Во всех случаях, согласно условию

«а», определения справедливы лишь при 0 ^ U < А/2, поскольку при U ^

^ А/2 функция F

(U) = 0. Заметим, что ограниченная по полосе частот

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

«окна» Хэмминга при значении параметра а = 1. Форма отдельных функ-

ций «окон» представлена на рис. 8.1, а соответствующие им функции р

—

на рис. 8.2.

8.2. ВЫВОД СВЕРТОЧНОГО АЛГОРИТМА

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

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

и преобразование Гильберта [формула

(6.15)],

можно заменить простой

сверткой исходных данных с фиксированной функцией ядра.

Пусть, как и в разд. 6.2, функцию двух переменных р необходимо под-

вергнуть обратному преобразованию Радона. Пока считаем переменную в

фиксированной и введем обозначение

Рв(П =

р(^0).

(8.11)

Далее для любого V имеем

1**&ур]('\0)=1#р-

](/-)

(8.12)

где

р'

— производная функции р

Согласно формуле (8.9), выражение в правой части (8.12) можно заме-

нить выражением вида

Ip'e

* PAW) = f Pi('W - ') if-

(8.13)

—

142

ГЛАВА 8

Видно,

что

р'

(

() = 0 при I t I ^ Е.

Если функция

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

ношением (8.10),

функция

удовлетворяет условиям разд.

8.1 («а» —

«в»),

то

производная р'

функции

существует

равна

р» = -4л f UF

(U)cos(2nUu)dU.

(8.14)

(Здесь был использован стандартный прием, который

анализе обычно

из-

вестен

как

метод дифференцирования

под

знаком интеграла.) Используя

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

правой части

(8.13)

по

частям

получим

[р'е

* Рл)(П = Г РвООрЖ' - П М. (8.15)

J -00

•

Определим функцию р

при

помощи соотношения вида

PAV

\в) = [р

* р'ЛО-

(8.16)

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

что

гильберт-образ

[Jtf

P](

аппроксимируется сверткой

( £ ', 0).

Последнее позволя-

ет определить оценку/*, даваемую сверточным алгоритмом

/• = -(\/2n)J$p

(8.17)

[формула (6.15)J.

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

нам

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

точного алгоритма

компактной форме.

Введем новый оператор

который назовем оператором свертки

по

первой переменной, связывающим функцию

двух переменных

функцию

q одной переменной

некоторой новой функцией

уЩ

двух переменных

согласно определению

[p*

qX/

0) =

*qW\

(8.18)

где функция

определена соотношением (8.11).

Математической абстракцией

при

реконструкции изображения

исполь-

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

(со

сворачивающей функцией

является

оценка

/*,

определяемая соотношением

J6(p*

q). (8.19)

При условии

что

ф) = -(1 2n)p

(U),

(8.20)

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

143

функции, определяемые соотношениями (8.17) и (8.19), тождественны друг

другу.

Подстановка последнего выражения в формулу (8.14) дает

[Л/2

q{u) = 2 UF

(U)cos(27iUu) dU.

(8.21)

На рис. 8.3 приведены графики функции q, вид которых определяется вы-

бором различных функции F

Резюмируя сказанное, отметим, что реконструкция изображения с ис-

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

вание Радона, выполняемое в два этапа:

а) сверткой по первой переменной;

б) обратным проецированием.

Функцию ядра на первом этапе обычно выбирают, исходя из выражения

(8.21),

причем F

— одна из функций, принадлежащая семейству функций

'0,00

3,00 4,00

Расстояние

в,оо

Рис.

8.3. Графики функций q

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

рис.

8.1.

Функции q имеют непрерывную производную. Ось и отмечена на графиках как ось расстояний,

ось д(и) — как ось величин. Заметим, что функция q симметрична [формула

(8.21)].

Ломаная ли-

ния на графике соединяет те значения функции q

при которых она берется в дискретной свертке

(8.23) в предположении, что d = \.

'tfjgo

Г^оо

?Joo

5^00 4^00 Б^ОО e,oo

Расстояние

«0

,00

2,00 3,00 4,00 5,00 в,1

Расстояние

Рис. 8.3 (продолжение).

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

145

«окна» и удовлетворяющая условиям разд. 8.1 («а» — «в»). Вопрос о выбо-

ре функции «окна» обсуждается в разд. 8.6.

8.3.

ПРИМЕНЕНИЕ СВЕРТОЧНОГО АЛГОРИТМА

В разд. 7.2 обсуждался вопрос о том, как применять оператор обратно-

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

функция р известна в точках с координатами (nd, тА), где —N ^ п ^ N,

О ^ т ^ М

—

1 и МА = 7г. Для использования развитого выше метода не-

обходимо рассчитать значения свертки р *

q в указанных точках по значе-

ниям функции р в тех же точках.

Комбинируя выражения (8.18) и (8.1), получим

\_р

*у q]{nd. тА) = р(Л mA)q(nd - /) d/. (8.22)

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

этому можно предположить, что р( 2, в) = 0 при I е

^ £; кроме того,

напомним использованное нами в разд. 7.2 предположение о том, что

Nd ^ £. Далее, аппроксимируя интеграл в правой части (8.22) интеграль-

ной суммой Римана, получим

(rul

тА) = d X Р("Ж mA)cj((n - n)d). (8.23)

п = -N

Сумму в (8.23) называют дискретной сверткой.

Теперь соединим вместе рассмотренные выше операции в общую про-

цедуру получения оценки /*

(г,

ф) по проекциям, измеренным вдоль парал-

лельных, эквидистантных лучей с одинаковыми углами наклона. Точнее,

мы предполагаем, что имеются оценки p(nd, тА) величин [&f](nd, тА)

для целых значений пит, лежащих в интервалах -N ^ п ^ N и 0 ^ т ^

^ М- 1, т.е.:

а) Для каждого значения т, лежащего в интервале 0 ^ т ^ М

—

1, по

формуле (8.23) вычисляем величины p

(n'd, тА) для

—

N ^ п' ^ N. Число

отсчетных точек для р

равно / = (2N + \)М. Связанный с р

вектор у

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

б) Оценку /*(г, ф) получают из соотношения (7.2), в котором функция р

заменена на р путем интерполяции значений p

[r Qos(mA - ф), тА] по ис-

ходным значениям p

(n'd, тА).

Здесь имеется ряд вычислительных особенностей, на которые хотелось

бы обратить внимание читателя.

Поскольку значения п и п' лежат в интервалах -N ^ п ^ N и

—

N ^ п' ^ N соответственно, то необходимо, чтобы значения функций q

были известны по меньшей мере в 4/V + 1 точках. Для системы регистра-

ции проекционных данных с фиксированной конфигурацией элементов то-

146

ГЛАВА в

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

зовать в сочетании с большим числом измеренных проекций. [В самом де-

ле,

обычно имеет* место свойство четности функции <?, т.е. q(

—

v) = q(v);

формула (8.21).] При этом функцию q необходимо вычислять лишь, в

IN + 1 точках. (См. также рис. 8.3 в некоторых точках для функции «ок-

на», изображенных на рис. 8.1.)

Для любого фиксированного т значения p

(n'd, тА) (—N ^ л' ^ N)

получают из величин

p(nd

тА) (—N ^ л ^ 7V). Следовательно, проекция,

полученная под одним ракурсом (т.е. состоящая из множества параллель-

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

курсом. Вычисления по формуле (8.23) под отдельным ракурсом могут на-

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

При этом не требуется каких-либо данных о других ракурсах.

Аналогично вклад в величину оценки /*(г, ф) каждой свернутой проек-

ции p

(nd

тА) для произвольного ракурса (при этом т фиксировано, л из-

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

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

считать значения p

[r cos(mA - ф)

тА] по величинам p

{nd, тА) при

-ДГ ^ л ^ N.

Таким образом, на практике операцию свертки, следующую за опера-

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

ординат (г, ф), при которых берется оценка /*(г, ф). Вместо этого может

быть предложен следующий алгоритм вычислений на ЭВМ.

Сначала задается последовательность f

... ,f

дискретизиро-

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

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

;

, ф.). Итерационный

процесс начинают с оценки

fo(rj,

<Pj)

= 0, (8.24)

при 1 ^ j ^ J. Для каждого т, лежащего в интервале О^/я^Л/- 1, из

m-й картины мы получаем (т + 1)-ю картинку с помощью следующего

двухступенчатого процесса:

а) Вычисление p

(n'd, тА) для

—

N ^ л' ^ N по значениям

p(nd

тА)

(—N ^ л ^ N), а также предварительный расчет значений q[(n'

—

n)d].

б) Вычисление для 1 ^

j'^

J значений

/m+ i0> Ф)) = fnfj-js ф

) + Ap

(r;COs(wA - фД гоД). (8.25)

Для расчетов по формуле (8.25) можно использовать процедуру интерполя-

ции данных.

Дискретизированное изображение /

, полученное в результате указанно-

го выше итерационного процесса, и является нашей оценкой /\ Заметим,

что,

как только величины f

m +

j(r , фЛ вычислены, потребность в значени-

ях/^'),

Фр отпадает, и поэтому в компьютере можно использовать одну и

ту же ячейку блока памяти для последовательной записи значений/

(г, фХ

Щ.Ф-i /•(/> Ф>

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

147

Отсюда читатель может сделать вывод

о том, что

описанный выше

способ вычислений более эффективен,

чем

применение

«в

лоб» формул

свертки

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

во

всех точках

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

(fp

Ф\)»

•••

(

<t>j)-

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

что нет

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

по

формуле (8.24)

раз

(т.е. по

одно-

му разу

для

каждого элиза

и для

каждого ракурса),

достаточно повто-

рить

эту

процедуру всего лишь

(2N + \)М раз (т.е. по

одному разу

для

каждого луча

и для

каждого ракурса). Поскольку обычно величина

одно-

го порядка

с (2N + I)

, то

этим достигается существенный выигрыш

объеме вычислений.

Единственный вопрос, который остался открытым

при

обсуждении

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

— это

вопрос

выборе сворачивающей функции.

Решение этого вопроса опирается

на

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

опе-

рации дискретизации исходных данных.

8.4. ПРЕОБРАЗОВАНИЕ ФУРЬЕ

Преобразование Фурье представляет собой оператор, который воздей-

ствует

на

функцию вещественного переменного,

результате чего получа-

ется комплексная функция

(т.е.

функция, представимая

форме

а + /£, где

—

вещественные функции,

а / =

V—

1).

частности,

мы

будем

ис-

пользовать запись ехр[/а]

для

обозначения комплексной величины cos а

sin а, поскольку любое комплексное число

а + ib

можно записать

экс-

поненциальной форме

ехр[/а],

где г

выбирают равным положительному

корню Va

+ b

а а

находят

из

соотношения

a tga = b. При

этом

назы-

вают модулем,

а a —

аргументом комплексного числа

а + ib.

Модуль

обычно обозначается

виде

\а + ib\, что

совпадает

обозначением абсо-

лютной величины

\а\

вещественного числа

а.

Если

ф —

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

то

ф(и) можно определить

виде

ф(и)

du = (

а(и)

duj + i(

Ь(и)

du),

(8.26)

где

а и b

— вещественные функции,

для

которых справедливо равенство

а(и)

+ ib(u) - ф(и).

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

зывает комплекснозначную функцию

вещественного переменного

дру-

гой комплекснозначной функцией

Уф

(также вещественного переменного)

помощью соотношения

/•00

[&ф](и)= <ttu)exp[-2niUu]du.

(8,27)

<х>

Аналогично одномерное обратное преобразование Фурье представляет

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

веществен-