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

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

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

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

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

Рассмотрим задачу восстановления рентгеновского изображения на примере рентгеновского томографа с параллельной схемой сканирования (см. рис. 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.