Традиционные методы вычислительной томографии ФЕДЕРАЛЬНОЕ АГЕНСТВО ПО ОБРАЗОВАНИЮ Федеральное образовательное учреждение высшего профессионального образования «ЮЖНЫЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ» Д.Н. Карпинский МЕТОДИЧЕСКИЕ УКАЗАНИЯ к разделу «Традиционные методы вычислительной томографии» спецкурса «Применение томографических методов в медицинской диагностике» для студентов специальности «Прикладная математика» Ростов-на-Дону 2007 Печатается по решению кафедры теории упругости факультета математики, механики и компьютерных наук ЮФУ, протокол N1 от 10 сентября 2007 года. Методические указания разработаны доктором физико-математических наук, профессором кафедры теории упругости Д.Н.Карпинским. 1. ВВЕДЕНИЕ Томография - одно из бурно развивающихся направлений в области получения и обработки информации. Томография позволяет заглянуть внутрь наблюдаемого объекта. Основная проблема томографии - как по получаемым в томографическом эксперименте проекционным данным (например, по рентгеновским снимкам) "увидеть" внутреннюю структуру анализируемого объекта. Область математики, в которой разрабатываются методы решения подобных задач, известна как "интегральная геометрия" [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. ПРЕОБРАЗОВАНИЕ РАДОНА 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. МЕТОД РАЗЛОЖЕНИЯ В РЯД ФУРЬЕ (МЕТОД А. КОРМАКА) В этом разделе рассмотрим восстановление функции изображения по ее проекциям, полученным при помощи внешнего источника излучения. Запишем искомую функцию в полярной системе координат . Тогда по переменной , , произвольная двумерная функция будет периодической и ее можно разложить в ряд Фурье , . (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) 4. РЕКОНСТРУКЦИЯ ТОМОГРАФИЧЕСКИХ ИЗОБРАЖЕНИЙ ПРИ АППРОКСИМАЦИИ ПРОЕКЦИЙ ОРТОГОНАЛЬНЫМИ ПОЛИНОМАМИ 4.1. Рассмотрим алгоритм реконструкции изображения, основанный на приближенном представлении проекционных данных в виде конечного ряда ортогональных полиномов. Пусть имеется полная ортонормированная последовательность функций . Тогда, если искомая функция квадратично интегрируема, то она может быть представлена в виде , (4.1) где , (4.2) а - действительная неотрицательная весовая функция, относительно которой функции в области задания взаимно ортогональны. Учитывая равенство (5.1), задачу реконструкции функции по ее радоновскому образу можно сформулировать как задачу нахождения коэффициентов по получаемым проекционным данным. Формально это означает, что требуется найти соотношение, например, типа (4.2), но которое определялось бы не функцией , а .Вид искомого соотношения зависит от конкретной ортогональной последовательности и определить его в общем случае не удается. В [5] приводится решение данной задачи для ортогонального базиса, составленного из функций , (4.3) где - полиномы Цернике, для которых выполняютсясоотношения , . (4.4) Опуская громоздкие промежуточные выкладки, приведем окончательные выражения и сопроводим их необходимыми пояснениями, вскрывающими их физическую сущность. Предварительно заметим, что если изучаемая функция задана в некоторой ограниченной области , то всегда эту область можно охватить окружностью с некоторым минимальным радиусом а , и, положив в тех точках , ,где соответствующий круг не пересекается с , рассматривать задачу о восстановлении функции в пределах данной окружности. Далее, произведя нормировку координат , на величину , можно перейти к случаю восстановления функции в пределах окружности единичного радиуса. Лишь при выполнении данного условия возможно использовать последовательность функций (4.3). Для реконструкции функции заданной в круге единичного радиуса, нужно по полученным проекционным данным рассчитать величины , (4.5) где - полиномы Чебышева второго рода. Затем в равенство (4.1) вместо подставить найденные значения ,а в качестве использовать (4.3). При таких условиях последующее суммирование всех членов получившегося ряда позволяет реконструировать искомую функцию, так что , (4.6) где и - полярные координаты в плоскости , . Чтобы разобраться, почему суммирование в (4.6) по индексу проводится от до , достаточно вспомнить, что все коэффициенты при равны нулю. Выбор полинома Чебышева приводит к тому, что коэффициенты обладают еще одним свойством: они также равны нулю, когда сумма их индексов является нечетной. Это следует непосредственно из формулы (4.5), если учесть два обстоятельства: 1) согласно (2.8) ; 2) полином Чебышева четного (нечетного) порядка является соответственно четной (нечетной) функцией своего аргумента. Объединяя оба условия, имеем , если или нечетно. (4.7) Полезно также вспомнить, что для используемых полиномов Чебышева второго рода, которые определяются формулой , (4.8) (4.9) так что эти полиномы ортогональны на отрезке [- 1, 1] относительно весовой функции . Учитывая (4.9), можно показать [5], что . (4.10) Сопоставляя (4.6) и (4.10), видим, что, как искомая функция , так и ее радоновский образ , выражаются через двойные суммы по индексам и , в которых используются одни и те же коэффициенты , но разные последовательности ортогональных функций. Пример Пусть , ее радоновский образ находится по (2.7) при и оказывается равным . Согласно (4.5), если то (из-за центральной симметрии функции), а для получаем значения коэффициентов разложения = = Выполняя суммирование в (4.6) с данными коэффициентами получим приближенное значение исходной функции изображения . 5. РЕГУЛЯРИЗАЦИЯ ФОРМУЛ ОБРАЩЕНИЯ Обычно вместо точной проекции известна искаженная проекция , (5.1) где описывает соответствующую случайную погрешность, проявляющуюся в данном случае в виде аддитивной добавки. Тогда задачу реконструкции можно переформулировать следующим образом: требуется по приближенным проекционным данным найти приближенную функцию , которая в каком-то смысле хорошо описывала бы искомую функцию . Непосредственная подстановка "зашумленных" проекционных данных [7] в указанный вычислительный алгоритм приводит к большим искажениям в . Дело в том, что задача реконструкции относится к так называемым некорректным задачам [8]. Физическая суть "некорректности" проявляется в том, что если пользоваться точным решением некорректной задачи, то даже при небольших искажениях в исходных данных это решение может существенно отличатся от искомой функции . Устранить это нежелательное явление можно, регуляризируя формулы обращения. В методах, основанных на преобразовании Радона (раздел 2) для этого достаточно "подавить" влияние высоких частот в , что можно, например, достичь умножением на регуляризующие функции . Обычно регуляризующие функции выбирают в следующем виде: ; (5.2) ; (5.3) (5.4) Постоянная называется параметром регуляризации и подбирается эмпирически при расчете. Чем больше интенсивность ожидаемых искажений, тем больше должно быть значение параметра . Формулы обращения преобразования Радона (2.25) с учетом регуляризации получаются путем замены на , а (2.32) такой же заменой в (2.29). Что касается метода ортогональных полиномов (раздел 4), то описанный выше алгоритм реконструкции функции является точным в том смысле, что если ее радоновский образ известен точно, то по нему, в принципе, можно найти точные значения всех коэффициентов и далее по формуле (4.6) осуществить точное восстановление искомой функции. Однако на практике реализовать подобное точное восстановление невозможно. Этому препятствуют, по крайней мере, две причины. Первая кроется в самой сущности обсуждаемого алгоритма, ибо, для того чтобы он был точным, необходимо согласно (4.6) в общем случае определить бесконечное число членов . Вторая связана с невозможностью точного измерения радоновского образа. В результате определяемые по нему коэффициенты будут отличаться от их точных значений . Таким образом, в реальном алгоритме восстановления участвует ограниченное число членов ряда (4.6). Для определенности в дальнейшем будем считать, что ограничение проводится по индексу , так что . Этого условия достаточно, так как в силу (4.7) оно однозначно определяет конечное число всех отличных от нуля коэффициентов . Изменяя порядок суммирования в (4.6) и делая его аналогичным (4.10), имеем . (5.11) Известно [5], что ограничение суммирования в (5.1) приводит к функции , хотя и отличной от , но это отличие, оцениваемое по среднеквадратической погрешности , (5.12) будет минимально, если коэффициенты в (4.1) рассчитываются по прежним формулам (4.2). Данный факт говорит о том, что вынужденное на практике ограничение числа определяемых коэффициентов не должно привести к изменению тех формул, по которым они рассчитываются. С увеличением числа - членов суммы (5.11) погрешность (5.12) монотонно уменьшается. Важно подчеркнуть, что это происходит только тогда, когда коэффициенты известны точно. Если же они определяются с некоторыми ошибками, то отмеченная зависимость нарушается. В этом случае конкретный характер поведения погрешности (5.12) с ростом числа М во многом определяется статистикой ошибок измерения. В результате уменьшение усредненной погрешности за счет увеличения числа членов суммы(5.11) может происходить только до некоторого предела, после которого она начинает увеличиваться. Более того, часто при бесконечном увеличении числа слагаемых погрешность стремится к бесконечности. Таким образом, вторая причина, связанная с неточностью определения коэффициентов , так же, как и первая, приводит к необходимости использовать при восстановлении ограниченное число членов ряда (5.1), но в отличие от первой она указывает на то, что это ограничение возможно осуществить оптимальным образом. В данном случае не требуется регуляризации в том виде, в каком она была введена ранее. Ее роль как «сознательного ограничителя точности в идеальных условиях» будет выполнять «сознательное», оптимальное ограничение числа членов аппроксимирующих полиномов для данного уровня шумовых флуктуаций. ЛИТЕРАТУРА 1. Гельфанд, И.М. Интегральная геометрия и связанные с ней вопросы теории представлений [Текст]: монография / И.М Гельфанд, М.И. Граев, Н.Я. Виленкин. - М.: Физматгиз, 1962. - 656 с. 2. Календер, В. Компьютерная томография. Основы, техника, качество изображения и области клинического использования [Текст]: монография / В. Календер. - М.: Техносфера, 2006, -344 с. 3. Терещенко С.А. Методы вычислительной томографии [Текст]: монография /С.А.Терещенко. – М.: Физматлит, 2004. - 319 с. 4. Наттерер Ф. Математические аспекты компьютерной томографии: Пер. с англ. [Текст]: монография /Ф. Наттерер. -М.: Мир, 1990.-288 с. 5. Хелгасон, С. Преобразование Радона: Пер. с англ. [Текст]: монография / С. Хелгасон. - М.: Мир, 1983. - 152 с. 6. Хермен, Г. Восстановление изображений по проекциям: основы реконструктивной томографии: Пер. с англ. [Текст]: монография / Г. Хермен. - М.: Мир, 1983. - 349 с. 7. Троицкий, И.Н. Статистическая теория томографии [Текст]: монография / И.Н.Троицкий. – М.: Радио и связь, 1989. - 240 с. 8. Тихонов, А.Н. Методы решения некорректных задач. [Текст]: монография / А.Н. Тихонов, В.Я. Арсенин. - М.:Наука, 1986. - 287 с. 9. Гельфанд, И.М. Обобщенные функции и действия над ними [Текст]: монография / И.М. Гельфанд, Г.Е. Шилов. - М.: Физматгиз, 1959. - 497 с. ПРИЛОЖЕНИЕ А Чтобы вычислить (2.9), воспользуемся соотношением [9] , (A1) где - простые корни уравнения , - их число. Пусть . Тогда , , , так что . Подстановка (А1) в (2.9) дает = = , (А2) где при переходе к последнему равенству было учтено, что . ПРИЛОЖЕНИЕ Б Убедиться в справедливости (2.24) можно, если воспользоваться (2.8) и под интеграл в (2.17) вместо подставить , затем сделать замену переменных