Лабораторная работа. Изучение методов компьютерной томографии

Изучение методов компьютерной томографии.

1.2 Общие принципы томографии1.2 Общие принципы томографии

1.3 Классификация методов вычислительной томографииКлассификация методов вычислительной томографии

Особенности томографии плазменных объектов

Глава 2. Постановка задачи и методы томографирования.

2.1 Основные определения и постановка задачи томографии

2.2 Восстановление сечений с использованием Фурье преобразований

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

2.4 Обратная проекция с фильтрацией свёрткойОбратная проекция с фильтрацией свёрткой

2.5 Итерационные методы восстановления

2.5.1 ОпределениеОпределение

2.5.5 Алгебраический метод восстановления (ART) или получевая коррекция

2.5.6 Итерационный метод наименьшик квадратов (ILST) или одновременная коррекция
Глава 3. Описание программы TOMOGRAPHY1.

3.1 ВведениеВведение

3.2 Описание программы

Описание программы

Практические задания к лабораторной работе

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

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

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

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

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

Метод заметного повышения информативности рентгенограммы предложил еще в 20-х годах французский врач Бокаж, а его идею вскоре реализовал на практике итальянский инженер Валлебона, назвавший созданный им аппарат томографом. Суть идеи непосредственно видна из рис. 1.1. Если во время съемки рентгенограммы перемещать в параллельных плоскостях любые два из трех участвующих в ней компонента (рентгеновская трубка, объект, фотопленка), оставляя неподвижным третий, то на пленке даст четкое изображение только один слой, а элементы остальных слоев размажутся. Регулируя скорости перемещений, можно выделять различные слои и тем самым последовательно изучать структуру трехмерного объекта.

Рис. 1.1. Схема томографии по Бокажу.

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

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

Рис. 1.2. Схема вычислительной томографии.

Далее, перемещая луч по определенному закону в выбранной плоскости (например, параллельно самому себе или веером с некоторым углом раствора), получаем набор луч-сумм, определяющий одномерную проекцию. Будем для определенности говорить о параллельном смещении зондирующих лучей; тогда одной проекции соответствует некоторый угол относительно выбранной системы отсчета. Поворачивая затем луч в той же плоскости на малый угол, повторяем процесс регистрации луч-сумм и получаем новую проекцию и т. д. Процесс сканирования по углу продолжается до тех пор, пока полный угол поворота не составит 180°. В результате в памяти ЭВМ накапливается необходимая исходная информация для реконструкции изображения в выделенной плоскости (сечении). При этом отсутствуют какие-либо помехи в виде размытых элементов изображения, соответствующих другим плоскостям, и теоретически можно ожидать весьма высокого контраста томограммы. Возможность количественной оценки плотностей структурных элементов при использовании ЭВМ очевидна. Нетрудно убедиться, что метод пригоден также при исследовании излучающих и рассеивающих объектов, хорошо сочетается с техникой голографической интерферометрии, ЯМР и ЭПР, поляризационных измерений и т. п.

Описанный принцип томографирования в начале 60-х годов одним из первых применил А. Кормак в модельных экспериментах, ориентированных на медицинскую рентгенодиагностику.

Исключительная гибкость метода вычислительной томографии порождает, с одной стороны, большое разнообразие применяемых при получении томограмм физических, технических и расчётно-математических решений, с другой стороны, обеспечивает методу чрезвычайно широкое поле приложений.

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

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

Итерационные методы восстановления изображения используют аппроксимацию восстанавливаемого объекта массивом ячеек равной плотности, представляющих собой неизвестные величины, связанные системой линейных алгебраических уравнений, свободными членами которых являются отсчеты на проекции. Решаются системы уравнений итерационными методами, что и дало название данному классу методов восстановления. В настоящее время известно несколько итерационных методов восстановления изображения. Отличаются они в основном последовательностью внесения поправок во время итерации. Среди них наиболее известны и употребительны три метода: алгебраический метод восстановления (ART), метод одновременного итерационного восстановления (SIRT) и итерационный метод наименьших квадратов (ILST).

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

Глава 2. Постановка задачи и методы томографирования.
2.1 Основные определения и постановка задачи томографии.
Обычно для восстановления трехмерной функции используют одномерные проекции ее двухмерных сечений при фиксированном , т.е. как бы восстанавливают послойно при разных двухмерные сечения объекта ("томография"- описание сечений).

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

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

Рентгеновская трубка излучает пучок лучей в плоскости исследуемого сечения объекта. Рассмотрим случай, когда рентгеновский пучок имеет вид параллельных лучей. После прохождения рентгеновского пучка через сечение ослабленные пропорционально плотности объекта лучи воспринимаются детекторами. Система трубка - детекторы вращается вокруг исследуемого сечения. С системой трубка - детектор связана подвижная система координат , центр которой совпадает с центром системы , а оси повернуты относительно неподвижной системы на угол . Направление лучей совпадает с координатой u (рис.2.1).

Рис. 2.1. Расположение объекта в системе трубка - детектор в случае параллельных лучей: 1 - объект; 2 - источник излучения; 3 - приемник.

	(2.1)

и наоборот

	(2.2)

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

(2.3)

или для подвижной системы координат при различных получим набор проекций:

(2.4)

где - угол, под которым получена проекция. Информация о проекциях снимается с детекторов.

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

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

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

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

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

,	(2.5)

и наоборот

,	(2.6)

	(2.7)
	(2.8)

(2.8)- преобразование Фурье; (2.7) - обратное преобразование Фурье.

Положим в уравнении (2.8) , тогда

(2.9)

Сравнив (2.9) с (2.4) видим, что Фурье-образ искомой функции в точках Фурье-плоскости, лежащих на прямой , есть одномерный Фурье-образ соответствующей проекции - это так называемая теорема о проекциях и сечениях (2.11):

т.е. так как

, то

	(2.10)
	(2.11)

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

элемент площади при переходе от декартовых координат к полярным

,тогда

	(2.12)

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

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

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

Рис. 2.3a Проекции снятые с объекта (объект - белое пятно)

view1,view2, view3 - проекции.

Рис. 2.3b Восстановление объекта (белого пятна) методом "обратной проекции"

view1,view2, view3 - проекции.
Восстановление производится путём обратного проектирования каждой проекции через плоскость, т.е. величина сигнала, соответствующая данной лучевой сумме, прикладывается ко всем точкам, которые образуют этот луч. После того, как это сделано для всех проекций, получается приближённая аппроксимация исходного объекта. Так как в этом методе проекции как бы растягиваются обратно через восстанавливаемое сечение, метод назван обратной проекцией. Для каждой точки изображения восстановленная плотность является суммой всех лучевых проекций, которые проходят через эту точку. Поэтому метод обратной проекции иногда называется методом суммирования или линейной суперпозиции. Математическое описание метода обратной проекции может быть представлено:

если в (2.12) принять , тогда

(2.13)

Или при измерении конечного числа проекций

(2.14)

где суммирование производится по всем углам проекции . Аргумент

(2.15)

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

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

Рис. 2.4 Восстановление точечного объекта методом "обратной проекции": а - объект (1) и его проекции (2); б - обратное проецирование; ясно видно формирование фона в пространстве, окружающем точку.

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

Предположим, что есть Фурье-образ некоторой функции ,

(2.16)

тогда произведение на Фурье-образ проекции можно рассматривать как Фурье-образ свертки самой проекции и

(2.17)

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

(2.18)

где - результат свертки измерений проекции с некоторой функцией , называемой ядром свертки. Как и в методе обратного проецированияМетод обратной проекции, суммирование идёт по М - количеству проекций, а коэффициент представляет угловое расстояние между соседними проекциями.

Точное вычисление невозможно ввиду расходимости этого интеграла, но это и не нужно, т.к. дискретизация отсчетов на проекциях через интервалы автоматически ограничивает область задания максимальным значением . Ядро в пределах от до вычисляется:

(2.19)

Если положить , то .

Рис. 2.5. Функция ядра свёртки (по обеим осям относительные единицы).

Процедура восстановления функции по методу обратной проекции с фильтрацией сверткой будет следующая:

1. Находим функцию как интеграл свертки проекции с функцией ядра свертки .

2. Определяем функцию в точках v.

3. Вычисляем искомую функцию , используя уравнение (2.18).
Удобство алгоритма обратной проекции с фильтрацией состоит в том, что вычислительный процесс восстановления может идти почти одновременно с регистрацией проекции: как только получена очередная проекция , осуществляется ее свертка с ядром и добавление данных в ячейки памяти, накапливающие результаты суммирования по формуле (2.18). Как только обработана последняя проекция, функция восстановлена.
2.5 Итерационные методы восстановления
2.5.1 Определение.
Термин “итерационный” относится к методу последовательных приближений, при котором выбирается произвольное начальное изображение; для него рассчитываются проекции, а затем в изображение вводятся поправки для лучшего согласования этих проекций с измеренными проекциями. Итерации повторяются до тех пор, пока не будет получена удовлетворительная сходимость.
2.5.5 Алгебраический метод восстановления (ART) или получевая коррекция.
Этот метод был независимо разработан Гордоном и др., который назвал его методом алгебраического восстановления (ART – Algebraic Reconstruction Technigues). Процедура восстановления следующая. В каждой итерации вычисляется сначала одна лучевая сумма при исходном значении плотности в ячейках . По ней и с учетом измеренной лучевой суммы определяется поправка, которая вводится во все точки, входящие в состав данного луча. Затем операция повторяется для второго луча, третьего и т.д. При этом поправки, введенные от предыдущей лучевой суммы, учитываются в каждом новом расчете. Эти операции продолжаются до тех пор, пока не будут обработаны все проекции, после чего итерация завершена. Если критерий получения решения не удовлетворён, то происходит переход к следующей итерации.
Рассмотрим элементарный пример реконструкции двумерного объекта, когда сам объект предельно прост, алгоритмический “механизм” совершенно прозрачен, а привлечения ЭВМ вообще не требуется.

Возмём в качестве такого объекта квадрат ABCD (рис. 2.7.), разделённый на 9 равных клеток (ячеек). Числа от 1 до 9, разбросанные по клеткам, соответствуют плотности или какой-нибудь другой характеристике, находимой томографически. Пусть известны 4 проекции, определяемые направлениями сторон квадрата AB и AD и его диагоналей AC и BD. Если в каждой проекции взять по 3 луч-суммы, то в первых двух случаях вклад внесут все клетки, а в двух других лишь по 7 клеток из 9. Таким образом, мы исходим из 12 значений луч-сумм и ищем 9 структурных элементов объекта, т.е. решаем переопределённую задачу.

Рис. 2.7. Тест-объект из девяти элементов.

(2.26)

где индекс ‘k’ нумерует клетки, i обозначает номер итерации, - исходные значения на рис. 2.7., то .

Перейдем к следующей проекции (лучи идут вдоль AD). Теперь для каждого луча следует скорректировать сумму чисел, получаемых после первой итерации, на известную луч-сумму данной проекции. Так, сумму в первом столбце (15) следует, очевидно, уменьшить на 2 и вычесть из каждого числа по 2/3; в третьем столбце нужно, наоборот, добавить в каждую клетку по 2/3; во втором столбце изменений нет. Видно, что в нашем примере вторая итерация фактически оказывается неинформативной: предыдущее однородное распределение лишь слегка деформируется, создавая небольшой градиент вдоль AB и совершенно не выявляя сложной структуры объекта (рис. 2.8 B). Погрешность восстановления, оцениваемая по (2.26), даже несколько возрастает по сравнению с первой итерацией:.

Рис. 2.8. Различные стадии (A-F) восстановления тест-объекта, изображенного на рис.2.7.

Далее можно вновь привлечь первую проекцию и начать, таким образом, второй цикл. На рис. 2.8 Е,F показаны результаты, получаемые после шестой () и восьмой () итерацией. Процесс можно было бы, конечно, продолжать и дальше, но уже и из проделанных выкладок ясно, что 7-8 итераций позволяют получить неплохой результат восстановления.