Методы восстановления изображений в рентгеновской томографии Текст научной статьи по специальности «Математика»

МЕТОДЫ ВОССТАНОВЛЕНИЯ ИЗОБРАЖЕНИЙ

В РЕНТГЕНОВСКОЙ ТОМОГРАФИИ И.П. Гуров, В.С. Сизиков, Д.С. Щекотин

Рассмотрено два типа задач реконструкции изображений в рентгеновской томографии: задачи с полными и неполными данными. Проанализированы аналитические методы (обращения Радона и Кормака) и численные методы (преобразования Фурье, свертки и обратной проекции, итераций и регуляризации) решения задачи реконструкции. К методу свертки и обратной проекции применен метод регуляризации Тихонова, а вместо итерационного методы Качмажа использован метод итераций Фридмана. Приведены результаты компьютерного моделирования.

Постановка задачи

Рассмотрим задачу восстановления рентгеновского изображения на примере рентгеновского томографа с параллельной схемой сканирования (см. рис. 1) [1-6].

Л ГУ

Рис. 1. Параллельная схема сканирования

На рис. 1 показано сечение Э исследуемого объекта, характеризуемое плотностью вещества, точнее, коэффициентом поглощения рентгеновских лучей /(х, у), где х, у - неподвижная относительно объекта система декартовых координат (с ней совмещена система координат х', у', необходимая далее). Расположенные на раме рентгеновские трубки излучают узконаправленные рентгеновские лучи интенсивности 10, которые, пройдя через вещество и испытав частичное поглощение, регистрируются соответствующими детекторами (приемниками). Такой эксперимент проводится для ряда значений угла поворота рамы ф е [0, п).

Введем вращающуюся (неподвижную относительно рамы) систему декартовых координат s, l. Тогда можно обозначить через I (s, ф) значения интенсивности излучения, принятые детекторами. Обозначим далее

q(s,<p) ^In1^. (1)

1 о

Отношение 1 (sф) принято называть прозрачностью, а функцию q(s, ф) -1 о

поглощением [6, с. 19], или тенью. Искомая функция f (х,y) связана с измеренными значениями q(s, ф) интегральным уравнением Радона [1-4, 6], или теневым уравнением:

Rf - jf (х, y) dl = q(s, ф), (2)

L( «,ф)

где R - оператор преобразования Радона (оператор теневого преобразования), а интегрирование ведется по лучу в виде прямой L(s, ф), уравнение которой

х cos ф + y sin ф = s. (3)

Интеграл jf (х, y) dl называется массой вещества на луче зрения [6, с. 19], а

L( ^ф)

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

Тень изображения при фиксированном угле ф обозначим через qф (s) и будем называть линейной интегральной проекцией изображения q под углом ф.

Заметим, что существует еще схема веерного сканирования (более распространенная схема), для которой также справедливо уравнение типа (2) [4].

Виды задач реконструкции

Задача реконструкции рентгеновского изображения сводится к решению интегрального уравнения типа (2) относительно f (х, y) по известным значениям q(s, ф). При этом существует два типа задач реконструкции: задача с полными данными и задача с неполными данными [4].

В задаче с полными данными предполагается, что исходные интегральные проекции q(s, ф) известны вдоль всех лучей, проходящих через исследуемое сечение D, т.е. при всех необходимых значениях s и ф .

Задачи с неполными данными могут быть следующих типов:

1. задача с ограниченным диапазоном углов,

2. внешняя задача, когда функция q (s, ф) задана лишь для | s | > a, где a > 0,

3. внутренняя задача, когда q(s, ф) задана лишь для | s | < a, где a > 0,

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

Что касается аналитических и численных методов решения уравнения (2), то известны следующие аналитические методы: согласно формуле обращения Радона [2, 4] и согласно формуле обращения Кормака [4]. Также существуют следующие численные методы: метод преобразования Фурье, метод свертки и обратной проекции, метод итераций и методы регуляризации.

Многие из этих методов исследованы достаточно подробно [1-6], однако остался ряд нерешенных вопросов. Например, в методе свертки и обратной проекции [4, 5] имеет место явная расходимость решения f (х, y) и, чтобы ее устранить, используется

предельная (максимальная) частота Фурье, хотя, как известно [7, с. 259-261], более эффективно устраняют расходимость методы сглаживающих окон и, тем более, метод регуляризации Тихонова. Это показано в работе [3, с. 33], однако введенный в ней стабилизирующий множитель можно значительно упростить, что и предлагается в данной работе.

Кроме этого, в качестве метода итераций в [4] использован метод Качмажа. Между тем, более эффективным представляется метод итеративной регуляризации Фридмана [7, с. 272].

Наконец, представляются совершенно новые возможности в решении задачи реконструкции рентгеновских изображений в связи с приведением «неудобного» уравнения (2) к стандартному уравнению - двухмерному интегральному уравнению Фредгольма I рода типа свертки [8]:

ÍÍ

/ (хy') dx 'dy' = ^ (x, y), (4)

V(x - x ')2 + (y - y f

где

n

S(x, y) = — I q(x cos ф + y sin ф, ф) dф. (5)

n J

0

Уравнение (4) уже исследовалось [2] и применялось на практике [1], однако до настоящего времени не рассмотрен вопрос о его решении в случае задачи с неполными данными.

Метод преобразования Фурье

Этот хорошо известный метод основан на непосредственной численной реализации теоремы Брейсуэлла о центральном или проекционном слое (проекционной теоремы) [9]:

q (ш, ф) = f (ш cos ф, ш sin ф), (6)

где

ТО

q(ш, ф) =| q(s, ф) e~iasds, (7)

— ТО

ТО ТО

М, ш2) = <x, y) (8)

-ТО-ТО

- преобразования Фурье (ПФ). Искомая функция в форме обратного ПФ равна

ТОТО

f (x, y) = 4Л7 ÍÍ f (ш 1, ш 2) ei(ШlX+Ш2 y} dШl dш 2. (9)

-ТО-ТО

Однако при численной (дискретной) реализации Фурье-алгоритм дает значительные искажения (артефакты) [4]. Это связано с тем, что в (9) требуется знание f (ш 1, ш 2), а из (6) мы получаем f (ш cos ф, ш sin ф). Поэтому при дискретизации задачи требуется интерполяция, что ведет к большему числу операций и потере точности [4].

Для устранения недостатков стандартного алгоритма ПФ разработаны усовершенствованные Фурье-алгоритмы [4]. Тем не менее, более популярным и эффективным является алгоритм свертки и обратной проекции.

Метод свертки и обратной проекции

—то—ТО

Данный метод, восходящий к работам [10, 11] и нашедший последующее развитие в работах [3-5] и др., является дальнейшей эволюцией Фурье-алгоритма и позволяет полностью избежать интерполяции за счет перехода от декартовых координат ю2 к полярным координатам ю, ф в пространстве частот. Искомое решение имеет вид [3, с. 32]:

П X

f (x, y) = —- I ёф I q(s, ф) p(x cos ф + y sin ф-s) ds , (10)

4п 2 J J

0 -x

где

X

p(t) =| | ю | в'ю'ёю (11)

— X

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

Однако решение (10)—(11) дает расходимость. Действительно, используя формулу Эйлера и записав (11) в виде

X

p(t) = 2| | ю |cos(ot)ёю , (12)

0

видим, что p(t) расходится при любом t е (—x, x). Чтобы устранить расходимость,

в [4, 5] вводится предельная частота ю max = —, где h — шаг дискретизации, в результате

h

ю 2 (ю t ^

^ mav • 9 ^ max

p(t) = юmax sinc(юmax t) -^f^inc2 ^ . (13)

2 V 2 У

Такой устойчивый алгоритм часто называют решением в пространстве S—¡ h [7, с.

259], но он является самым грубым из устойчивых алгоритмов. Более эффективное подавление высоких гармоник Фурье обеспечивают методы сглаживающих окон [7, с. 259—260]. А наиболее эффективным является метод регуляризации Тихонова [7, с. 260— 261].

Метод регуляризации Тихонова

Применение данного метода в работе [3, с. 33] позволило получить устойчивый (регуляризованный) вариант решения (10)—(11) в форме:

П X

fa (Х У) = 1 ёф I q(S ф) pа (Х c0s ф + У sin ф - s) ds , (14)

0 —X

где

X

pa (t) = ||ю|Га (|ю|) еШёю , (15)

—X

а > 0 — параметр регуляризации,

Wa(|»i) = H 2(Я2,('ю|) 2p (16)

H2(| ю |) + аю2p

— стабилизирующий множитель, причем

H(| ю |) = —, (17)

— ю

р = 1,2,3,... - порядок регуляризации. Обычно р = 1, а а выбирают способом невязки или способом подбора [6, с. 194], [7, с. 242-249].

Однако выражение (16) является довольно громоздким. Более простым, но столь же эффективным, как следует из работы [12] (см. также [6, с. 170]), является определение стабилизирующего множителя в виде

^а(|ш|) = 1 2р . (18)

1 + аш 2 р

Метод итераций Фридмана

Основное достоинство методов итераций заключается в том, что они пригодны для восстановления томографических изображений в случае задачи с неполными данными. В качестве примера рассмотрим метод итераций Фридмана [7, с. 272]. Метод заключается в том, что вначале выполняется инициализация:

~ 2

/0 еП, 0 <р<--,

' р И2>

а затем выполняются итерации при г = 1, 2, ...:

/ = И// - И / ],

/ +1 = / + в 4 ,

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

информацию, П - пространство изображений, И - оператор Радона, а И * -сопряженный ему оператор.

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

сходимость гарантируется только при условии 0 <в<2/||И||2, а это требует знания

нормы оператора И . Эту норму обычно трудно или невозможно вычислить аналитически. Поэтому на практике берут достаточно малое значение в. Кроме того, определенную сложность представляет определение числа итераций (процесс итераций сходится лишь в случае точных значений q(s, ф), а при наличии шумов он расходится из-за некорректности задачи). Обычно число итераций выбирают по невязке [7, с. 273274] или по поправке [7, с. 274-275]. В данной работе число итераций п выбиралось способом подбора - путем выбора такого числа итераций, при котором изображение восстанавливается наилучшим образом с точки зрения зрительного восприятия. Такой способ, конечно, содержит элемент субъективности, однако он основывается на априорной информации о решении и, как показывает решение различных модельных примеров, является весьма эффективным.

Моделирование

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

На рис. 2 показаны примеры исходных изображений размером 512 х 512 пикселей. После решения прямой задачи томографии получены тени изображений, показанные на рис. 3. После восстановления изображений методом свертки и обратной проекции и методом итераций Фридмана получены восстановленные изображения (рис. 4 и 5).

Итак, в данной работе показана эффективность использования метода регуляризации Тихонова и метода итераций Фридмана для реконструкции рентгеновских изображений.

4 \у X

--►

Рис. 2. Исходные изображения.

Рис. 3. Тени изображений

41 у X у

-►

Рис. 4. Первое восстановленное изображение: методом свертки и обратной проекции (слева) и методом итераций Фридмана (справа).

Рис. 5. Второе восстановленное изображение: методом свертки и обратной проекции (слева) и методом итераций Фридмана (справа).

Литература

1. Тихонов А. Н., Арсенин В. Я., Рубашов И. Б., Тимонов А. А. Первый советский компьютерный томограф // Природа. 1984. № 4. С. 11-21.

2. Тихонов А. Н., Арсенин В. Я., Тимонов А. А. Математические задачи компьютерной томографии. М.: Наука, 1987.

3. Троицкий И. Н. Статистическая теория томографии. М.: Радио и связь, 1989.

4. Наттерер Ф. Математические аспекты компьютерной томографии. М.: Мир, 1990.

5. Суинделл Б., Уэбб С. Рентгеновская трансмиссионная компьютерная томография // Физика визуализации изображений в медицине. Под ред. С. Уэбба. М.: Мир, 1991. Т.1. С.138-173.

6. Сизиков В. С. Математические методы обработки результатов измерений. СПб: Политехника, 2001.

7. Верлань А. Ф., Сизиков В. С. Интегральные уравнения: методы, алгоритмы, программы. Киев: Наук. думка, 1986.

8. Тихонов А. Н., Арсенин В. Я., Рубашов Н. Б., Тимонов А. А. О постановке основных задач вычислительной томографии. М.: Препринт ИПМ АН СССР, 1982. № 141.

9. Bracewell R. N., Riddle A. C. Strip integration in radio astronomy // Austral. J. Phys., 1956. Vol. 9. P. 198-217.

10. Bracewell R. N., Riddle A. C. Inversion of fan-beam scans in radio astronomy // Astrophys. J. 1967. Vol. 150. P. 427-434.

11. Ramachandran G. N., Lakshminarayanan A. V. Three-dimensional reconstruction from radiographs and electron micrographs: application of convolutions instead of Fourier transforms // Proc. Nat. Acad. Sci. US. 1971. Vol. 68. P. 2236-2240.

12. Сизиков В. С. Использование регуляризации для устойчивого вычисления преобразования Фурье // Ж. вычисл. матем. и матем. физ. 1998. Т. 38, № 3. С. 376386.