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

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

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

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

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

изображениям г, д и е соответственно. Значения порога и весовых коэффиииентов при сглажива-

нии приведены в тексте (разд. 11.4).

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

199

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

точному хронометражу, поскольку, как было указано в разд. 6.5, затраты

машинного времени зависят от многих факторов, не обязательно имеющих

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

ции.

Для подтверждения этого положения мы применили к стандартным

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

мини-ЭВМ типа Eclipse

S/200

(машины, изготовленной фирмой Data Ge-

neral для фирмы General Electric и снабженной специальным томографиче-

ским процессором). Реконструкция потребовала лишь 90 с, что намного

меньше величины, приведенной в табл. 10.2 (847 с) и относящейся к вычис-

лениям на большой ЭВМ типа Cyber 173 с использованием системы про-

грамм SNARK 77. Поскольку процессоры обеспечивают высокое быстро-

действие, то для наших целей можно создать недорогой специальный то-

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

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

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

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

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

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

больших объемах), что их запоминание для осуществления операции по-

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

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

рис. 3.3, е (разд. 3.4.). При этом сверточный алгоритм для веерного пучка,

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

но,

оказывается проще для реализации с использованием более быстро-

действующих вычислительных устройств.

Выбор функции «окна» дает существенную разницу в качестве рекон-

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

бражение является все же «лучшим». Значение параметра а = 1,0 соот-

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

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

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

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

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

установку уровней черного и белого). Изображения внутричерепной обла-

сти при а = 0,8 и а = 0,54 выглядят лучше, чем изображение при а = 1,0,

а в целом не отличаются друг от друга. При нелинейном сглаживании был

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

фактическому усилению отдельных артефактов при а = 1,0 и в меньшей

степени при а = 0,8. Из всего этого можно сделать вывод, что выбор

функции «окна» и нелинейного сглаживания зависит от того, какая инфор-

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

200

ГЛАВА 10

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

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

операции повторного разбиения), настолько это известно автору, даио в

[127].

Эти

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

время спустя ие появилась работа (100).

Подробный вывод формулы обратного преобразования Радона для веерного

пучка (10.3), а также (10.9) дай в [77].

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

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

[29],

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

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

[153],

в которой рассмотрена

функция точечного отклика). В [29] приведены также вопросы, ие рассматриваемые

в данной книге, например вопрос о нахождении оптимальной функции «окна». Мате-

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

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

[145].

Преимущества использования метода линейной интерполяции по сравнению с

интерполяцией по ближайшим значениям были показаны, например, в работе [77].

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

на работе [76], в которой приведена библиография более ранних работ, иа которых

она сама базировалась (например, [40]).

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

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

0,01 с [46].

Алгебраические алгоритмы

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

В этой, а также в двух последующих главах будут рассмотрены алго-

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

ния.

Алгебраические алгоритмы реконструкции образуют большую группу

алгоритмов реконструкции. Их название исторически сложилось абсолют-

но случайно: в них нет чего-либо более «алгебраического», чем в методе,

который мы рассмотрим в следующей главе. Характерные особенности ал-

гебраических алгоритмов реконструкции требуют тщательного анализа, ко-

торый производится ниже.

11.1.

ЧТО ТАКОЕ АЛГЕБРАИЧЕСКИЙ

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

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

значенные для дискретной реконструкции изображения. Как уже отмеча-

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

изображения дг, для которого справедливо соотношение у = Rx + е

где

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

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

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

типа тех, которые рассматривались в разд. 6.4. Запись

дг*

принята для обо-

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

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

итерационными. Другими словами, с их помощью получают такую после-

довательность векторов х*°\ х<

\ х*

\ ... , которая сходится к оценке дг*.

Это означает, что для 1 ^ j < J имеется гарантия того, что величина

х*-

0-я компонента вектора на £-й итерации) произвольно близка к х* при ус-

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

дг**

из х<*) называется шагом итерации.

Величины х*

к + !)

в алгебраических алгоритмах получают из дг

(Аг)

, анали-

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

(6.23)].

На практике упо-

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

значим через i

величину [A:(mod/) + 1], т.е. #

= 1, i

= 2, ... ,/,_, = /,

// = 1, iy

j = 2 а через г, — У-мерный вектор-столбец (1 < # < /),

у-компонента которого равна /у,

. Другими словами, /} представляет собой

транспонированную /-ю строку R.

202

ГЛАВА II

Каждый к-й шаг итерации в алгебраических алгоритмах можно описать

функцией о^, аргументом которой является /-мерный вектор и которая со-

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

ным вектором (на языке математики это можно записать как а^: R

х R

—

, где R — множество вещественных чисел). Кроме того,

{k+

= a

). (ll.i)

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

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

, а

, .... Чтобы полу-

чить (к + 1)-ю итерацию, мы обращаемся к к-й итерации,

-Vi

строке про-

екционной матрицы R и /

-й компоненте вектора измерений у. Подобный

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

вектор х**

'* можно запомнить в том же блоке памяти ЭВМ, в который

занесены значения дг

(Л)

и необходимость в которых отпадает после k-ft ите-

рации. Отметим, что изложенный в конце разд. 8.3 сверточный алгоритм

также обладает свойством эффективности по объему памяти, чего нельзя

сказать об итерационных процедурах, описанных в следующей главе.

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

бом выбора последовательности функций о^. Ниже мы проиллюстрируем

предыдущий анализ на очень простом примере.

Один из способов выбора функций а

таков: для /-мерных векторов х и

/ и для любого действительного числа z имеем

(\\

Г,

г) =

* +

при

<'-

О * о,

при <Г, О = 0. (11.2)

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

торов, т.е.

<*,*> = I'j-V U1.3)

Заметим, что в данном случае а^ одна и та же для всех к.

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

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

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

<^, r

) Ф 0, то имеет место соотношение вида

1= 1

т.е. i

-t приближенное равенство становится точным после А>го шага ите-

рации, что следует из соотношений (11.1) — (11.3). Правую часть выраже-

ния (11.4) можно записать в следующем виде:

АЛГЕБРАИЧЕСКИЕ АЛГОРИТМЫ РЕКОНСТРУКЦИИ ИЗОБРАЖЕНИЙ 203

= (г-

(k)s

Уйк

~ ^

^ <Г Г■ > = V

что и утверждается в (11.4).

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

с по-

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

чений

л***

становится весьма простой: мы просто умножаем *<*) на вектор

. На практике данная процедура обновления значений x

ik)

с точки зрения

вычислений весьма экономична.

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

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

(6.17)].

При

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

Этот вывод имеет два следствия. Во-первых, большинство компонент век-

тора r

равно нулю. При разбиении изображения на ( х t элизов самое

большее (2(

—

1) элизов может иметь пересечения с прямой линией. Таким

образом, из (

компонент вектора r

только не более (2f — 1) компонент

(а обычно только порядка f ) имеют ненулевые значения. Во-вторых, по-

ложение и величина ненулевых компонент г

1к

легко можно определить, ис-

ходя из относительного геометрического расположения /^-го луча относи-

тельно (t х ()-сетки. Таким образом, нет необходимости запоминать в

ЭВМ всю проекционную матрицу /?. За один цикл вычислений необходимо

иметь лишь одну строку этой матрицы, всю информацию с которой легко

рассчитать.

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

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

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

Предположим, что имеется набор чисел у,, ...

, таких, что / = 0, ес-

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

(11.3) потребует только V операций умножения, причем U < J. Аналогич-

но величину </, /> можно рассчитать за V операций умножения. Имея вы-

численное значение [(z - </, *»/</, />], выполненное с помощью 2V опера-

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

производиться также путем последующих U операций умножения. Это

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

которых у = j

при некотором и (1 ^ и ^ £/), а изменение требует допол-

нительного умножения х на фиксированную величину /..

Можно видеть [формула

(11.2)],

что отдельные шаги вычислений доста-

точно просты для эффективной реализации на ЭВМ.

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

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

итерации х^

к + !)

из х

{к)

. Предположим, что в дополнение к требованию

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

дующие условия на способ проведения А-й итерации.

204

ГЛАВА 11

а) только те элизы, которые пересекаются

I'^-M

лучом, должны иметь

измененные значения плотности;

б) изменение значения плотности в элизе должно быть пропорциональ-

ным разности у

1к

—

ifc

, дг***). т.е. «ошибке» в i

-M приближенном ра-

венстве перед Ar-й итерацией;

в) изменение

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

кшГ

Указанные условия, аналогичные условиям, налагаемым на операцию

дискретного обратного проецирования (разд. 7.3), дают единственный спо-

соб определения функций и приводят к соотношению (11.2). В ранних рабо-

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

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

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

Прежде чем перейти к детальному описанию конкретных алгебраиче-

ских алгоритмов, сделаем два замечания.

Во-первых, для точного определения последовательности х®\ х*

х&\

... необходимо точно задать вектор начальных значений х*°\ выбор

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

ритмов. Более подробно на этом вопросе мы остановимся ниже.

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

струкции, связанный с решением системы неравенств Nx ^ д (разд. 6.4), то

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

отношением (11.1), не всегда является адекватной. В этом случае может по-

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

принципе не изменит вычислительных требований. Этот вопрос будет де-

тально обсужден в разд. 11.3.

11.2.

РЕЛАКСАЦИОННЫЕ МЕТОДЫ РЕШЕНИЯ СИСТЕМ

УРАВНЕНИЙ И НЕРАВЕНСТВ

В данном разделе дано математическое обоснование алгебраических ал-

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

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

равенств. Другими словами, необходимо найти У-мерный вектор х, для ко-

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

<л„х> <q

1 <i < />, (11.5)

где л. — заданные 7-мерные векторы, a q

— заданные вещественные чис-

ла.

Неравенства (11.5) взяты из соотношения (6.42).

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

мере одну ненулевую компоненту. Физический смысл этого предположения

состоит в том, что во внимание не принимаются те измерения, в которых

не дает вклада ни один из элизов. Причина введения данного предположе-

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

чаев,

когда (л,., л

> = 0, аналогично тому, что было сделано в формуле

(112).

АЛГЕБРАИЧЕСКИЕ АЛГОРИТМЫ РЕКОНСТРУКЦИИ ИЗОБРАЖЕНИЙ 205

Введем некоторое множество векторов /V (1 ^ / ^ Я) и величину N.

Для любого i в интервале 1 ^ / ^ Р имеем

= {х|<п„х><^} (11.6)

ftNt.

(И-7)

Другими словами,

Щ —

множество векторов, удовлетворяющих / нера-

венствам из общего числа Р неравенств в выражении (11.5), а N— мно-

жество векторов, полностью удовлетворяющих всем Р неравенствам. Зада-

ча при этом заключается в нахождении одного элемента множества N.

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

заданным значениям л,, ... , п

и д

... , д

находить х из множества N.

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

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

(

. _ jx, если <г, х> < z,

* *' '

~ \х + K

\(z -

<Г,

х»/<Г. f>]f в противном случае, (11.8)

где

Х**>

— вещественное число, которое в контексте данной книги названо

релаксационным

параметром.

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

релаксацион-

ным методом

решения системы неравенств

(или просто релаксационным

методом для неравенств), а именно

<0)

произвольное,

k+l)

= <x

\ (11.9)

где функции а

определяются соотношением (11.8), a i

= [Ar(modP) + 1].

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

вида

0< e

< А<*> < е

< 2. (11.10)

для всех к

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

ность величин х^°\ л*

*, х*

\ ... , которые сходятся к значению вектора во

множестве N при единственном условии, что множество N непустое. Дока-

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

Теперь дадим геометрическую интерпретацию релаксационного метода

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

онных параметров.

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

Я,

= {х|<Я,х> = <?<}. (11.11)

206

ГЛАВА 11

Рис. 11.1. Рисунок, демонстрируюший применение релаксационного метода для

частного случая, когда / = J = 2(Х<*> = 1 для всех к).

При этом значения я, и л

выбираются в соответствии с выражением (11.12), a q

и q

— в соот-

ветствии с выражением (11.13). (Рисунок заимствован из работы [71] с разрешения авторов.

есть множество векторов

JC,

для которых /'-е неравенство переходит в точ-

ное равенство. Заметим, что H

является подмножеством множества N

Каждая из величин f/ в математике называется гиперплоскостью. При

J = 3 гиперплоскость переходит в обыкновенную плоскость в обычном

трехмерном пространстве, а при 7=2 — в прямую линию на плоскости.

Указанный двумерный случай иллюстрируется рис. 11.1.

Здесь мы имеем две гиперплоскости Н

}

и //

для которых

«.=[']■

-,«[]].

01.12,

<?,

24. q

= 30.

(1113)

АЛГЕБРАИЧЕСКИЕ АЛГОРИТМЫ РЕКОНСТРУКЦИИ ИЗОБРАЖЕНИЙ 207

Установим следующий очевидный геометрический факт: вектор n

(пред-

ставленный в виде линии, соединяющей на плоскости начало координат с

точкой п

(

) перпендикулярен линии Я. Этот вывод можно обобщить на

случай пространств произвольной размерности. Читатель, знакомый с ос*

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

ношения (11.11).

Поскольку в релаксационном методе для неравенств, описываемом вы-

ражением (11.9), роль параметра / в формуле (11.8) выполняет

ik>

то ока-

зывается, что если итерации х^

к + !)

и х^ отличаются друг от друга, то

разность

JC**

+ !

> - х^

к)

(которая умножается на n

) представляет собой

вектор, перпендикулярный гиперплоскости H

Дадим теперь геометрическую интерпретацию величинам Л^. В двумер-

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

Л^

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

одну сторону от линии H

а также на ней самой. Например, на рис. 11.1 Л^

и N

— полуплоскости, включающие в себя и начало координат. Их пере-

крывающаяся часть заштрихована на рис. 11.1. Подобное геометрическое

представление также достаточно общо, и поэтому можно говорить о том,

что

Л^

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

по одну сторону от гиперплоскости Н^ а также на ней самой. Величина 7V,

как это определено соотношением (11.7), задается пересечением двух полу-

плоскостей.

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

ний (11.8) и (11.9), оценка х^

к)

не меняется на /г-й итерации (т.е. x

ik + !)

= х***), если

л***

принадлежит полуплоскости N

. Если х^

к)

не принадле-

жит полуплоскости N

, то оценка «перемещается» перпендикулярно гра-

ничной полуплоскости Hj

(т.е. отрезок х^

к + !)

—

х^

к)

перпендикулярен

а степень ее «перемещения» зависит от величины параметра \

{к

Если Х

(

*> = 1, то в этом случае, исходя из соображений, аналогичных

учитываемым в соотношении (11.3), можно показать, что

ikt

k + !)

> =

* т.е. значение х^

к + ])

лежит на гиперплоскости H

. Достаточно

просто можно доказать и другие положения. На Аг-й итерации в релаксаци-

онном методе для неравенств при х^

\ не принадлежащем полуплоскости

i/(9

перемещение оценки от х^ до х^

к + !)

происходит вдоль перпендикуля-

ра,

опущенного на гиперплоскость H

, и обладает следующими геометри-

ческими свойствами:

а) если К

{к)

< 0, то оценка удаляется от Н,

б) если Х

(

*> = 0, то перемещение отсутствует;

в) если 0 < Х

(Аг)

< 1, то перемещение оценки происходит в направлении

на гиперплоскость H

, но не достигает ее;

г) если Х

(Л)

= 1, то перемещение происходит до точного совпадения

оценки с H

д) если 1 < Х<*> < 2, то перемещение продолжается и после достижения

оценкой гиперплоскости //

#А

, однако оценка х^

к ч 1)

оказывается ближе к

, чем предыдущая оценка х

{к)

;

е) если Х

(

*> = 2, то оценка х

{к + 1)

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

х(к) относительно гиперплоскости /У

/А

;