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

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

228

ГЛАВА И

предложен в [73], а алгоритм реконструкции изображения объекта, имеющего лишь

два уровня плотности, дан в работе [62]. Любопытная демонстрация возможностей

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

нием дополнительной матрицы (не обсуждавшейся в данной книге) может рассмат-

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

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

Интересное теоретическое исследование о порядке, в котором необходимо брать

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

опубликовано в [58]. Данная работа и библиография к ие& интересна также в том

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

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

ций,

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

звана в [71] непрерывным алгебраическим алгоритмом реконструкции.

С другими модификациями алгебраических алгоритмов реконструкции, которым

не нашлось места в данной книге, читатель может ознакомиться в работах [50, 71].

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

конструкции, поскольку, как это было доказано в

[105],

они максимизируют энтро-

пию [формула

(6.44)].

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

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

описана в работе

[121].

Алгоритмы реконструкции

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

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

предназначенные

для

минимизации квадратной функции общего вида

включающие

себя

как

частные случаи

ряд

различных оптимизационных

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

в гл. 4.

Указанные процедуры отличаются

по

своей

сущности

от

алгебраических алгоритмов реконструкций

и в

литературе

носят название «алгоритмы реконструкции

одновременными

итерациями» (SIRT).

12.1.

МАТЕМАТИЧЕСКОЕ ВВЕДЕНИЕ

Вначале напомним некоторые понятия. Говорят,

что

матрица

симметрична, если

она

имеет квадратную форму

и для

всех

ее

элементов

nijj выполняется условие

- = т^

Говорят также,

что

матрица

положительно определена, если

она

симметрична

и,

кроме того,

для

любого вектора (столбца)

дг,

имеющего

по

крайней мере одну ненулевую

компоненту, выполняется условие х

Мх

> 0.

Заметим,

что для

каждой

положительно определенной матрицы

существует обратная

ей

матрица

А/~

. Если матрица симметрична

и для

любого вектора

имеет место

условие х

Мх

то

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

Проблему, которую

мы

собираемся разрешить, можно сформулировать

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

D и

неотрицательно определенные матрицы XV

и W

необходимо найти

вектор дг, который минимизирует величину

на

множестве

К, т.е.

\\0'

х\\

(12.1)

где

{х\к(х) минимально

(12.2)

при условии

к(х) = (у - RxfWtto - Rx) + (х - x

)

(x - х

(12.3)

Отметим,

что

правая часть соотношения (12.3) идентична функции

(6.39),

условие минимальности которой совпадает

обобщенным

критерием квадратичной оптимизации, если

230

ГЛАВА 12

= аА

(12.4)

= hB + сС

. (12.5)

Выражение (6.39) представляло собой обобщение соотношений

(6.33) — (6.35), (6.37),

(6.38). Теперь мы рассмотрим величины W

и х

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

ниже табл. 12.1 £/

и Uj являются единичными матрицами

соответствующих размеров, 6

и 9

— матрицами соответствующих

размеров с нулевыми элементами, 9 — вектором со всеми нулевыми, а

х — ненулевыми компонентами х.

В гл. 6 детально не обсуждалась сущность фигурирующей в

соотношении (6.38) величины В; теперь этот недостаток будет исправлен.

Хотелось бы только отметить, что В — неотрицательно определенная

матрица, откуда следует, что W

во всех выражениях [за исключением

(6.34),

когда W

= 9

] представляет собой положительно определенную

матрицу. Вот почему предполагалось, что W

в формуле (12.3)

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

матрицу с нулевыми элементами.

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

оптимизации к решению совместной системы уравнений [пример решения

подобной задачи встречался в разд. 11.3, в котором минимизируемая нами

функция действительно представляла собой частный случай функции (12.3)

при W

= Up x

= ц

]. Мы здесь исследуем два случая, когда

матрица W

положительно определена и когда W

— нулевая матрица.

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

определена

в соответствии с условием (12.5). Причины, по которым мы записываем

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

известна, однако критерий оптимальности относится к величине С""

Примером служит соотношение (6.33), в котором W

= ^^аИ^- кова-

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

С — матрица из большого числа J х J элементов (где / — число базисных

Таблица

12.1

Формула

(6.33)

у-\

»х

(6.34)

(6.35)

(6-37)

(6.38)

Л),

ЬВ + Uj

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

изображений), то в тех случаях, когда этого можно избежать, обращать

матрицу С нежелательно. Более того, необходимо стремиться

разработке

такого алгоритма минимизации к(х), при котором бы величина С

-1

не

использовалась.

Для разработки подобного алгоритма предположим, что С — положи-

тельно определенная матрица (с > 0), а В — неотрицательно определенная

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

таблице, за исключением соотношения (6.34), в котором W

= 9

данные

предположения справедливы]. Из линейной алгебры хорошо известно, что

неотрицательно определенная матрица имеет свой

квадратный

корень, т.е.

существует матрица С

,/2

, такая, что

i/2

= с (12 6)

Кроме того, матрица С

1/2

оказывается также положительно определенной;

поэтому существует и обратная ей матрица

С~

,/2

причем

С"

1/2

С~

1/2

«С

-1

Нам нет необходимости вычислять матрицу С

1/2

в рассматрива-

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

которых базируются алгоритмы.

Введем новую векторную переменную и с помощью следующего

соотношения:

и = С

1/2

(х - х

). (12.7)

Тогда

х = С

112

и + х

, (12.8)

и из выражений (12.3) — (12.5) следует, что

/с(х) = к(С

1/2

и 4- х

)

= аЬ ~ R(C

ll2

u + х

)]

Л1> - R(C

l/2

u + х

)]

(С

112

и)

(ЬВ + сС-

ХС

112

и)

= u

laC

1/2

ARC

112

+ ЬС

ВС

+ cU

Ли

-2au

ll2

A(y - Rx

)

+ а(у - Rx

)

A(y - Rx

). (12.9)

Вводя обозначения

Р =

l/2

ARC

l/2

+ ЬС

1,2

ВС

1/2

4- cUj, (12.10)

z = aC

ll2

A(y - Rx

), (12.11)

Ци) = WPu - u

z, <12 12)

замечаем, что при выполнении условия (12.8)

232

ГЛАВА 12

) = ifc(x) - НУ - R*o)

A(y - Ях

). (12.13)

Отсюда следует, что величина

х*

минимизирует^*) тогда и только тогда,

когда х* =

1/2

и*

+ х

, где w* минимизирует Л(м). Таким образом,

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

выражении (12.2), можно свести к нахождению значений м, которые

минимизируют функцию

в соотношении (12.12).

Рассмотрим далее случай, когда fV

= ©у, и предположим, что а Ф О;

при W

- 6

минимизация никогда не достигается. Вновь введем

переменную м, на этот раз определенную как

u^D^x.

(12.14)

Тогда

х =

DM,

(12.15)

и из условий (12.3) и (12.4) следует, что

к(х) = k(Du)

= а(у- RDu)

A(y - RDu)

= u

(aDR

ARD)u - 2au

Ay + av

Ay. (12.16)

Вводя обозначения

Р = DR

ARD, (12.17)

z = DR

Ay, (12.18)

и используя определение (12.12) для функции И(и), замечаем, что если

условие (12.15) выполняется, то справедливо также соотношение

h(u) = (\/2а)к(х) - ±у

Ау. (12.19)

Из него следует, что

минимизирует функцию к(х) тогда и только тогда,

когда х* = Du*

где и* минимизирует функцию Л(м). Поэтому в данном

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

согласно (12.2),

можно свести к нахождению значений м, минимизирующих функцию h в

выражении (12.12).

В обоих рассмотренных случаях мы получаем из (12.12) систему

уравнений, используя при этом следующий результат. Для любой

неотрицательно определенной матрицы Р размером

и произвольного

У-мерного вектора z вектор и минимизирует функцию (12.12) тогда и

только, тогда, когда

Ри = z. (12.20)

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

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

выше случаях, поскольку матрица Р [определенная либо выражением

(12.10),

либо (12.17)] является неотрицательно определенной, как это

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

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

что если система уравнений (12.20) не имеет решений, то множество К

пустое. Однако в рассматриваемых нами случаях существование решения

(12.20) гарантируется. Если Р определена соотношением (12.10), то

р — положительно определенная матрица и поэтому имеет обратную

матрицу. Немногим более сложно показать (детали мы опускаем), что

система уравнений (12.20) также имеет решение, если Р и z определяются

соотношениями (12.17) и (12.18) соответственно.

Как уже отмечалось, матрица W

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

Р — положительно определенная матрица; таким образом, при этом

имеется единственное значение и = P~

z, минимизирующее функцию h(u)

и, следовательно, единственное значение дг*, которое минимизирует к(х).

Данное значение дг* образует единственный элемент множества К, и

поэтому вторичный критерий оптимальности (12.1) становится излишним.

Предположим, что и* — решение с минимальной нормой системы

уравнений (12.20) в случае, когда W

= 6

. Тогда среди векторов,

минимизирующих функцию /?(")> вектор и* обладает наименьшей нормой.

Определив величину дг* с помощью соотношения

дг*

= Du

, замечаем, что

среди векторов дг, минимизирующих функцию к(х), вектор дг* обладает

наименьшим значением Ш~

х*1 = lw*l [см. комментарии, относящиеся к

соотношению (12.19)]. Таким образом, мы считаем, что вектор

дг*

найден.

Подводя итоги, можно отметить, что проблема квадратичной оптими-

зации, описываемая соотношениями (12.1) — (12.3), сведена здесь к

проблеме решения системы уравнений (12.20). Если W

— положительно

определенная матрица, то система (12.20) имеет единственное решение, а

если W

= В, то при этом потребуется решение с минимальной нормой.

В следующем разделе будет обсуждаться -класс алгоритмов решения

системы уравнений (12.20).

12.2.

МЕТОД РИЧАРДСОНА

РЕШЕНИЯ СИСТЕМ УРАВНЕНИЙ

В разд. 11.2 уже рассматривался метод нахождения решения о

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

безусловно, применим и к системе (12.20). В данном разделе будут

рассмотрены и некоторые другие методы.

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

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

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

релаксации, равной единице, удовлетворение одного из уравнений может

234 ГЛАВА 12

приводить

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

соответствующего данному уравнению. Набор подобных полос будет

визуально восприниматься

как

сильно зашумленное изображение, похожее

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

на

рис. 11.2,г.

алгебраических алгоритмах реконструкции

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

параметра релаксации,

а в

альтернативном подходе

это

будет достигаться

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

по

всем лучам одновременно.

По

этой

причине последний метод считаете*; относящимся

алгоритмам

реконструкции

одновременными итерациями (SIRT). Различия между

SIRT

подобными

им

алгоритмами касаются деталей, которые ввиду

отсутствия места здесь

не

рассматриваются.

Одним

из

способов коррекции итерации

для

одновременной коррекции

погрешностей

во

всех уравнениях является итерация вида

„<*+!>

= „«' +

Я<*>(г

Р//<*>). (12.21)

Отметим,

что

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

случае, когда

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

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

размерности

z и Ри^

будут отличаться

от

размерности

и^

Итерационный метод,

котором

к-я

итерация осуществляется

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

выражением (12.21),

численном анализе получил название

метода Ричардсона. Прежде

чем

исследовать указанный метод

на

сходимость, познакомимся

с ним на

двух примерах, рассмотренных

предыдущем разделе,

что

позволит

нам

получить представление

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

проблеме реконструкции изобра-

жений.

Если матрица

положительно определенная,

то

матрицы

Р и z

задаются выражениями (12.10)

(12.11) соответственно. При

к ^ 0

опреде-

лим

как

<*>

= С

'Ч

к)

+ х

Ц2.22)

Тогда

из

соотношений (12.10), (12.11)

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

что

{k+u

i:2

ik+u

+ x

= с

V" + х"Ъ - JV*')] + *о

к)

+ х

+ A<*»[C"

z - С

РС- "V

х"" +

lk)

[aCR

A(y-

)

(aCR

+ ЬСВ +

cUj)(x

ik)

- x

)].

откуда

<*+i»

<*>

x»*[aCR

A(y

+ (bcB +

cUj)(x

v"")]-

(12.24)

-*o)]

(12.23)

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

Из приведенного

предыдущем разделе анализа,

частности

из

коммента-

риев к соотношению (12.13), следует, что если последовательность

i/°\ и

(1)

и&\

. . .

сходится

решению системы уравнений (12.20),

то

последовательность

х®\ х*

*, х^

сходится

единственному значению

которое минимизирует функцию

к{х).

Отметим,.

что

алгоритм,

использующий итерационную процедуру вида (12.24),

не

содержит

величины

и, а

непосредственно формирует последовательность векторов

х**\ которая сходится

(по

предположению)

искомому вектору

изображения

х*.

Отметим также,

что

данный алгоритм

не

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

-1

или С

,/2

несмотря

на то что

С~

содержится

функции

к(х), а в

процессе математического

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

,/2

Для демонстрации справедливости рассмотренной выше процедуры

рассмотрим функцию

к(х),

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

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

выражением

(11.26).

При

этом

а = Д Ь = 0, с = 1, А = U

В = G

, С = Uj и

= ц

> а

соотношение (12.24) принимает

вид

(*+п

(М

R\y

ik}

)

+ (/i

{k)

(12.25)

Данной итерационной процедуре можно дать прямую интерпретацию,

исходя

из

теории реконструкции изображений. Текущая оценка вектора

изображения дг<*> изменяется

при его

сложении

другим 7-мерным

вектором, имеющим

две

компоненты, одна

из

которых пропорциональна

величине

{y - Rx^

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

проекцией (разд.

7.3)

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

проекцией,

связанной

текущей, Ar-й оценкой.

При

обратном проецировании указанной

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

для

коррекции текущей, существующей оценки

так,

чтобы

ее

проекция

стремилась

зарегистрированной проекции. Вторая компонента х-вектора

равна

—

х^ и

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

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

и его

текущим значением, которое можно

использовать

для

получения более близкой

наблюдаемой

эксперименте

величине оценки вектора изображения перед

его

измерениями. Относитель-

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

г,

физический

смысл которого будет раскрыт

разд.

11.3.

Таким образом,

из

приведенного анализа следует,

что

метод Ричардсона

этом случае дает

итерационную процедуру, функционирующую достаточно понятным обра-

зом.

В случае если

= 0

то Р и z

задаются

помощью выражении

(12.17)

(12.18) соответственно.

При к ^ 0

определим х^

к)

как

{к)

= Du

П2.26)

Тогда

из

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

(12.18),

также (12.21) получаем

236

ГЛАВА

= D[u

ik)

Ри

{к)

У]

= Du

{k)

+ А

(М

Л>> - D

ARDu

{k)

l (12.27)

откуда

<*+ о

<*>

дю^дт^ -

{k)

(12.28)

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

последовательность и®\ м*

\ и&\ . . . сходится к решению уравнения с

минимальной нормой Л(м), то последовательность х®\ х^

\ х&\ . . . схо-

дится к оценке x*

которая минимизирует функцию к(х) и которая такова,

что для любого другого значения дг, минимизирующего к(х)

выполняется

условие lD~

x*l ^

lD~

xl.

Отметим, что и в этом случае алгоритм,

включающий в себя процедуру итерации вида (12.28), не содержит величин

и, а непосредственно формирует последовательность векторов

изображений хР\ которые по предположению сходятся к искомому

вектору изображения.

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

Ричардсона к проблеме реконструкции изображения, необходимо вернуться

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

и№ и такую последовательность параметров релаксации Х**\ чтобы с

помощью соотношения (12.21) получить последовательность векторов,

которая сходилась бы к решению с минимальной нормой системы

уравнений Ри

—

Имеется целый ряд способов выбора параметров Х<*\ обеспечивающих

сходимость (12.21) к решению системы уравнений (12.20). Ниже мы

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

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

выходит за рамки данной книги. В конце раздела будет рассмотрен

немного более сложный способ подоЧюа параметров Х

(

*\ кроме одного,

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

решению за конечное число итераций.

Вещественная величина р называется собственным значением матрицы

Р, если существует такой ненулевой вектор и, что

Ри = ри. (12.29)

Нетрудно заметить, что неотрицательно определенная матрица имеет

лишь неотрицательные собственные значения. Величины ^макс^

Рмш/*

мы

будем использовать для обозначения наибольших и наименьших

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

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

определенной матрицы Р существуют стандартные методы.

Рассмотрим один из них. Для любого значения «

<0)

с помощью

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

алгоритма (12.21) формируется последовательность и

\ м

\ ir

\ . . .,

которая сходится к решению системы уравнений (12.20) при условии, что

для всех к ^ 0 имеем

(

*> = К (12.30)

где

0 < X <

2/ip^^).

(12.31)

Таким образом, выбор «оптимального» параметра

ведется по формуле

X = l/b^ + РминР). (12.32)

Полученное соотношение в сочетании с теоремой о минимальной

норме, приведенной в разд. 11.2, показывает, что если величина ы

(0)

выбирается в виде линейной комбинации из столбцов матрицы Р, то

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

(

определяемым с помощью соотношений (2.30) и (2.31), сходится

решению

с минимальной нормой системы уравнений Ри = z.

Несколько более сложен способ выбора параметров >S

\ если заранее

известно желаемое число проводимых итераций. Например, если нам

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

взять следующую последовательность параметров Х^:

= Z/fcWP

"мш/* - OWc* " Р

мт

Р*<*ЫЫк) -

1)/8)1,

(12.33)

нежели постоянное значение в формуле (2.32), причем

Р(0) = 1, р(1) = 4, р(2) = 2, р(3) = 3.

12.3.

СГЛАЖИВАЮЩИЕ МАТРИЦЫ

На практике матрицы В и С, фигурирующие в выражении (12.5),

представляют собой так называемые сглаживающие матрицы. В

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

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

Заметим, что обе матрицы В и С имеют размерность J х У, и поэтому

они отображают векторы изображения в векторы изображения. Предполо-

жим,

что эти векторы соответствуют (л х п)-разбиению изображения на

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

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

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

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

нечно большим. В этом случае значение плотности картинки в каждом эли-

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

и соседних элизах. Процедура сглаживания определяется тремя сглаживаю-

щими весовыми коэффициентами, обозначаемыми в разд. 11.4 через и>,, w

и*