МЕТОД РАЗЛОЖЕНИЯ В РЯД ФУРЬЕ (МЕТОД А. КОРМАКА)

ВВЕДЕНИЕ

 

Томография - одно из бурно развивающихся направлений в области получения и обработки информации. Томография позволяет заглянуть внутрь наблюдаемого объекта. Основная проблема томографии - как по получаемым в томографическом эксперименте проекционным данным (например, по рентгеновским снимкам) "увидеть" внутреннюю структуру анализируемого объекта. Область математики, в которой разрабатываются методы решения подобных задач, известна как "интегральная геометрия" [1].

Хронология развития вычислительной томографии:

1895 г. – открытие рентгеновских лучей;

1917 г. – преобразование Радона;

1920 г. – рентгенограмма в медицине;

1930 г. – линейная томография, вращательная томография;

1942 г. – РВТ в радиоастрономии;

1961 г. – сверточный алгоритм;

1964 г. – алгоритм РВТ А. Кормака;

1972 г. – серийный томограф Г. Хаунсфилда;

1977 г. – учебный курс по вычислительной томографии в университете штата Нью-Йорк;

1979 г. – Нобелевская премия А. Кормаку и Г. Хаунсфилду.

1.2 В настоящее время существуют следующие виды томографии:

1) рентгеновская томография;

2) радионуклеидная томография;

3) ЯМР – томография;

4) ультразвуковая томография;

5) оптическая томография;

6) протонно-ионная томография;

7) томография в радиодиапазоне;

8) ЭПР - томография.

Особенно важное значение методы томографии имеют для медицинской диагностики [2].

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

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

Рисунок 1. К выводу уравнения переноса излучения (1.1).

Стационарное уравнение переноса излучения в чисто поглощающей неоднородной среде, описывающее процесс излучения в веществе, представляет собой баланс частиц или энергии и имеет вид

 

(1.1)

 

Решением уравнения (2.1) будет закон Бугера-Ламберта-Бэра для неоднородной поглощающей среды, который составляет основу расчетов ТВТ.

 

, (1.2)

 

где - интенсивность источника излучения.

Рассмотрим теперь закон распространения излучения при действии внутренних источников излучения (самоизлучающие объекты).

Рисунок 2. К выводу закона переноса излучения при действии внутреннего источника.

 

Пусть точечный источник излучает в телесный угол с интенсивностью в веществе с распределением линейного коэффициента ослабления вдоль прямой, соединяющей источник с небольшой площадкой , наклоненной под углом к этой прямой. Тогда для интенсивности , приходящейся на площадку , получаем [3]

 

. (1.3)

 

Выражение (1.3) учитывает четыре основных фактора: пространственное распределение источника излучения, геометрическое ослабление, ослабление излучения в веществе и наклон площадки детектора. Формула (1.3) лежит в основе ЭВТ.

 

ПРЕОБРАЗОВАНИЕ РАДОНА

 

2.1 Рассмотрим задачу восстановления двумерного распределения коэффициента ослабления при просвечивании объекта излучением внешнего источника. Источник излучения проходит дискретно вдоль объекта. Синхронно с источником с другой стороны объекта движется детектор излучения. Набор отсчетов, полученный таким образом, определяет одномерную функцию, называемую проекцией. Затем система «Источник-детектор» поворачивается относительно объекта на некоторый угол , и снимает новый набор отсчетов, определяющий следующую проекцию. По полученному набору одномерных проекций необходимо восстановить двумерное распределение . Такую схему измерений называют круговой геометрией измерений, а проекции называют параллельными проекциями.

Рисунок 3. Схема кругового сканирования с параллельными проекциями.

 

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

Рисунок 4. К выводу формул преобразования Радона.

 

Известно, что всякая прямая может быть описана уравнением

 

, (2.1)

 

где - расстояние от начала координат до этой прямой; - угол, образованный с осью перпендикуляром, опущенным из начала координат на эту прямую.

Произвольная прямая однозначно задается двумя параметрами и . Поэтому и результат интегрирования функции по некоторой прямой будет зависеть от этих же параметров, т.е. . Предположим, что функция интегрируется по всевозможным прямым. Подобное интегрирование можно также рассматривать как некоторое преобразование, которое данной функции на плоскости ставит в соответствие функцию на множестве всех прямых, задаваемую интегралами от вдоль прямых. Это преобразование называют преобразованием Радона [4,5], а функцию часто называют образом функции в пространстве Радона или проекцией, которая в обозначениях (1.2) имеет вид

 

. (2.2)

 

Задача ставится следующим образом: функция неизвестна, но известна функция , являющаяся образом в пространстве Радона; требуется по функции определить . Другими словами решение поставленной задачи сводится к отысканию явной формулы обращения или к поиску преобразования, обратного преобразованию Радона. Впервые формула обращения была получена в статье Иоганна Радона, опубликованной в 1917 году в Трудах Саксонской академии наук. Однако эта работа была незаслуженно забыта и формула обращения была открыта заново в 1961 году.

Согласно определению радоновского образа и с учетом того, что интеграл от заданной функции вдоль прямой равен интегралу по всей плоскости произведения этой функции на - функцию, аргументом которой является левая часть уравнения (2.3), имеем [6,7]

 

. (2.3)

 

Интегрирование, осуществляемое по двум переменным, можно свести к интегрированию по одной переменной. Для этого введем еще одну прямоугольную систему координат , повернутую относительно на угол . Вспомним, что при переходе от одной из этих систем координат к другой координаты меняются следующим образом:

(2.4)

(2.5)

 

Сделаем в (2.3) замену переменных (2.4)

 

=

= (2.6)

 

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

 

. (2.7)

 

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

 

. (2.8)

 

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

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

Пример 1.

Пусть . Подставим это выражение в (2.6) и получим (см. Приложение А)

 

=

= . (2.9)

 

Из (2.9) следует, что если функция отлична от нуля в точке , то функция, описывающая ее образ в пространстве Радона , отлична от нуля на линии

 

, (2.10)

где .

Рисунок 5. - функция (а) и ее радоновский образ (б)

Пример 2.

Пусть . Подставляя это выражение в (2.6), получим

 

. (2.11)

Рисунок 6. Функция (а) и ее радоновский образ (б)

 

Область, где принимает максимальные значения, представляет собой линию, которая определяется выражением (2.10).

Пример 3.

При (2.12)

получаем

 

(2.13)

 

 

Рисунок 7. Функция (а) и ее радоновский образ (б)

 

2.2 В случае самоизлучающего объекта основной задачей ЭВТ является задача восстановления двумерного распределения источников излучения . Для простоты будем считать, что область, в которой распределены источники излучения, целиком расположена в области поглощения излучения, характеризующейся функцией распределения коэффициента ослабления . Обычно при измерениях с помощью ЭВТ, также как и при ТВТ, используют круговую схему с параллельными проекциями.

 

Рисунок 8. Круговая геометрия измерений в ЭВТ.

 

В [3] показано, что для ЭВТ с постоянным коэффициентом ослабления экспоненциальное преобразование Радона в декартовых координатах имеет вид

 

, (2.14)

 

а в полярных координатах

 

. (2.15)

 

Выражение (2.15) можно переписать в другом виде

 

. (2.16)

 

2.3 Выражения (2.3), (2.6) позволяет по функции найти ее радоновский образ . Существует соотношение, определяющее аналогичную связь между преобразованием Фурье этих функций. Пусть - одномерное преобразование Фурье функции по переменной , а - двумерное преобразование Фурье функции по переменным . Согласно определению

 

, (2.17)

. (2.18)

В трехмерном пространстве введем прямоугольную систему координат, у которой по одной оси отложены значения , а по двум другим – значения и .

 

Рисунок 9. Центральное сечение двумерного преобразования Фурье

 

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

Теорема.

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

 

. (2.19)

 

Для доказательства (2.19) подставим в (2.17) вместо выражение (2.6) и сделаем замену переменных, аналогичную (2.4), полагая в (2.4) . Тогда получаем

 

=

. (2.20)

 

Сравнивая последний интеграл в (2.20) с (2.18), видим, что они равны, если в (2.20) под и понимать соответственно и . Следовательно, последний интеграл в (2.20) равен , что и доказывает сформулированную теорему. Легко убедиться, что теорема о центральном сечении справедлива и для случая, когда верно равенство (2.7).

 

2.4 Рассмотрим теперь формулы обращения и алгоритмы реконструкции, основанные на теореме о центральном сечении. Известно, что по двумерному преобразованию Фурье можно найти :

. (2.21)

 

Сделаем в (2.21) замену переменных, перейдя в плоскости к полярным координатам , так что , . Тогда (2.21) принимает вид:

 

. (2.22)

 

Теперь воспользуемся равенством (2.19) и вместо подставим в (2.22) функцию , после чего получим

 

(2.23)

 

Равенство (2.23) и является искомой формулой обращения, позволяющей с учетом (2.17) по найти функцию . Однако привлечение этого равенства для обработки данных томографических экспериментов оказывается не очень удобным из-за используемой в нем области интегрирования. Принимая во внимание равенство

 

, (2.24)

 

получим окончательное выражение для обращения преобразования Радона (см. Приложение Б)

. (2.25)

 

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

(2.25) в несколько ином виде. Обозначим

 

(2.26)

 

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

 

. (2.27)

 

Тогда (2.25) принимает вид

 

, (2.28)

 

где при вычислении интеграла по величина должна быть заменена в соответствии с (2.26) на . В целом, алгоритм обращения преобразования Радона можно интерпретировать как последовательность операций:

1) для данного радоновского образа определяется его преобразование Фурье ;

2) функция умножается на ;

3) вычисляется обратное преобразование Фурье произведения и тем самым определяется функция ;

4) аргументу функции присваивается значение (2.26);

5) проводится интегрирование функции по углу .

Рассмотрим теперь иной вид формулы обращения по сравнению с (2.25). Обозначим через импульсную реакцию фильтра с частотной характеристикой . Связь между этими функциями устанавливается прямым и обратным преобразованием Фурье

 

(2.29)

(2.30)

 

Заметим, что функция обладает свойством .

Подставим в (2.25) вместо правую часть (2.30), а вместо - (2.17). Тогда получим

 

(2.31)

 

Интегрирование по дает , а последующее интегрирование по приводит к выражению

(2.32)

 

Выражение (2.32) отличается от (2.25) тем, что в последнем участвует преобразование Фурье радоновского образа, а в (2.32) сам радоновский образ. Алгоритм (2.32) можно представить как совокупность трех последовательных операций:

1) вычисляется свертка данного радоновского образа с функцией ;

2) аргументу функции , описывающей получаемую свертку, присваивается значение (2.26);

3) проводится интегрирование функции по углу .

2.5 Обращение экспоненциального преобразования Радона (2.14) – (2.16) представляет существенно более сложную задачу. Ограничимся здесь рассмотрением только случая радиально-симметричной функции . Тогда экспоненциальное преобразование Радона превращается в экспоненциальное преобразование Абеля [2]

 

= = .

 

В [2] показано, что обратное экспоненциальное преобразование Абеля имеет вид

 

=

. (2.33)

 

МЕТОД РАЗЛОЖЕНИЯ В РЯД ФУРЬЕ (МЕТОД А. КОРМАКА)

 

В этом разделе рассмотрим восстановление функции изображения по ее проекциям, полученным при помощи внешнего источника излучения. Запишем искомую функцию в полярной системе координат . Тогда по переменной , , произвольная двумерная функция будет периодической и ее можно разложить в ряд Фурье

 

, . (3.1)

 

Аналогично разложим в ряд Фурье по переменной проекцию

 

, . (3.2)

 

В полярной системе координат (2.3) имеет вид

 

, (3.3)

 

Далее найдем гармонику

 

=

= , (3.4)

 

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

 

,

 

где - функция Хевисайда, , . Следовательно,

 

, =

= (3.5)

 

где - многочлен Чебышева 1-го рода порядка . Выражение (3.5) представляет собой интегральное уравнение относительно неизвестной функции . В [3] показано, что решение (3.5) имеет вид:

. (3.6)

 

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

Для радиально-симметричной функции в полярной системе координат преобразование Радона превращается в частный случай преобразования Абеля

 

= =

= . (3.7)

 

В [3] показано, что решение интегрального уравнения (3.7) имеет вид

 

. (3.8)