Добавил:
Опубликованный материал нарушает ваши авторские права? Сообщите нам.
Вуз: Предмет: Файл:

Симонов Томографические измерителные информационные системы 2011

.pdf
Скачиваний:
16
Добавлен:
12.11.2022
Размер:
9.04 Mб
Скачать

Рассмотрим пример, когда для восстановления используется только три проекции (рис. 3.40), хотя для практических целей потребуется гораздо большее их число. Восстановление производится путем обратного проецирования каждой проекции через плоскость, т. е. величина физического сигнала, соответствующая данной луче-

вой сумме P(l,θ), прикладывается ко всем точкам, которые обра-

зуют этот луч. После того, как это сделано для всех проекций, получается приближенная аппроксимация исследуемого объекта. Так как в этом методе проекции как бы растягиваются обратно через восстанавливаемое сечение, метод и назван обратной проекцией. Для каждой точки изображения восстановленная μ(х,у) является суммой всех лучевых проекций, которые проходят эту точку. Поэтому метод обратной проекции называется методом суммирования или методом линейной суперпозиции.

P (l, θ)

2

Y Y

2

l

l

 

P (l, θ)

X

X

1

2

P (l, θ)

l

 

 

 

Рис. 3.40. Восстановление объекта методом обратной проекции: 1 – объект; 2 – его проекции; справа – восстановленный объект

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

М

 

μ(х, у) = P(x cosθj + y sin θj , θj )Δθj ,

(3.126)

j=1

281

где x cosθj + y sin θj = l в соответствии с (3.122) и суммирование

проводится

по

всем

углам

проекции

θj.

Аргумент

l = x cosθj + y sin θj

соответствуют только тем лучам, которые про-

ходят через точку (х,у), коэффициент Δθj представляет угловое рас-

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

щью уравнения (3.126), не идентична истинной μ . Как видно из

рис. 3.40, восстановленная картина содержит значительный ложный сигнал. Это происходит потому, что точки за пределами объекта получают часть интенсивности спроецированного обратно сигнала. Кроме того, точки внутри объекта получают интенсивность составляющих сигнала от соседних точек, в результате чего небольшие перепады μ не различаются.

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

Предположим, что объект исследования состоит из одной точки. Тогда результат восстановления по проекциям методом суммирования (3.126) будет представлять собой не точку, а многолучевую звезду, центр которой находится в восстанавливаемой точке

(рис. 3.41).

1

2

Рис. 3.41. Восстановление точечного объекта методом обратной проекции: 1 – объект; 2 – его проекция; справа – восстановленный объект

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

282

вать фон, пропорциональный 1/r, где r – расстояние от точки. Фон является основным источником погрешностей, которые уменьшают достоинства этого метода. Один из путей уменьшения этого фона является использование определенного фильтра, который бы убирал ложные сигналы от проекций.

В практических приложениях регистрируемые данные соответствуют оценкам функции P(l,θ) для обширной совокупности дискретных значений параметров l и θ, а реконструируемое изображение имеет вид двумерного массива чисел. Будем считать, что отсчеты функции P(l,θ) берутся равномерно по l и θ. Рассмотрим случай регистрации проекций для M значений угла θ, взятых с шагом Δθ, и N значений расстояния l, взятых с шагом l . Определим

целые числа N + и N следующим образом:

N

+

= (N 1) 2,

 

 

 

 

нечетно,

N

 

N

= −(N 1) 2

 

 

 

 

 

 

 

 

N + = N 2 1,

 

 

 

N = −N 2

N четно.

 

 

 

 

 

 

 

 

 

Для того чтобы совокупность лучей

{(n l, mΔθ): N n N +, 1 m M}

(3.127)

(3.128)

полностью покрыла круг радиусом T (т. е. область Ω), положим

l =

T

 

π

 

 

, Δθ =

 

.

(3.129)

N +

М

Р(n l, mΔθ) называется дискретными проекциями параллель-

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

 

 

 

 

Обозначим через θm величину mΔθ, а через (ν, θm ) одно-

мерное преобразование Фурье по l

проекции при угле θm

 

 

Т

 

 

 

 

(ν,θm ) = Р(l,θ)ei2πνldl ,

(3.130)

Т

где ν =12 l – пространственная частота.

В соответствии с дискретным преобразованием Фурье имеем

283

 

Σ (ν,θm )= (ν + γ l,θm ),

(3.131)

γ=−∞

где Σ (ν,θm ) – результат суммирования совокупности функций, получаемых из функции (ν) путем сдвига ее аргумента.

Учитывая (3.124), что P(l, θm ) = 0 при l >T , т. е. функция P(l, θm ) не является функцией с ограниченной шириной спектра, следовательно выборочные значения P(n l, θm ) и Σ (ν,θm ) под-

вержены эффекту наложения, который был рассмотрен в гл. 2. Влияние этого эффекта можно уменьшить фильтрацией функции

P(l, θm ) , осуществив ее до дискретизации. В реальности усред-

няющее действие такого рода осуществляется за счет апертурных ограничений (конечной ширины окна единичного детектора) блока детекторов. Даже в том случае, когда наложение мало, из выборочного характера данных следует, что функцию (ν,θm ) можно оце-

нить лишь внутри конечной полосы спектра шириной ν <12 l .

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

Рассмотрим проекционные теоремы.

Обобщенная проекционная теорема. Если h(l) – произвольная функция одного переменного, для которой существуют нижеследующие интегралы, то при любыхуглахθ справедливо равенство

Т

Р(l,θ)h(l )dl = ∫∫μ(x, y)h(x cosθ + уsin θ)dx dy . (3.132)

Т Ω

Доказательство. Рассмотрим левую часть равенства (3.132). Подставим в нее выражение для P(l,θ) из равенства (3.120) и, воспользовавшись выражением (3.122), сделаем замену переменных, перейдя от координат (l, t ) к прямоугольным координатам (x, y) .

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

284

екцией при угле θ эквивалентна определенной операции над исходным объектом μ(x, y). Если операция на μ(x, y) обратима, как например, преобразование Фурье, то из этого следует способ нахождения μ(x, y) по заданной P(l,θ).

Теорема о проекционном слое или центральном сечении.

Пусть (ν,θ) одномерное преобразование Фурье функции P(l,θ) по первому аргументу, определяемое формулой (3.130). Тогда

(ν,θ) = μˆ (νcosθ, νsin θ),

(3.133)

где μˆ – одномерное преобразование фурье-функции μ(x, y) . Доказательство. В качестве функции h(l) в выражении (3.132)

возьмем функцию exp(i 2πνl ) . Результат следует из определений

одномерного преобразования Фурье по (3.130).

В соответствии с физической интерпретацией этот результат называется теоремой о центральном сечении. Действительно, пусть θ – заданный фиксированный угол. Тогда в теореме утверждается, что при изменении ν величина (ν,θ) равна значению μ при ра-

диусе ν и угле θ в фурье-пространстве. Перейдем к выводу формул реконструкции.

Взяв обратное преобразование Фурье равенства (3.133), получим

π ∞

μ(х, у) = ∫ ∫ (ν,θ)exp(i 2πν(x cosθ + y sin θ)) ν dνdθ . (3.134)

0 −∞

Подстановка выражения (3.130) для (ν,θ) в (3.134) дает фор-

мулу реконструкции функции μ(х,у) по заданным P(l,θ).

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

μ(х,у), чем ее аппроксимация μ(х, у) с ограниченной шириной спектра. Преобразование фурье-функции μ(х, у) можно записать в виде произведения фурье-функции μˆ на функцию окна W (ν) :

ˆ

(νcosθ, νsin θ) W (ν) .

(3.135)

μ(νcosθ, νsin θ) = μˆ

285

 

Переписав (3.134) для функции μ(х, у), получим

π 12 l

μ(х, у) = ∫ ∫

0 12

(ν,θ)W (ν)exp i2πν(xcosθ+ ysinθ)

 

ν

 

dνdθ.

 

 

 

 

 

 

 

 

l

Заменим теперь (ν,θ)

 

 

(3.136)

 

ее выражением из (3.130) и заменим

порядок интегрирования по l

 

и ν . В результате получим

 

π Т

 

 

 

 

 

μ(х, у) = ∫ ∫ P(l,θ) g (x cosθ+ y sin θ−l )dl dθ,

(3.137)

0 Т

 

 

 

 

 

где

 

 

 

 

 

1 2

l

 

 

 

g (l ) =

 

ν

 

W (ν) exp(i 2πνl )dν ,

(3.138)

 

 

 

 

 

1 2

l

 

 

 

хcosθ+ y sin θ = l ' – луч, проходящий через точку P (см. рис. 3.39).

Выражения (3.137) и (3.138) составляют основу метода реконструкции – метода свертки и обратного проектирования для параллельных лучей.

Поясним название этого метода. Представим (3.137) в виде следующей последовательности операций:

 

 

Т

 

 

 

 

р(l ',θ) = P(l,θ) g (l 'l) dl ,

(3.139)

 

 

Т

 

 

 

 

π

 

 

 

μ(х, у) = p(x cosθ+ y sin θ, θ) dθ.

(3.140)

 

 

0

 

 

Промежуточная

величина p(l ',θ) ,

определяемая

формулой

(3.139),

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

l проекции

p(l,θ) при

угле θ и функции

g (l ) , заданной выражением (3.138). Функция

p(l ',θ)

называется свернутой проекцией при угле θ, а g (l ) – сво-

рачивающей функцией. Формула (3.140) отражает формальное применение прямого метода обратной проекции, описываемого

286

уравнением (3.126), существенной отличительной особенностью которого является операции над исходной проекцией p(l,θ) в виде

фильтрации сверткой (ибо свертка по физической сути является операцией фильтрации) с последующей фильтрацией в виде сворачивающей функции g (l ) со своим фильтрующим окном W (ν) .

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

gˆ (ν) =

 

ν

 

W (ν) ,

(3.141)

 

 

где gˆ (ν) – фурье-преобразование функции g (l ) .

Операция, описываемая выражением (3.140), называется обратным проецированием.

Таким образом, применение двойной операции фильтрации над исходной проекцией p(l,θ) в виде свертки, которая как бы усред-

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

μ(x, y) в виде μ(x, y) методом обратного проецирования.

Рассмотрим третий этап алгоритма – этап адаптации формулы обращения к дискретным зашумленным данным.

Для аппроксимации функции изображения μ(mx x, my y), где

тх, ту – масштабные коэффициенты дискретизации изображения; x , y – дискреты изображения по проекциям p(n l, θm ) , где

1 m M , N n N + (3.127)–(3.129), требуется алгоритм реконструкции, удобный для реализации в спецпроцессоре или универсальный ЭВМ.

Можно вычислить интеграл от обратной проекции (3.140) с помощью формулы трапеции:

287

μ(mx x, my

M

 

y) Δθp(mx x cosθm + my y sin θm , θm ). (3.142)

 

m=1

 

Для каждого значения θm необходимо найти значения сверну-

той проекции

p(l ',θm ) при Мх×Му значениях l, где Мх×Му – мат-

рица изображения исследуемого объекта.

 

Один из способов нахождения значений p(l ',θm )

состоит в вы-

числении свертки отдельно для каждого значения l '

при помощи

аппроксимации (3.139) по правилу трапеции для выборочных значений g (l 'n l ) сворачивающей функции. Такой подход слишком

трудоемок, поскольку величина Мх×Му обычно составляет от 5122 до 20482. Целесообразнее сначала оценить р(n l, θm ) при

N n N + , а затем посредством интерполяции по этим оценкам найти требуемые N значений функции р. Таким образом, аппрок-

симация свертки (3.139) осуществляется с помощью двух операций над дискретными данными: дискретной свертки, результат которой обозначим через рg , и последующей интерполяции рg , результат

которой обозначим через Ри . Эти операции можно записать в виде:

Рg (n' l, θm ) =

 

N +

 

l

P(n l, θm ) g ((n'n) l ) ,

(3.143)

 

n'=N

 

где N n' N + ;

 

 

 

Ри (l ',θm ) =

 

N +

 

l

Pg (n' l ) J (l 'n' l ) ,

(3.144)

n'=N

где J(l) – функция интерполяции, а число слагаемых (точек интерполяции) зависит от ширины интервала (детектора), J(l) отлична от нуля.

В качестве примера можно рассмотреть функцию Jл(l), соответствующую линейной интерполяции между соседними выборочными значениями:

288

1

(1

 

l

 

l ), l

l

 

 

 

 

Jл (l ) =

 

 

 

 

(3.145)

 

l

 

 

 

 

 

0,

 

 

 

 

l

l.

 

В этом случае сумма в (3.144) состоит всего из двух слагаемых. Необходимо заметить, что в дискретной свертке (3.143) используются выборочные значения g (n l ) , взятые через равные интер-

валы. Выбрав конкретную функцию окна, эти значения можно вычислить один раз и хранить в памяти ЭВМ.

В качестве функции окна можно взять

1− α

 

ν

 

ν

c

,

 

ν

 

 

≤ ν

c

,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

W (ν) =

 

 

 

 

 

 

 

 

 

 

 

 

 

 

(3.146)

 

 

 

 

 

ν

 

 

> νc ,

0,

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

где νс – частота среза, а коэффициент регуляризации α принимает значения от 0 до 1 (принципы регуляризации формул обращения рассматривались в главе 2). Положив νc =1 / (2 l) и подставив

(3.146) в (3.138), получим

g (n l ) =

3 2α

 

,

n = 0 ;

 

 

 

 

 

 

12( l )2

 

 

 

 

g (n

l ) = −

 

α

,

n – четно;

 

π2 (n l )2

 

g (n

l ) = −

 

1−α

,

n – нечетно.

(3.147)

π2 (n l )2

 

 

 

 

 

 

Последовательность значений g (n

l ) , соответствующая α = 0,

известна под названием дискретного сворачивающего ядра Лак- шминараянана–Рамачандрана [3].

После выбора g (n l ) и J (l ) с помощью алгоритма свертки и обратного проецирования вычисляются функции (3.143), (3.144) и (3.142), причем в (3.142) вместо Р используются значения Ри .

Рассмотрим соотношение между функцией Р(l ',θ) , определенной в (3.139), и функцией Ри (l ',θ), заданной выражениями (3.143) и (3.144). Это можно сделать, проанализировав одномерные преоб-

289

разования Фурье этих функций, которые обозначим соответственно через (ν,θ) и и (ν,θ) .

Из формулы (3.139) и теоремы о свертке для преобразований Фурье следует, что

(ν,θ) = (ν,θ) gˆ (ν) .

(3.148)

Преобразование Фурье выражений (3.143) (3.144) в результате дискретной свертки запишется

 

ˆ

(ν) .

(3.149)

и (ν,θ) = (ν,θ) g(ν) J

ˆ

 

 

 

Значение выражения (3.149) в том, что оно определяет факторы, влияющие на выбор интервала дискретизации l , функции окна W(ν), которая определяет сворачивающуюся функцию g(l), и интерполирующей функции J(l).

Поясним соотношение между этими параметрами на примере (рис. 3.42). В верхней части рисунка представлена зависимость

Σ(ν,θ) от ν при фиксированной θ. Заметим, что Σ(0,θ) является суперпозицией сдвинутых копий функции (ν,θ) .

Так как μ(х, у) = 0 вне ограниченной области Ω, ее фурье-образ не ограничен по протяженности. По этой причине фурье-преобра- зование (ν,θ) функции Р(l,θ) неограниченно и содержит другие компоненты функции Σ(ν,θ), а именно (ν + γ l ) при γ ≠ 0. Влияние эффекта наложения можно ослабить, уменьшив величину дискретизации (апертуру детектора) l или функцию Р(l,θ) до ее дискретизации, пропустив через фильтр нижних частот.

На рис. 3.42 показан график функции gˆΣ(ν)

при трех различ-

ных видах функции окна W(ν) в соответствии с выражением

(3.141). Кривая 1 соответствует функции окна

 

1,

 

ν

 

≤ ν

 

;

 

 

 

 

 

W (ν) =

 

 

 

 

c

 

(3.150)

 

 

 

 

 

0,

 

ν

 

> νc ,

 

 

 

 

 

 

 

 

 

где νc =1 2 l .

 

 

 

 

 

 

 

290

Соседние файлы в предмете [НЕСОРТИРОВАННОЕ]