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

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

Неитерационные алгоритмы

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

разложения функции в ряд

Рассмотренные в двух предыдущих главах алгоритмы реконструкции

изображений с разложением функций в ряд, а именно релаксационный ме-

тод и метод Ричардсона, содержат итерационные процедуры. Указанные

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

большинство обсуждавшихся в двух последних главах алгоритмов было

ориентировано на базисные изображения в виде матриц изображений.

Принципиально любые неитерационные алгоритмы можно применять к ре-

шению системы уравнений, которые возникают в задаче реконструкции

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

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

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

венно возможными с точки зрения реальных возможностей вычислитель-

ной техники.

В данной главе будут рассмотрены другие методы выбора системы ба-

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

но более простая система уравнений, и поэтому привлечение итерационных

алгоритмов решений здесь не требуется.

13.1.

РАЗЛОЖЕНИЕ ИЗОБРАЖЕНИЯ НА КОЛЬЦЕВЫЕ ГАРМОНИКИ

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

интересующие нас свойства матрицы проекций R.

Рассмотрим приведенный в разд. 12.1 критерий квадратичной оптими-

зации для случая, когда W, = С/, и И^

= 6

. Последний в точности совпа-

дает с методом наименьших квадратов [выражение

(6.34)].

Для сохранения

общности анализа предположим также, что D = £/

, хотя, как нетрудно

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

венные черты и для любой другой диагональной обращаемой матрицы D.

Как было показано в разд. 12.1, искомое решение х* является решением с

минимальной нормой для системы уравнений вида

= R

y; (13.1)

[соотношения (12.17), (12.18) и (12.20)]. Совместная система уравнений

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

249

(13.1) называется системой нормальных уравнений, связанной с задачей

минимизации нормы

Ну

- ДдгН

Поскольку матрица R

R имеет размеры У х У, то, даже если она и

имеет обратную матрицу, ее обращение представляет в общем случае весь-

ма сложную проблему. Обращение матриц с размерами

стандартны-

ми методами (например, методом

исключения

Гаусса) требует на ЭВМ

приблизительно У

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

на обратную матрицу с размерами

необходимо У

умножений.

В нашем

случае (при работе со стандартным фантомом,

где У

13 225)

проб-

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

следних главах были рассмотрены итерационные методы, дающие прием-

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

тому, что большинство элементов матрицы R равно нулю.

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

базисных функций (для соответствующих конфигураций системы регистра-

ции данных), чтобы матрица

R имела простую структуру и тем самым

позволяла существенно упростить вычисления.

Предположим, что имеются такие матрицы Р_

, **-у+ р • • .» Ру раз-

мером U х U [те

= U(2V + 1), а выбор обозначений станет понятен

ниже],

что R

R можно записать в виде

блочно-диагоналъной

матрицы,

т.е.

Р-.

®L

Р-,

®t

(13.2)

где О

—

нулевая матрица размером U х U. При этом обратная

R мат-

рица (конечно, если она существует) имеет вид

Р-1

О,

©i/

0с/

(13.3)

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

риц P

по отдельности, что требует всего

(2 К

+ 1)С/

/3 операций перемно-

жения. Если U « Vy, то при этом достигается существенная экономия вы-

числений. Кроме того, при этом для умножения вектора на обратную R

матрицу необходимы лишь (2V + \)U

операций умножения, которые

можно производить поблочно. В итоге можно сказать, что если

R пред-

ставлена в виде матрицы (13.2), то задачу можно свести к 2V + 1 частным

задачам и тем самым существенно уменьшить объем вычислений на ЭВМ.

Для иллюстрации этого утверждения нами был использован метод Гаусса,

250

ГЛАВА

однако данный вывод в принципе не зависит от способа решения системы

(13.1).

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

ставить матрицу R

R в форме (13.2), необходимо не только тщательно вы-

брать базисные функции изображений, но и сделать так, чтобы система ре-

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

мо,

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

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

лучей относительно общего центра. Отметим, что все описанные в разд.

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

ям.

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

кие положительные целые величины М и N, для которых справедливо вы-

ражение

/ = M(2N + 1) (13.4)

где

вп.

= 0

,о + т(2п/М1 (13.6)

а п и т — те единственные целые числа, которые в интервалах

-NKn^N,

0 <т< М - 1, (13.7)

и удовлетворяют соотношению

i = m + пМ + NM + 1. (13.8)

[ср.

с выражением

(6.8)].

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

ния,

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

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

ласти. В данном разделе будут описаны базисные изображения, принимаю-

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

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

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

гармоническими функциями угла

<р

(разд. 8.4, особенно та часть, которая

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

Проводя более детальный анализ, предположим, что d — положитель-

ная вещественная величина, характеризующая ширину кольца, а (У и

V — положительные целые числа, описывающие номер кольца и номер

пространственной гармоники соответственно. Выберем V таким, чтобы

Ud ^ Е [формула

(6.1)].

Пусть число базисных изображений J определяет-

ся как

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

251

J = U(2V + 1). (13.9)

Для 1 < j < J существует единственное значение м, лежащее в интер-

вале 0 ^ и < V - 1, а также единственное значение v (- V < t; < Г), для

которых справедлива формула

7 = и + wl/ + Ft/ + 1. (13.10)

Определим при 1 ^ у < J базисное изображение bj с помощью соотноше-

ния

Ъ{г, ф) =

(r)T

(<t>\

(13.11)

где для и, лежащего в интервале 0 < и < £/ - 1, имеем

ад-Я

i«*

*'<(« +и*

(13

12)

(0 для остальных г,

Кроме того, при

—

V ^ u ^ К

■С

cos(,0) при^О,

(13ЛЗ)

причем значения и и v определяются в соответствии с выражением (13.10).

Почему такая система базисных функций изображений является разум-

ной? Нам необходимо убедиться в том, что любое интересующее нас изо-

бражение можно аппроксимировать линейной комбинацией базисных изо-

бражений bj. Заметим, что, если выбрать величину d достаточно малой

для всех

<р

и и, лежащих в интервале 0 < и ^ U

—

1 и удовлетворяющих

условию ud < г < (и + \)d, разумным было бы предполагать функцию

/(г, ф) приблизительно равной

Ш) = /((« + iK

ф).

(13.14)

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

ция/ принимает постоянные значения вдоль линий, идущих радиально. Та-

ким образом, можно считать, что в м-м

кольце

/ равна функции f

с чисто

азимутальной зависимостью от

<р.

С учетом условия (13.12) при 0 < г < Е

оно приводит к следующей зависимости:

/(г,<£)-

{г)/

{ф).

(13.15)

Учитывая, что f

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

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

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

252

ГЛАВА U

ла функции Г

, определенных формулой (13.13). Усечение упомянутого ря-

да Фурье с сохранением гармоник v в интервале - V < v <

дает прибли-

женное значение функции /

. Комбинируя приближенные выражения для

fif*

Ф)

[соотношение (13.15)] и для усеченного ряда Фурье, получим аппрок-

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

функции

в виде линейной комбинации элементар-

ных изображений bj, определяемых соотношением (13.11).

В последующем анализе для любого целого числа, такого, что

1 < j < Л обозначим через

те единственные, которые лежат в ин-

тервалах 0 < Uj < С/

—

1, -К^^Ки удовлетворяют при и = Uj и

v -

соотношению (13.10). Там, где это не вызывает путаницы, индексу

будет опущен.

Для определения элементов г^ , матрицы проекций R необходимо вы-

числить значения &bj во всех точках ( t

пт

Для упрощения обозначений при -/V £ п ^ N введем функцию а

Л>

определенную как

farctg (z/

) = ) *

C-ir/2,

если

'„

*0;

если

<«

= 0, и г > 0;

если

= 0 и г < 0.

(13.16)

В данных обозначениях выражение (6.4) примет вил

№Ж„>

Ял.

J Wn + z

+ a„(z))

<fe.

(13.17)

Пусть для —N

O^NHO^U^

I/ имеем

(VM)

- /

, если /

< ud,

-Ч

если ud < /„.

(13.18)

С использованием соотношений (13.12) и (13.18) выражение (13.17) можно

переписать в виде

l*bjWn^n,J= f

TAO

+ aJLz))dz

+ I "' W„.

+ a„(z))dz. (13.19)

Дальнейшего упрощения данного выражения можно достичь с помощью

метода, весьма эффективного в наших исследованиях. Утверждается, что

имеет место соотношение вида

№]('„,

Ящт

)

= 2Т£0

Ячт

) f

"'"

cos(ra

(z)) dz. (13.20)

Важность этого соотношения будет» яснее, если при 0 ^ и ^ U

—

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

253

—

V < v < К и — N < и < N определить функцию С

„ в виде

С„>) = 2 Г'"*'со5(ш

(2)) dz. (13.21)

При этом соотношение (13.20) можно переписать в виде

rij=T

^y (13.22)

Таким образом, (i, у)-й элемент матрицы проекций можно представить в

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

чениями i и v (но не ы), а другой — значениями

и л (но не /и).

Доказательство того, что соотношение (13.20) является следствием

(13.19),

необходимо дать отдельно для случаев v > 0 и v ^ 0, однако ниже

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

случая аналогично. Итак, имеем

sin(i?0

v<x

(z))

dz = siitt^J cos(ixx

(z)) dz

+ cosO^J Г "' sin(i;a„(z)) dz. <

)

Рассуждая аналогично и замечая, что а„(—z) = — c*

(z), получаем

I sin(i;^

tni

+ va

(z)) dz = sin(t>0

) J

cos(v<x

(z))

-cos^J sin(r;a

(z)) rfz. (13.24)

Справедливость утверждения (13.20) доказывается путем сложения левых и

правых частей соотношений (13.23) и (13.24).

Получив выражение (13.22) для любого элемента

г,-

матрицы проекций

/?,

мы можем вычислить обобщенный элемент

матрицы R

т.е.

m = 0 л = -N

= £ G.

.>)G^n)(Y V-.«)V..jY (13.25)

«=-N \т = 0 /

Теперь получим следующее свойство ортогональности, элементарное,

но громоздкое доказательство которого здесь опущено.

254

ГЛАВА 13

Пусть W

—

произвольное положительное целое, а у — произвольные

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

и>

и y

со-

отношением вида

= у +

(2n/W). (13.26)

Пусть также L — произвольное целое число, а а и b — неотрицательные

целые числа, для которых а + b < W. Тогда

L+

w -

X sin ay

cos by„

X cos ay

cos by

и- = /.

L+H'- 1

X sin ay

sin />y

w = L

Из приведенных выше соотношений нетрудно видеть, что при М ^

^2K+1H-K^I;^K

имеет место равенство

м-

1(Ш,

))

=//,л/. оз.зо)

т = 0

где h

= 1/2, если v Ф Q

и h

= 1, если

= 0. При выполнении неравенст-

ва - V ^ v

= v

^ К имеем

м-

I ^(«и^Д,.J = 0. (13.31)

ги = 0

Отсюда следует, что матрица R

R представима в форме (13.2) в случае, ес-

ли число ракурсов М не превышает числа пространственных гармоник

IV + 1. Для величины v, лежащей в интервале —V^v^ V, матрица P

(13.2) имеет размеры U х С/, причем ее (к, ( )-й элемент имеет вид

<РД., = Л,М I G

.(/i)G^

,(n). (13.32)

п= -Л*

= 0, (13.27)

при Ф fc,

при= /7 и а Ф 0 (13.28)

при

= Ь = 0,

■{:

0 при а Ф b или а = fr = 0,

уу для остальных г,

(13.29)

Вычисление всех коэффициентов

)

представляет собой трудоемкую и

длительную процедуру для разумно выбранных значений U, V, М и N, од-

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

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

255

счет использования свойства симметрии матрицы P

, поскольку G

u v

(n) =

я С

(n) [выражение (13.21)]. В любом случае достигается одна и та же

цель — сведение первоначально сформулированной задачи решения систе-

мы J уравнений с J неизвестными (13.1) к задаче независимого решения

2V + 1 систем из U уравнений с U неизвестными.

13.2.

ПОЛИНОМИАЛЬНОЕ РАЗЛОЖЕНИЕ

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

носится к выводу соотношения (13.31), которое позволяет привести матри-

цу R

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

тогональности. Один из способов доказательства этого свойства состоит

в следующем: для фиксированных значений v и п полагаем, что V

v п

есть

Л/-мерный вектор с компонентами

)\ тогда соотношение*(13.31)

устанавливает, что если v

Ф v

, то векторы V

v t s

и V

v t r

ортогональны

при любых значениях s и t (т.е. их скалярное произведение равно нулю).

Рассмотрим выражения (13.25), (13.32) и предположим, что сумма вида

I Gu

(n)G

UqtVq

{n) (13.33)

и- N

также удовлетворяет условию ортогональности. Говоря более точно, бу-

дем предполагать, что всякий раз, как только v

= v

значения суммы

(13.33) отличны от нуля, если и Ф u

. В случае когда последнее утвержде-

ние справедливо, матрицу R

МОЖНО

считать диагональной, а нахожде-

ние обратной ей матрицы становится тривиальной задачей. В данном раз-

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

которые обеспечивают получение подобных результатов.

К сожалению, требуемое для этого математическое описание выходит

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

ражения. При рассмотрении данного вопроса мы будем придерживаться

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

ла.

Предположим, что U и V — положительные целые числа, причем, как и

раньше, V — номер пространственной гармоники, а смысл U будет поня-

тен ниже. Число J базисных функций равно U(2V + 1), а величины и и v

(точнее, Uj и t>.) определяются по формуле (13.10). Опишем базисное изо-

бражение Ь. (1 ^ j ^ J) соотношением вида

где при 0 ^ и ^ U - \ и -V ^ v ^ V величина T

определена выражени-

ем (13.13), и, кроме того,

256

ГЛАВА 13

*u»={o

для остальных

/£)

,1|

С|,|,и((^)

) при г < £,

(13.35)

где в свою очередь Q

]v] и

— полином w-й степени. Для полноты картины

приведем точную формулу для

и отбросим детали вывода, так как в

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

^^=i^-^fX

]

Y-

(,336)

Система базисных изображений (13.34) отличается от используемой в

(13.11) системы лишь выбором радиальной составляющей. В соотношении

(13.12) функция R

выбирается в форме «ящика» (прямоугольной функции),

принимающей значения 1 в некоторой небольшой области и 0 вне ее. В то

же время R

v и

по определению (13.35) является полиномом степени

+ 2и. Хотя последнее утверждение приведено без доказательства, чи-

татель, вероятно, знает о том, что любую

функцию

/ двух полярных коор-

динат, записанную в форме полинома степени w от переменных г cos

г sin ф, можно представить в виде линейной комбинации трех базисных

изображений bj

для которых выполняется условие \v\ + 1и ^ w. Други-

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

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

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

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

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

функций bj обладает свойством ортогональности в случае, когда регистра-

ция исходных данных производится по определенной программе. Этот во-

прос рассмотрен ниже.

Для произвольного значения f , лежащего в интервале

— Е

^ ( ^ £,

имеем

\*Ъ$/

в) = (2£/(Н 4- 2и + 1)Х1 - {'/E)

)

JL'tE)W\

(13.37)

где для любых целых неотрицательных значений w функция U

представ-

ляет собой полином Чебышева w-й степени второго рода, определяемый

выражением

17„(cos ф) = (sin(w + l)tf>)/sin ф. (13.38)

Выражение (13.37) приведено без вывода, однако необходимо отметить,

что запись полинома Q

в форме (13.36) была дана как следствие соот-

ношения (13.37).

На первый взгляд неясно, почему упомянутое соотношение (13.37) явля-

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

ортогональности векторов, од-

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

257

надо при анализе полинома (13.38) можно видеть, что при достаточно тща-

тельном выборе величины £

правая часть выражения (13.37) может су-

щественным образом упроститься.

Например, если в интервале — N ^ п ^ N выбрать

/„= Ecosy

, (13.39)

где величина у

определена формулой (13.26) при произвольном у и

W = 27V + 1, то

[ibJ(/„,

0) = (2E/(\v\ + 2и + l))sin((M + 2i/ + 1)у„Ш0). (13.40)

Вычисления по полученной формуле (13.40) для значений 0, задаваемых

п m

(—N ^ п ^ N, 0 ^ т ^ М

—

1), точно приводят нас к выражению

(13.22),

где С

и v

определено в виде

UtV

(n)

= (2£/(|г 2и + l))sin((|r| + 2и + 1)у

). (13.41)

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

нять, что матрица R

R имеет форму (13.2) при условии, что число ракур-

сов М не меньше числа пространственных гармоник 2V + 1. При усло-

вии что - V ^ v ^

матрица P

в (13.2) имеет размеры U х С/, каждый

(£, I )-й элемент которой задается соотношением (13.32), а величина

Mf v

— соотношением (13.41).

Однако соотношение (13.32) можно существенно упростить, если в него

подставить (13.41); тогда

4h,.ME

'*■'

| + 2А- 1)(|г| + 2/ - 1)

х £ sin((N + 2А - l)r„)sin((|f| + 2/ - 1)у„). (13.42)

Если теперь в формуле (13.29) положить W = 2Л/ + 1, L = -N,

а =

I и I

+ 2* - 1 и fc =

+ 2 е - 1 (1 ^ к ^ С/, 1 ^ ( < С/), то по-

лучим

I sin(((r| 4- 2* - Dr

)sin((|r| + If - \)у

)

-I

0 при / Ф /с, (13.43)

(2N + 1)/2 при / = /с,

при условии что 2\v\ + 2* + 2/

-2<27V + 1. Комбинируя выражения

(13.42) и (13.43), получим