ФИЗИКО-МАТЕМАТИЧЕСКИЕ ОСНОВЫ КОМПЬЮТЕРНОЙ ТОМОГРАФИИ

 

Филонин О.В.

 

NB

  1. На данной странице приведены только основные расчетные соотношения, необходимые для выполнения расчетных работ.
  2. Приведенный здесь материал, не может заменить систематического курса лекций по КТ, созданного профессором, д.т.н. Филониным О.В., содержащего 26 лекций. Рабочая программа полного курса, приведена на соответствующей странице, данного ресурса.
  3. Полный курс лекций можно приобрести, связавшись с автором.

 

ЛЕКЦИЯ 1.

 

КОМПЬЮТЕРНАЯ ТОМОГРАФИЯ — СОВРЕМЕННЫЙ МЕТОД ДИАГНОСТИКИ

 

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

Компьютерная (вычислительная, реконструктивная) томография не имеет ничего общего с описанным выше томографическим методом диагностики. Методы компьютерной (вычислительной) томографии основаны на теореме Радона, доказанной в 1917 году. Долгое время казалось, что эта теорема имеет чисто академический интерес. Суть теоремы Радона заключается в том, что если имеется некоторая неизвестная многомерная функция , где  в общем виде многомерный параметр, и если известны интегралы от этой функции вдоль некоторых направлений

 

,

 

то имеется возможность вычислить с  помощью формулы (обращения) Радона (см. ниже) саму функцию . Сложность для реальных вычислений с помощью формулы Радона заключается в том, что как будет показано ниже, в знаменателе стоит преобразование Гильберта. Эта процедура приводит к “разрывам” в результирующих вычислениях, особенно при дискретных представлениях исходных интегральных функций . В настоящее время разработано достаточно много эффективных методов, основанных на определенных приближениях и позволяющих «обходить» указанные сложности. Но к сожалению эти методы требуют для своей реализации довольно мощные многопроцессорные специализированные вычислительные системы.

Поэтому и сегодня алгоритмы, основанные на «чистых»  прямых и обратных преобразования Радона используются достаточно редко, особенно в массовых системах диагностики. 

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

Указанные методы достаточно просто реализовывались в компьютерных системах и позволяли получать точность реконструкции 0,1% ¸1 %, что недоступно для любых других методов диагностики. Самым главным в компьютерной томографии является то обстоятельство, что с ее помощью удается поучать высококачественную информацию о внутреннем - пространственном строении объекта. Вначале такая информация представлялась в виде поперечных  срезов, но с развитием вычислительных средств сегодня, мы имеем возможность знать, распределение любого интересующего нас параметра ко всему объему исследуемого объекта.

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

  • начиная от тела человека и др. медико-биологических объектов,
  • изделий из композиционных материалов,
  • диагностика параметров ДВС, ГТД,
  • кончая такими экзотическими объектами, как плазма в установках ТОКАМАК и СТЕЛАРАТОР.

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

Компьютерная томография развивается по двум, относительно самостоятельным направлениям:

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

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

 

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

 

ФИЗИКО-МАТЕМАТИЧЕСКАЯ ПОСТАНОВКА ЗАДАЧИ КОМПЬЮТЕРНОЙ ТОМОГРАФИИ.

 

А) Формализованная математическая постановка задачи

 

 

Предположим, что имеется некоторая двумерная функция, заданная в неподвижной системе координат как f(x,y), (см. рис.1). Совместим с осью ОХ  полярную систему координат . Понятно, что искомую функцию можно будет определить и в полярной системе, как f(r). Задача, таким образом, сводится к тому, чтобы вычислить неизвестную пока функцию f(x,y) или, в более обобщенном представлении: . Для того, чтобы обсуждать в дальнейшем методы реконструкции , рассмотрим следующие геометрические построения:

  1. Совместим с началом координат системы XOY начало координат другой декартовой системы X`O`Y`, последняя система координат может вращаться в одной плоскости с O` и O. Координаты любой точки f(x,y) связаны с системой X`O`Y` соотношением:

 

 

 

           (1)      

          

 

Выражение (1) удобно записать, введя матрицу поворота в виде:

 

        (2)

 

 

 

  1. В системе координат (X`O`Y`), выберем любое направление (прямую), параллельное оси O`Y` и пересекающее искомую функцию f(x,y). Точки 1 и 2 являются граничными точками Î f(x,y). Вся совокупность точек внутри интервала (1,2) принадлежит искомой функции f(x,y), и определяет ее значение в этом интервале. Любую из точек интервала (1,2) можно также определить с помощью радиуса вектора , “конец” которого скользит вдоль отрезка (1-2).Таким образом, любая точка на отрезке (1-2) определяется параметрами . Выбранное направление  в декартовых системах К и К`, однозначно определено углами a, b и координатой xo'. В полярной системе координат выбранное направление можно определить с помощью нормального вектора , соединяющего точку О и отрезок (1-2). Понятно, что этот вектор «попадает» в точку xo', координаты которой в полярной системе определяются числами .
  2. Введем понятие “функция - проекции”:

 

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

  • вычислим интеграл вдоль выбранного направления 1-2, то есть:

                (3)

 

 

выражение (3) даёт интегральное значение – радоновский образ вдоль выбранного луча - хорды (1-2)  в координате xo',

  • рассмотрим бесконечное множество лучей, параллельных и распределенных вдоль оси O`X` в пределах . Тогда выражение (3) можно будет переписать в виде:

      (4)

 

 

Примечание: в выражении (4) произведена замена пределов  интегрирования, в предположении, что искомая функция f(x,y) (или ) гладкая, непрерывная и быстро убывает вне контура, ограничивающего область ее существования, изображенного на рис.1. (Другими словами, выбор граничного контура, определяющего конечную область носителя, практически полностью определяет искомую функцию).

Используя соотношение (2), выражение (4) можно переписать в виде:

                     (5)

 

В принципе выражение (5) и определяет функцию проекции, заметим, что a, b изменяются в пределах , следовательно, согласно (5) мы имеет дело с бесконечным набором одномерных функций проекций, каждая из которых является интегралом вдоль некоторого направления от искомой двумерной функции f(x,y).

NB Гораздо удобнее определить функцию проекции в полярной системе координат, причем целесообразно, исходя из соображений регистрации реальных сигналов выбрать другую полярную систему координат ), полярная ось которой ; но эта ось не должна принадлежать контуру определяющему , или , см. рис. 2.

В первом приближении функцию проекции можно определить:

 

               (6)

 

 

 

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

 

Поэтому целесообразно ввести в выражение (6) d-функцию и представить функцию проекции в одном из следующих видов:

       

 

          (7)

 

 

 

 

 

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

  • “Внешняя” система координат (X`O`Y`, рис.1 и ОР, см. рис. 2) должна “вращаться” в пределах или относительно системы XOY.
  • для любого угла поворота  см. рис. 2, необходимо вычислить одномерные интегралы вдоль направления 1-2, таким образом, получается набор одномерных функций проекций от искомой функции   или .

 

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

Для данного направления - l (см. рис. 2), линия - (1-2), в декартовой системе координат К` функцию проекции очевидно можно записать в виде:

   (8)

 

 

Заметим, что выражение (5) соответствует фиксированному углу a. Под радоновским образом функции  понимается непрерывная совокупность значений , при изменении параметра , т.е.:

 

 

  (9)

 

 

 

Выражения (8), (9) можно записать в полярной системе координат (см. рис. 3), следующим образом:

 

  (10)  

 

                                

    (11)    

 

 

Заметим, что выражение:

 

,

 

как известно, определяет взаимосвязь между полярной и декартовой системами координат и может быть выражено соотношениями:

 

,

 

 

т.е. описывает «вращение» системы k` относительно k, см. рис (1), (2).

 

ЛЕКЦИЯ 2.

 

ФИЗИЧЕСКИЕ ЯВЛЕНИЯ И ПРОЦЕССЫ, ИСПОЛЬЗУЕМЫЕ ПРИ ФОРМИРОВАНИИ ФУНКЦИЙ ПРОЕКЦИЙ В КОМПЬЮТЕРНОЙ ТОМОГРАФИИ

 

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

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

В то же время может использоваться и «внешнее излучение»: световое, рентгеновское, гамма - кванты, потоки нейтронов. Поток такого излучения, прошедшего через исследуемый объект в данном направлении может рассматриваться как функция проекции, при определенных условиях и способах регистрации.

Напомним, что оптическое излучение атомов формируется вследствие дискретных переходов электронов внешних оболочек, уровней (s, p) при наибольшем для данного атома значении n (главного квантового числа). Спектр такого излучения линейчатый, характерный только для данного атома. Другим видом атомного излучения является излучение, генерируемое вследствие  колебательного движения атомов, находящихся в молекулярных структурах, или в узлах кристаллической решетки. Спектры таких видов излучения являются сплошными, и «расположены» в инфракрасной или красной областях спектра.

Рентгеновское излучение формируется, либо вследствие дискретных переходов атомных электронов глубинных оболочек, при этом спектр такого излучения можно считать линейчатым, образующий K, L, M серии, либо вследствие торможения набегающего (внешнего электрона данным атомом) его электронной оболочкой или полем ядра. В последнем случае спектр излучения сплошной, в некоторых диапазонах энергий он может быть описан аналитически, например, формулой Шиффа.

Ядерное излучение можно разделить:

  • на электромагнитное, минимальная энергия кванта такого излучения 0,511 МэВ, и такие кванты получили название g - квантов,
  • электронное излучение,
  • потоки нейтронов, протонов, a - частиц, осколки ядер.

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

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

 

Формирование проекционных данных с помощью оптического излучения объектов

 

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

Наибольшее применение методы и средства оптической томографии нашли при исследовании параметров таких самосветящихся объектов, как: потоки высокотемпературных газов  (анализ процессов воспламенения и горения в камерах сгорания двигателей внутреннего сгорания, газотурбинных двигателей, ракетных двигателей и т.д.), пламена, низкотемпературная плазма, дуговые разряды, высокотемпературная плазма в таких установках, как ТОКАМАК и СТЕЛАРАТОР.

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

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

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

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

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

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

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

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

  1. можно предположить, что электроны по скоростям распределены по закону Максвелла:

                  (12)

 

 

Такое же предположение можно сделать и относительно лёгких ионов, поэтому в выражении (12) приняты следующие обозначения:

Ni — концентрация i-ой компоненты частиц,

Mi —масса частиц i-ой компоненты.

  1. в возбужденном состоянии распределение частиц, очевидно, определится формулой Больцмана:

                               (13)

 

 

здесь N0, g0 - концентрация и статистическая масса для нормального состояния.

NS, gS - концентрация и статистическая масса в возбужденном состоянии.

eS - энергия возбужденного состояния.

в) распределение атомов по стадиям ионизации определяется формулой Саха:

            (14)

 

Ne, Ni, Na - концентрации соответственно электронов, ионов возбужденных атомов.

Vj - энергия ионизации.

 

Давление в дуге значительно превышает атмосферное, поэтому предположение о термодинамическом равновесии оказывается достаточно правомочным. Следовательно, можно считать, что температуры - электронного, ионного, атомного газа равны между собой, т.е. Te@Ti@Ta.

 

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

 

Измерения показывают, что при токах (200 - 600) А (наиболее типичная дуга), Т = (11 000 - 16 000) К, при токах (20-60) А в воздухе при использовании угольных электродов Т» 6300 К.

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

  • интенсивность излучения: 

            (15)

 

 

здесь DE - излучаемая энергия через единичную площадку DS , при единичном телесном угле DW, в единичном диапазоне длин волн Dl, за единицу времени Dt.

  • излучательная способность:

        (16)

 

где DV - единичный объем.

Бесспорно, что должно быть справедливо соотношение: , здесь Dl - толщина излучаемого слоя, см. рис. 4.

 

Количество световых квантов с энергиями hn, излучаемое единичным объемом в геометрии 2p стерадиан, пропорционально вероятности перехода атома из состояния n в состояние m. Таким образом, выражение , определяет число атомов nm, имеющих такие переходы или находящиеся в состоянии m .

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

           (17)

 

Следовательно, интенсивность, излучаемая из плазменного цилиндра с площадью DS=1 и длиной Dl можно определить из соотношения:

 

 

 

           (18)

 

 

 

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

      (19)

 

Здесь  - статистическая сумма по возбуждённым состояниям, em - энергия возбужденного уровня m. Интенсивность излучения сплошного спектра в расчете на единичные параметры позволяет определить коэффициент излучения в виде:

 

              (20)

 

 

Если поглощение не велико, то можно записать приближенное соотношение:

 

                            (21)

 

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

Формула Крамерса справедлива до некоторой граничной частоты ng, что позволяет записать соотношения:

 

                 (22)

 

 

 

для частот n>ng

 

   (23)

 

 

 

где:     Ne  - концентрация электронов

Ni  - концентрация ионов 

ng  - граничная частота

z    - заряд иона до захвата электрона

s    - заряд иона после захвата электрона

 

 

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

 

         (23a)

 

 

Длина l  - в данном случае, определяется при томографической реконструкции параметров: T(x,y) ~  J(x,y) (см. полный конспект лекций).

Заметим, что статистическая сумма U0(1) по всем состояниям излучающего атома является функцией температуры - U0(1) (T). Если условие термодинамического равновесия можно считать выполненным, то n0(T) - концентрация частиц в основном состоянии может быть определена из системы уравнений, включающих — формулу Саха (14), закон Дальтона P = k T N, условие квазинейтральности: (n+ = n). Следовательно, экспериментально измерив, интенсивность J спектральной линии, и вычислив l, можно определить T(xi,yi).

Наибольшая сложность здесь заключается в определении величины  , но ее можно “обойти”, воспользовавшись эталонным источником излучения, абсолютная температура которого и вид спектра известны нам apriori. Такими источниками, например, являются ленточные лампы, которые градуируются по яркостным температурам (при l=0,655 мкм).

 

Относительная погрешность при измерении методом абсолютных температур равна: 

     (24)

 

 

 

 

Заметим, что в приближении Абеля погрешность составляет (8 -10) %, при томографической реконструкции относительная погрешность уменьшается до значений (0,4 -1) %.

 

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

 

ЛЕКЦИЯ 3.

 

ПРОХОЖДЕНИЕ РЕНТГЕНОВСКОГО И ГАММА - ИЗЛУЧЕНИЙ ЧЕРЕЗ ВЕЩЕСТВО

 

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

 

      (25)

 

 

где:

J0 - интенсивность падающего плоскопараллельного пучка,

J(x) - интенсивность пучка, прошедшего “толщину” образца,

mx - линейный коэффициент ослабления излучения, показывает при какой толщине данного материала, исходная интенсивность уменьшается в “е” раз,

Dx -толщина барьера (см. рис. 5).

 

 

При прохождении пучка g-квантов через вещество наблюдается множество эффектов взаимодействия их с атомами и электронами материала образца, однако наибольший вклад в ослабление падающего потока дают : фотоэффект, комптоновское рассеяние и эффект образования пар. Заметим, что при прохождении g - кванта через материал образца, он может полностью потерять свою энергию — полное поглощение, или ее часть ее — при эффектах рассеяния. В том и другом случае они «считаются выбывшими» из плоскопараллельного пучка. Таким образом, та часть фотонов, которые не провзаимодействовали с атомами и электронами материала образца, либо слабо провзаимодействовали, т.е. отклонилась на небольшие углы, образует поле излучения за барьером, которое несет информацию о характеристиках  материала, например, плотности, наличии локальных включений и т.д.

Рассмотрим процессы взаимодействия g-квантов с атомами вещества более подробно.

 

Фотоэффект — фотоэлектрическое поглощение g-кванта атомным электроном, при этом энергия падающего фотона расходуется  на “вырывание” электрона из атома. Энергия этого электрона       будет равна:

 

,                             (26)

 

где  -  энергия g - кванта,  энергия связи электрона в атоме. Из (26) ясно, что фотоэффект возможен при значениях энергии: , поэтому вероятность фотоэффекта  или сечение взаимодействия, имеет резкие скачки при энергиях , равных энергиям ионизации K-, L-, M-, ... оболочек. Заметим, что сечение взаимодействия одного g - кванта с одним атомом - s пропорционально вероятности этого процесса. Можно показать, что линейный коэффициент ослабления m=sN (где N - число атомов в 1 см3), тогда величина  приобретает смысл среднего свободного пробега g - кванта в веществе. Расчеты показывают, что фотоэффект происходит главным образом на К- оболочках.

Можно показать, что сечение фотоэффекта на К- оболочке:

                     (27)

 

 

где А — const, z — заряд ядра.

 

Несложный анализ (27) говорит о том, что для тяжёлых элементов, например Рв, фотоэффект имеет заметное значение даже при энергиях .

 

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

Угловое распределение фотоэлектронов зависит от их энергии. При малых энергиях (десятки пико - электрон-вольт) фотоэлектроны испускаются в направлении, перпендикулярном пучку g-квантов, с ростом энергии угол уменьшается, например: .

 

Комптоновское рассеяние

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

 

Воспользовавшись законами сохранения энергии и импульса, можно получить выражения для энергии g-кванта после рассеяния - , энергии электрона после взаимодействия -, так называемого Комптон - электрона и углами q и j (см рис. 6), запишем выражения:

     (27’)

 

 

 

Дифференциальное сечение рассеяния для Комптон - эффекта , (где  - телесный угол, в котором рассеиваются фотоны), можно рассчитать с помощью формулы Клейна-Нишина -Тамма:

    (28)

 

При малых энергиях, т. е.  формула (28) упрощается до выражения: 

     (29)

 

 

Эффект образования пар

 

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

,             (30)

 

здесь  - полные энергии электрона и позитрона,

 

EA- энергия ядра отдачи (или электрона отдачи), ее величину можно оценить с помощью формул:

,

 

или

 

  (31)

 

 

здесь М - масса ядра,  А - порядковый номер элемента. Сечение взаимодействия – для пар, определяется как:

     (32)

 

 

Сечение образования пар  довольно сложно зависит от энергии (см. рис.7). При больших энергиях, его величину можно оценить с помощью выражения (32)

 

Энергия g - кванта почти поровну распределяется между электроном и позитроном, которые разлетаются под углом p по отношению друг к другу.

 

 

С точки зрения регистрации гамма квантов по эффекту образования пар, надо иметь в виду, что позитроны  через очень короткое время аннигилируют,  испуская два кванта с энергиями 0,511 МэВ и углом разлета p.

На рис.7 изображены зависимости коэффициентов ослабления рассмотренных эффектов mф, mк, mп в зависимости от энергии гамма - квантов для типичных элементов Al, Pb с малыми и большими значениями A и Z.

 

Нетрудно показать, что результирующий коэффициент линейного ослабления определится как:

        

    (33)

 

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

       (34)

 

r -  плотность материала; если материал неоднороден, то r(x, y, z,). Понятно, что для "послойной"  неоднородности материала справедливо выражение:

 

            ,  (35)

 

где  -  массовые коэффициенты ослабления в материале i-го слоя, mi -относительные массовые количества.

 

Если распределение плотности материала произвольно по объему, то в выражении (35) очевидно следует перейти к пределу при .

 

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

 

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

 

или:

   (36)

 

 

 

б) для упрощения расчетов, внекоторых случаях, фактор накопления можно представить:

 

      (37)

 

 

где  a1; a2; A2 = (1-A1) - коэффициенты, не зависящие от mх. Индекс ¥ означает, что выражение (37) получено в приближении, что Dх - толщина барьера (образца) стремится к ¥.

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

 

       (38)

 

 

Заметим, что для органических соединений фактор накопления В=(1¸2), для элементов, имеющих малые и средние значения z, В=(2,2¸6).

 

 

ЛЕКЦИЯ 4.

Ядерный магнитный резонанс

 

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

Рассмотрим основные положения и  характеристики атомных ядер. Напомним, что ядро состоит из протонов и нейтронов (нуклонов), они отличаются двумя различными зарядовыми состояниями (p + e = n). Количество протонов  в ядре называется порядковым номером Z, количество нейтронов в ядре обозначают N. Для легких ядер N/Z » 1,  к концу периодической таблицы Менделеева это отношение возрастает N/Z » 1,6. Полное число нуклонов A=N+Z (массовое число), ядра, имеющее одно и то же Z при разных А, называют изотопами, ядра, которые имеют одно и то же А при разных Z называются изобарами, соответственно обозначения ZXA;AXZ,  где х - символ химического элемента. Сегодня известно около 300 устойчивых и свыше 1000 неустойчивых ядер. Ядра не имеют ярко выраженной границы, примерно оценить радиус можно с помощью эмпирической формулы R=R0A1/3; R0~(1,3¸1,7)´ 10-10 м. Объем ядра пропорционален числу нуклонов, а плотность ядерного вещества rя» 1011 кг/см3.

Энергия связи нуклонов в ядре равна разности энергий протонов и нейтронов в ядре и их энергий в свободном состоянии, т.е.:

 

    (39)

 

Нетрудно понять, что E(A,Z)<0 и численно определяется той работой, которую необходимо совершить, чтобы расщепить ядро на свободные нуклоны. Очень часто под энергией связи понимают положительное число Е. Под дефектом массы понимают D=М-А или , где М - масса атома, А - массовое число. Если ядро имеет наименьшее значение энергии  Еmin из возможных квантовых значений энергий, то говорят,  что оно находится в основном состоянии, если Е > Emin, то ядро возбуждено, случай Е=0 соответствует диссоциации ядра на составляющие n, p. Между нуклонами действуют ядерные и кулоновские силы, причем ядерные силы имеют неэлектрическую природу, т.к. установлено, что силы взаимодействия n « p, p«p, n«n одинаковы.

 

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

Рассмотрим магнитные и электрические свойства ядер. Нуклоны в ядре обладают орбитальным и спиновым магнитным и механическим моментами, являясь фермионами, их спин . Орбитальный и спиновый магнитные моменты соответственно равны: 

 

 

 

      (40)

 

Где l, s - соответственно орбитальные и спиновые квантовые числа, а  соответствующие гиромагнитные отношения:

 ;

 (41)

 

 

В соотношениях (41) используется постоянная:

 — ядерный магнетон.

 

Заметим, что в ядерной физике спином ядра называется его полный момент импульса, он геометрически складывается из полных моментов нуклонов, составляющих ядро. Внутренние квантовые числа нуклонов складываются алгебраически и дают суммарное число (0,1,2,3...) при четном А; если А - нечетное, то суммарное квантовое число будет полуцелым:  (1/2, 3/2, 5/2...). Понятно, что в первом случае ядра подчиняются статистике Бозе - Эйнштейна, во втором — Ферми-Дирака. Полный момент импульса для каждого нуклона определяется с помощью соотношения: . Спин ядра  и его магнитный момент  соответственно равны:

 

          (42)

 

 

Соответствнно полное внутренне квантовое число можно определить:

 

    (42*)

 

 

В (42) g - величина, аналогичная фактору Ланде для атома.

 

Поведение ядер во внешнем магнитном поле определяется значением их магнитного момента. Если ядро состоит из четного числа протонов и четного числа нейтронов, то оно называется четно-четным, если эти числа нечетные, то и ядра нечетно-нечетные. Магнитные моменты четно-нечетных ядер равны нулю. Магнитные моменты четно-нечетных и нечетно-четных ядер могут быть рассмотрены в простейшей однонуклонной модели ядра. В ней предполагается, что магнитные моменты таких ядер обусловлены движением одного “валентного” нуклона около остальной части  ядра, состоящего из четного числа нуклонов, их векторная сумма моментов (орбитальных и спиновых) равна нулю. В этом случае можно записать:

                                 (42**)

 

 

Модули магнитных моментов определяются из выражений:

     (43)

        (43*)

 

 

 

 

 

Электрический заряд нуклонов (протонов) распределен в ядре асимметрично. Мерой отклонения такого распределения от сферически симметричного является квадрупольный электрический заряд момента ядра . Если представить распределение заряда в виде эллипсоида вращения, то квадрупольный момент определяется:

                       (44)

 

где а, b - полуоси эллипсоида.

 

Дипольный электрический момент ядра в его основном состоянии равен нулю. При наложении внешнего магнитного поля происходит квантование спина ядра, и каждый уровень расщепляется на (2J+1) подуровня (Зеемановское расщепление ядерных уровней).

 

Определение: Ядерным парамагнитным резонансом называется избирательное поглощение электромагнитного излучения веществом, связанное с переходом его ядер между различными зеемановскими подуровнями энергии. Резонансные частоты для переходов подчиняются правилам отбора для магнитного (внутреннего) квантового числа , и соответственно равны:

,     (45)

 

 

здесь: 

g - фактор расщепления (аналогичен фактору Ланде)

 

 - ядерный магнетон,

 

H - напряжение внешнего магнитного поля.

Частоты ЯПР при одном и том же значении напряженности внешнего магнитного поля H в 10 000 раз меньше частот ЭПР и лежат в области  Гц при  Эрстед.

Если квадрупольный момент ядра  то ядро взаимодействует с внутримолекулярным и внутрикристаллическим полем, в результате такого взаимодействия наблюдается штарковское расщепление уровней. Например, если внутреннее (статическое) поле обладает  аксиальной симметрией, то энергетический спектр штарковских линий (подуровней) имеет вид:                  

       (46)

 

 

здесь:

 - магнитное квантовое число ядра,

 

J - его спин,

 

 - вторая производная от внутреннего поля (потенциала поля) в области ядра,

 

OZ - ось симметрии этого поля.

 

Ядерным квадрупольным резонансом называется избирательное поглощение энергии внешнего электромагнитного поля, связанное с переходами его ядер между штарковскими подуровнями. Если эти переходы удовлетворяют правилу отбора  для J-1, J-2, ... , -J, то условие возникновения ядерного квадрупольного резонанса  с частотой  можно определиь из выражения:

 

          (47)

 

 

 

Спектр (47) состоит из серии равноотстоящих линий. Наибольшая частота спектра соответствует  и равна:

           

                                       (48)

 

 

Именно ЯКМР (nqr) является наиболее эффективным способом определений структур молекул,  кристаллов, дает возможность получить разрешение до 10-3 м, в медико-биологической диагностике при помощи компьютерно – томографической реконструкции.

 

 

ЛЕКЦИЯ 5.

 

Методы и алгоритмы реконструкции 2-D томографии

 

Методы реконструктивной (вычислительной) компьютерной томографии условно можно разделить на три группы:

 

I. ART (Algebraic Reconstruction Technology) — методы алгебраической реконструкции.

 

II. CRT (Convolution Reconstruction Technology) — методы реконструкции с использованием интеграла свертки).

 

III. FRT (Fourier Reconstruction Technology) — способы реконструкции с использованием методов Фурье-синтеза.

 

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

 

Алгебраические методы реконструкции

 

             Алгебраические методы восстановления функций-сечений имеет самую длительную историю, алгоритмы реконструкции у них достаточно просты, в основном сводятся обычно к решению системы линейных уравнений, они наиболее просты для понимания. Главный их недостаток — громоздкость вычисления. Действительно, для того, чтобы восстановить матрицу 256 x 256 при 256 уровнях квантования амплитуды функции проекции, необходимо произвести » операций умножения, что занимает достаточно большое время.

  Суть алгебраических методов реконструкции сводится к следующему: предположим, что необходимо двумерное дискретное сечение, состоящее из m-строк и n-столбцов

 

           

 

(для упрощения полагаем m=n), каждый элемент которого имеет значение функции амплитуды , (см. рис. 8).

 

Таким образом, неизвестная функция f(x,y) представляется в виде дискретного изображения, содержащего (n x m) элементов, размеры каждого элемента (Dl x Dl). Для того, чтобы получить дискретные функции проекции, имеющие одинаковое количество отсчетов для любого угла qi , восстанавливаемую функцию достаточно поместить в круг, значения отсчетов в сегментах S0 положить либо равным нулю, либо считать их достаточно малыми. Иногда, для функций, обладающих свойствами аксиальной симметрии, особенно на периферии, например плазменный факел, поток высокотемпературного газа на срезе сопла и т.д., целесообразно воспользоваться внутренним кругом диаметром nDl, что естественно сократит число исходных уравнений.

 

Функции проекции получаются суммированием m ij  по данному лучевому направлению qi. «Размер» отсчета на проекции очевидно должен быть равен элементу реконструируемой функции, а именно Dl. Это накладывает некоторые дополнения в построение исходных уравнений, а именно: необходимо учесть при данном значении qi  вклад каждого элемента   mij  при проецировании его лучом шириной Dl. Проще всего это можно сделать введением коэффициента aij . учитывающего соотношение площадей пресекаемой лучом и площади D l2  данного элемента mij. Таким образом, для каждой функции-проекции можно записать уравнение в виде

                    (49)

 

здесь aij  - известный коэффициент, определяющий i, j -ячейки в i-й отсчет проекции, а  mij  - неизвестное значение функции в соответствующем элементе. Если число проекций будет равно n×m=n, то мы получим n2 линейных уравнений. Решая эту систему известными способами, например, с помощью метода обращения, можно восстановить неизвестные значения   m ij , т.е.:  f(x,y).

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

Перечисленные сложности в некоторой степени могут быть устранены с помощью итеративных методов реконструкции. Стратегия такого подхода заключается в следующем: производится последовательная коррекция приближенных значений  m ij в искомой функции f(x,y) для их согласования со значениями измеренных лучевых сумм в проекциях. Обычно здесь придерживаются следующего порядка: сначала задается условный начальный набор значений mi0 , например постоянная по изображению плотность, затем по этим значениям рассчитываются проекции. Если полученная расчетом лучевая сумма меньше ее измеренного значения, то значение функции в каждой ячейке, дающей вклад  в данную лучевую сумму, увеличивается на определенную величину. После проведения такой процедуры для всех ячеек и всех лучей, первую итерацию считают выполненной.

            Количество последовательных  итераций определяется степенью необходимости точности и скоростью убывания ошибки. Обычно весовые коэффициенты aij  не вычисляются по строгим  формулам и принимаются равными либо 1 либо 0 в зависимости от того, лежит ли центр ячейки i,j в полосе Dl или нет.

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

 (50)

 

здесь  значение функции в i, j ячейке до очередной итерации,

 - значение функции там же после итерации,

 - коррекция, “вносимая” в i, j ячейку от j-го луча.

Обычно используется два варианта коррекций : аддитивный и мультипликативный. В первом случае коррекция пропорциональна величине коэффициента  a ij , то есть можно использовать соотношение:

       (51)

 

 

здесь  разность между измеренной и рассчитанной после k-ой итерации лучевых сумм.

В случае мультипликативной коррекции величина  пропорциональна , полученной при предыдущей итерации, т.е.:

 

           (52)

 

 

В зависимости от организации процедуры итераций обычно выделяют три наиболее употребительных метода:

- одновременная коррекция (ILST - Iterative Least Squares Technique)

 

- поточечная коррекция (SIRT - Simultaneous Iterative Reconstruction Technique)

 

- получевая коррекция (ART - Algebraic Reconstruction Technique).

 

В технике расчета ILST коррекция рассчитывается одновременно вначале итерации, а затем корректируются сразу все ячейки, при этом число необходимых умножений на одну итерацию составляет pn3. При методе SIRT каждая ячейка корректируется одновременно для всех лучей, проходящих через нее, причем при коррекции последующих точек учитываются результаты предыдущих коррекций. Число необходимых умножений на одну итерацию здесь составляет p2n4 , что оказывается большим, чем в методе обращения матрицы. При реконструкции методом ART коррекция осуществляется одновременно для всех точек, дающий вклад в отдельный луч, затем процедура повторяется для следующего луча и т.д., при этом учитываются результаты предыдущих итераций (число умножений в этом случае составляет pn3), в смысле затрат машинного времени данный способ самый быстрый из алгебраических способов восстановления. Несмотря на приближенный характер решений получаемых при каждой итерации, при достаточном числе итераций, алгебраические методы не уступают в точности аналитическим методам, но требуют больших временных затрат.

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

 

ЛЕКЦИЯ 6.

 

Методы реконструкции, использующие интеграл свертки

 

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

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

 

1. N-мерное непрерывное преобразование Фурье

 

Рассмотрим некоторую функцию непрерывных переменных f (x1,x2...xn), совокупность которых удобно представить в виде n-мерного вектора, т.е. x1, x2 ... xn Þ , соответственно саму функцию обозначить как: .

 

N-мерное преобразование Фурье соответственно обозначается:  и определяется выражением:

 

                        (53)

 

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

 

Соотношение (53) называют прямым преобразованием Фурье, оно “переводит” функцию из пространства сигналов в «пространство частот».

Обратное преобразование Фурье имеет вид:

               (54)

 

 

 

Отметим важное свойство пары преобразований (53), (54), используемое в компьютерной томографии:

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

 

Заметим, что любое ортогональное преобразование не изменяет координаты, а лишь определяет их значение относительно “повернутой системы” К` (см. выражения (1) и (2)), аналогичный поворот координат происходит в пространстве частот.

 

 

2. Дискретное преобразование Фурье (ДПФ)

 

 

Теорема 1 Ширина спектра функции  в пространстве сигналов считается ограниченной, если существует N переменных в частотной области , таких, что:

 при  (i=1,2,3,...,N).

 

Совокупность  переменных  называют векторной шириной спектра.

 

 

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

 

Обозначим через  дискретную функцию соответствующую , имеющую шаг отсчетов , т.е.:

 (55)

 

 

Дискретное преобразование Фурье в этом случае будет иметь вид:

,   (56)

 

где:

 

 

 

 

Аналогичное обратное ДПФ: 

       (57)

 

 

 

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

           

где:

        (58)

 

 

 

Свойство 1  Если количество отсчетов, отличных от нуля для функции  конечно, то ее Фурье-образ может быть представлен конечным количеством отсчетов.

 

Обозначим, декартовы отсчеты функции  через  значения: , тогда дискретная функция  может быть определена:

           

      (59)

 

 

 

здесь ki - целые числа, такие, что выполняются условия:

 

       если Mi - четное

 

          если Mi - нечетное

 

 

Таким образом:

(60)

 

и

 

 (61)

                                  

 

Понятно, что функция  является периодической по параметру ki с периодом , следовательно, вектор  можно определить как:

           

 

 

 

и тогда выражения (60) и (61) можно записать в виде: 

                            (62)

 

      (63)

 

 

 

Важное примечание: Выражения (62) и (63) являются парой N-мерных ДПФ. Эффективным способом вычисления (62) и (63) являются алгоритмы быстрого преобразования Фурье БПФ.

 

3.Теорема о проекциях и сечениях

 

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

 

,                  (64)

 

где  - соответствующие гиперплоскости,  - ортогональная матрица поворота, следовательно, функция проекции может быть определена:

                 (64’)

 

 

Теорема. Фурье - образ проекции , т.е. проекция в частотной области представляет собой сечение Фурье-образа , получаемое при .

 

Поясним эту теорему на примере двумерной функции, (см. рис. 9). Если от исходной функции, получить проекции под заданными углами qi, то Фурье-образ от любой проекции Fi(Pi) является центральным сечением двумерного Фурье-образа F (wx,wy),  вычисленного  от функции f(x,y) под тем же углом qi ( в данном случае угол отсчитывается от оси w1(wx) .

 

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

 

 

 

Заметим, что, используя “напрямую” теорему о сечениях, можно построить алгоритм FT-реконструкции, который должен содержать следующие операции:

 

1. получение одномерных функций проекций;

 

2. расчет одномерных Фурье - спектров или функций проекций;

 

3. ”формирование” двумерного Фурье-спектра, который оказывается определенным в полярной системе отсчета;

 

4. ”пересчет” полярных отсчетов в декартовы;

 

5. вычисление обратного двумерного преобразования Фурье.

 

Однако алгоритм, реализованный по пп.1-5 малоэффективен, т.к. возникают чисто алгоритмические сложности при переходе от полярных к декартовым координатам, а вычисление обратного 2-хмерного ДПФ должно содержать n2 одномерных преобразований для матрицы n ´ n.

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

 

Из геометрии получения исходных данных понятно, что вращение объекта относительно оси, проходящей через 0’ относительно неподвижного источника излучения и системы детекторов определяемых в системе координат X0Y, эквивалентно вращению системы векторов и источника излучения вокруг той же оси относительно неподвижного объекта. Искомой функцией в том и другом случаях является f(x,y), которая связана с f(x`,y`) матрицей поворота. Если для данного направления q выбрать полярную систему координат, чтобы ось  , то координата каждой точки будет определена с помощью вектора  , задаваемого относительно оси 0r (или 0`r`) углом (j-q) и модулем .

Используя понятие массового коэффициента ослабления m(x`,y`) и плотности материала в данной точке r (x`,y`) искомую функцию f(x`,y`), а следовательно и f(x,y) можно определить как:

                   (65)

 

Для каждой регистрируемой проекции под углом q интенсивность прошедшего через объект излучения, очевидно, может быть определена как:

   (66)

 

 

Элементарно преобразовав выражение (66) под “теневой” проекцией будем понимать функцию вида:

 

         (67)

 

 

Вычислим одномерные преобразования Фурье для каждого угла q, понимая, что координата x` принимает значения от -l до +l (или от -r до +r) , т.е. для данного значения угла q:

       (68)

 

 

Заметим, что без учета влияния рассеянного излучения диаметрально противоположные проекции одинаковы, следовательно, для того, чтобы получить, двумерный Фурье-спектр, угол q достаточно изменять в пределах от 0 до p, но l при этом изменяется от -r до +r (для Фурье-образа угол q сохраняется, R=2p/r). Но в силу симметрии рассматриваемого подхода в Фурье плоскости удобнее изменить пределы интегрирования, т.е. 0£R£¥, 0<2p. Тогда обратное двумерное преобразование Фурье будет иметь вид:

   (64)

 

 

Из рис.9 видно, что для выбранной геометрии значение l в теневой проекции g(l;q), которое определяет вклад точки  очевидно определится как , переопределим теневую проекцию как:

        (70)

 

 

Выражение (69) по сути дела определяет некоторый алгоритм реконструкции исходной функции, полученный в полярных отсчетах (или координатах). Нетрудно видеть, что восстановление изображения f(r,j) имеет специфические “смазы” вблизи точки, соответствующей оси объекта и возрастающую размываемость к периферии.

 

Эти артефакты (искажения) заложены самой геометрией реконструкции. Действительно плотность значений функции f(r,j) в центре максимальна и убывает с ростом r. Таким образом, для того, чтобы функция f(r,j) максимально соответствовала искомой f(x,y) или f(x`,y`) необходимо модифицировать метод, например, изменить функцию проекции g(rCos(j-q)) так, чтобы интеграл в (70) максимально приближался к функции распределения плотности. Это можно реализовать следующим образом:

 

 

или

 (71)

 

 

NB  Выражение (71) записано для одной проекции под данным углом q, об этом говорит знак “;” при описании параметров.

 

  (72)

 

 

Выражение (69) также можно переписать подобным образом, изменив соответствующие пределы интегрирования:

 

.  (73)

 

 

Интегралы (72) и (73) стремятся описать , но отличаются друг от друга множителем .

 

Следовательно, для того, чтобы скорректировать функцию проекции , введем “новую” функцию , которая должна быть определена следующим образом:

           (74)

 

тогда, очевидно можно записать:

                  (75)

 

 

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

Так, исходные данные определяются функцией g(l,q). Получим исходные данные перехода к g`(l,q); для этого выполним обратное преобразование Фурье, использовав выражение (74), то есть:

           (76)

 

 

Важные замечания:

 

 

 

 

Предположим, что это  удалось сделать, т.е. получилось вычислить , или:

           

 

 

Тогда можно воспользоваться теоремой о свертке, а именно:

               (78)

 

 

Таким образом, мы получили выражение для f(r,j) максимально приближенное к искомой функции f(x`,y`).

 

На первый взгляд, кажется, что задача решена, действительно, достаточно вычислить интеграл свертки (78), осуществить обратное проецирование (75) и в полярных координатах мы получим f(r,j)Þf(x,y), но вся сложность в том, что интеграл (77) вычислить невозможно, т.к. при R®¥, функция  терпит разрыв.

 

Существует несколько способов «обхода» данного обстоятельства, в частности бесконечные пределы интегрирования можно заменить конечными, которые определяются, например, шириной спектра реконструируемого изображения ±А/2, т.е.:

                 (79)

 

 

Следовательно, функциюможно представить в виде ряда Фурье с периодом А, т.е.:

 

;               (80)

 

.    (81)

 

 

 

Если A достаточно велико  и, следовательно, вычислив интеграл в (81) имеем: 

        (82)

 

 

 

 

 

Таким образом, если имеем дискретный набор исходных данных g(l;q), при l=ma, то вычисление интеграла в (78) заменяем суммированием:

,  (83)

 

 

представляющее собой, по сути, дискретную свертку. Или, воспользовавшись выражением (82), перепишем (83) в виде: 

       (84)

 

 

 

Хотя из выполненных выше операций это не очень очевидно, тем не менее, можно показать, что формула суммирования (84) не является точныи приближением для интеграла (78) в данной точке.

 

Но в то же время если:

 

Лекция 7.

 

 

Восстановление искомых функций по конечному числу проекций

 

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

а) Конечное число проекций, получаемых при равных углах поворотом q i  при изменении q i в пределах
     0
£ q £ p;  0 £ q £ 2p .

 

б) Ограниченное количество функций проекций при углах конвергенции: .

 

в) Получение функций проекций при неравных углах поворота, т. е.:

             

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

Проанализируем возможности реконструкции для таких случаев. Считаем, что искомую функцию необходимо определить в декартовых координатах f(x,y), а ее Фурье - образ может быть определен в полярных координатах:

                              (85)

             

 

              Предположим, что эта функция аналитична в параметрах (w,q). Если это не совсем так, то можно воспользоваться свойствами симметрии f (r,(q - j)) и F (w,q).

 

Действительно, для большинства объектов реальные части Фурье преобразования удовлетворяют условию:

 

                                     (86)

 

В этом случае мнимые части, также удовлетворяют условию:

                                  

                                    (87)

 

таким образом, Фурье-образ может быть представлен набором Фурье - образов проекций, т. е.:

 

                 (88)

 

Для конечного числа проекций последнее выражение можно записать в виде:

 

                   (89)

 

 

Если данные F(R,q) определены в ограниченной области q, т. е. , то коэффициенты am, bm, cn, dn в (89) могут быть определены следующим образом:

·         Предположим, что данные (отсчеты) получены в точках отсчетов 0,1,2,...L вдоль «Фурье» - радиуса w в виде набора , это дает возможность получить коэффициенты am, bm, cn, dn из условия:

          

                 (90)

             

 

              Понятно, что условия (90) разрешимы, если число угловых отсчетов равно (2М+1) и 2(N+1) соответственно. Следовательно, области Фурье - спектров, после вычисления am, bm, cn, dn, могут быть определены следующие значения составляющих спектров:

 

                  (91)

             

 

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

             

              Проанализируем более эффективные алгоритмы реконструкции, так же основанные на алгоритме свертки и обратного проецирования. Полагаем, что в качестве реконструируемой функции f(x,y) выбрана функция распределения плотности r (r,j) - в полярных координатах. Обратное преобразование Радона для r (r,j) можно представить в виде:

 

                  (92)

 

Представим действие «обратного» оператора Радона R-1 в виде последовательности простых операторов:

 

- обозначим взятие частной производной от функции проекции по переменной l:

                       

                                      (93)

 

 

- воспользуемся преобразованием Гильберта для переменной l

                      

                     (94)

 

 

- введем оператор обратного проецирования:

 

                      (95)

 

 

- тогда обратное преобразование Радона можно записать в виде последовательного действия трех операторов:

                                              

                    (95*)

 

Примечание:

              Интеграл в выражении (94) является несобственным, действительно, при (l = t) он расходится. Его можно вычислить, преобразовав в интеграл в смысле главного значения Коши:

                       

                                   (96)

 

             

 

              Вычислить данный интеграл можно, если представить преобразование Гильберта в (94) в виде свертки двух функций и функции  , т. е.:

 

                               (97)

 

             

 

              «Напрямую» операцию (97) выполнить невозможно, поэтому аппроксимируем функцию c(l) некоторой функцией cА(l), так, чтобы выполнялось условие:

                                  

                                   (98)

             

             

 

              Такой подход (98) носит название метод регуляризации, а множество - функций называется семейством регуляризирующих функций.

 

Примечание:

              Методы регуляризации наиболее детально проработаны для задач компьютерной томографии академиком Тихоновым А. Н. и профессором Арсениным В. Я. Полученные в результате алгоритмы реконструкции оказались весьма эффективными.

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

                            (99)

 

где

 

                                     (100)

                                            (101)

 

 

В выражении (99) величина:  - оператор обратного преобразования Фурье от произведения спектров.

 

Примечания:

 

 

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

                                              

                                             (102)

 

Это можно реализовать следующим образом:

                                (103)

 

Здесь   называют функцией «окна», которая очевидно должна обладать следующими свойствами:

             

                на участке  монотонно невозрастающая функция

                

 

 

где

             

Вычислим:

             

                                 (104)

 

 

Тогда с учетом выражения (103) запишем для функции  в виде:

 

                              (105)

 

 

Вычислим предел аппроксимирующей функции:

 

 

 

                 (106)

 

 

Следовательно, мы показали, что условие (98) выполняется, т. е.:

 

                                   (107)

 

Заметим, что искомая функция f(x,y) (или ) определена в конечной области - сечении объекта, следовательно, область ее существования можно задать в виде круга, ограниченного окружностью радиуса R, т. е.:

                       

                                       (108)

 

 

С учетом этого проинтегрируем правую часть (107) по частям, то есть:

                        .             (109)

 

Воспользуемся теперь соотношением  (105) и вычислим производную от функции cА под знаком интеграла, т.е.:

              .             (110)

 

Введем следующие обозначения:

                                   

                                               (111)

                      

 

              .        (112)

             

 

              Отметим, что формула (112) в соответствии с выражением (97) представляет собой аппроксимацию преобразования Гильберта при А®¥.

 

Воспользуемся оператором обратного проектирования, запишем:

              ,      (113)

 

 

здесь - оператор  определяется выражениями (111), (112), а значок ~ говорит о приближенности аппроксимации.

             

              Таким образом, мы проиллюстрировали один из возможных путей аппроксимации обратного преобразования Радона, который сводится к двум процедурам

  1. Функция проекции сворачивается с функцией h(l) определенной (111) и (112)
  2. Выполняется процедура обратного проецирования.

 

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

              Следовательно (94) эквивалентно значений свернутых проекций вдоль семейства проектируемых прямых под всеми углами и суммированию значений размазанных проекций в точках задания функции .

              Конкретный вид функции «окна» W(w) в выражении (103) нуждается в пояснении. Существует достаточно много «окон», удовлетворяющих поставленным выше условиям. Такие «окна» хорошо известны из теории цифровой обработки сигналов.

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

              В компьютерной томографии, когда число проекций достаточно велико 360, 720 и т. д. И число отсчетов в каждой проекции 512, 1024, 2048, наибольшее распространение получили следующие виды окон.

1) Прямоугольное «окно» (Rect.)

 

 

             

 

                 (114)

 

             

Такое окно дает ядро свертки Рамачандрама - Лакшминараньянана

Рис. 10

 
(115)

 

                                  

2) «окно» Шеппа - Логана

              (116)

 

 

ядро свертки, соответственно имеет вид:

                  (117)

 

 

3) «окно» Ханна

                       (118)

 

 

 

                  (119)

 

 

 

4) «окно» Хемминга

                                  (120)

 

 

 

 

соответственно, ядро свертки

 

 

                 (121)

 

 

 

 

 

 

 

 

 

5) «окно» Парсена

 

 

                         (122)

 

 

соответственно ядро

                 (123)

 

 

 

Здесь величина: .

 

              Кроме рассмотренных ядер свертки известны окна Баттерворта, Хорна и др. Их сравнение показывает, что лучшим разрешением для не зашумленных данных обладает ядро Рамачандрама, так как прямоугольное окно дает наиболее узкий центральный максимум ядра свертки. В то же время высокий уровень боковых лепестков преобразования Фурье от прямоугольного «окна» приводит к усилению шумов и появлению артефактов в области высокого контраста. Использование других окон диктуется необходимостью подавления шумовых составляющих, которые всегда присутствуют в проекционных данных. Однако при этом теряется разрешающая способность, за счет расширения центрального максимума ядра свертки. Поэтому в каждом конкретном случае в зависимости от характерных особенностей восстанавливаемых изображений подбираются оптимальные варианты.

           

Приближенная оценка показывает, что число умножений для алгоритма свертки и обратного проецирования составляет примерно N3, в то время как для алгоритма Фурье реконструкции число умножений порядка N2 (log2 N)2 (без учета интерполяционных процедур). Таким образом, сверточные алгоритмы оказались более «быстродействующими», имеют высокую устойчивую способность, обладают более простыми интерполяционными процедурами в сравнении с FT, ART - методами.

 

Phynist3d LabOFPT

 

© Copyright 2005-2011 LABOFPT - Look Team

All rights reserved.

Design – Phylonin O. V.