CC BY
63
УДК 620.179
ИСПОЛЬЗОВАНИЕ АЛГОРИТМА ОБРАТНОГО ПРОЕЦИРОВАНИЯ ДЛЯ РЕКОНСТРУКЦИИ СЕЧЕНИЯ ПРИ ИСПОЛЬЗОВАНИИ ОБРАТНО РАССЕЯННОГО ИЗЛУЧЕНИЯ
Б.И. Капранов, Х.-М.В.А. Кренинг, М.А. Мудров
Томский политехнический университет E-mail: [email protected]
Приведены результаты использования алгоритма обратного проецирования для целей математической реконструкции сечения объекта по интегральным данным, получаемым в геометрии одностороннего доступа с использованием обратно рассеянного рентгеновского излучения. Установлено, что разработанная геометрия сканирования позволяет получать интегральные данные по суммарному рассеянию с линии сканирования (альбедные суммы), являющиеся аналогом лучевых сумм в трансмиссионной томографии.
Ключевые слова:
Лучевая сумма, альбедная сумма, обратная проекция, комптоновское рассеяние, реконструкция, томограмма, артефакты.
Основным алгоритмом математической реконструкции распределения коэффициента линейного ослабления в сечении объекта в трансмиссионной томографии является алгоритм обратного проецирования [1, 2]. Он основан на получении в результате томографического эксперимента массива проекционных данных, представляющих собой лучевые суммы R (5, ф) (преобразование Радона), определяемые как интеграл от искомого распределения f (х, у) вдоль траектории луча l:
R(s,cp)= s*cos<p-y'sin<p,s*sm(p + y'cos<p dy',
где а - радиус окружности, охватывающей исследуемое поперечное сечение.
Проекции вычисляются под всевозможными углами ф и для тех значений s, при которых двумерная функция f (x, у) отлична от нуля. На практике величина s ограничивается физическими размерами исследуемого объекта, а угол ф изменяется в пределах от 0° до 9max. Чаще всего берут 9max = 360 .
По набору проекционных данных производится реконструкция распределения функции ф (x, у) в плоскости.
Суть алгоритма обратного проецирования в том, что оценку функции f (x, ) в любой точке находят путем сложения лучевых сумм для всех лучей, проходящих через данную точку.
Для реализации алгоритма обратного проецирования находят функцию F (r, ф) двух полярных переменных, значения которой в любой точке
г = у/х2 + у2, (р = arctg у/
вычисляются как
F{r,(p)= ^R{r -cos(9 -<p),Q)dQ
Капранов Борис Иванович,
д-р техн. наук, профессор кафедры физических методов и приборов контроля качества Института неразрушающего контроля ТПУ.
E-mail: [email protected] Область научных интересов: радиационный (альбедная томография) и ультразвуковой методы контроля. Крёнинг Ханс-Михаэль Вильгельм Адольф, руководитель Международной научно-образовательной лаборатории неразрушающего контроля Института неразрушающего контроля ТПУ. E-mail: [email protected] Область научных интересов: неразрушающий контроль. Мудров Михаил Анатольевич, инженер Аттестационного регионального центра специалистов НК Института неразрушающего контроля ТПУ. E-mail: [email protected] Область научных интересов: радиационный (альбедная томография) и ультразвуковой методы контроля.
В комптоновской томографии для получения информации о распределении плотности в сечении объекта используется рассеянное излучение.
Интенсивность рассеянного излучения зависит от плотности и от ослабляющих свойств материала. Существующие томографы на рассеянном излучении отображают распределение интенсивностей по локальным объемам (voxel). Типичные размеры воксела, используемые в томографе ComScan составляют 0,4*0,4*1,2 мм. Перемещение воксела по сечению контролируемого изделия синхронизируется с разверткой на экране монитора. В результате формируется сечение объекта.
Основными недостатками таких систем являются:
1. Низкая пространственная разрешающая способность, определяемая поперечными размерами воксела.
2. Слабая статистика потоков излучения, выходящих из малых объемов, что требует использования мощных источников рентгеновского излучения.
Одним из методов устранения этих недостатков комптоновской томографии является использование интегральных измерений потоков рассеянного излучения из всей линии рассеяния (луча), а не из отдельной точки линии (воксела), т. е. измерение альбедной суммы, являющейся аналогом лучевой суммы в трансмиссионной томографии. Статистика такого интегрального измерения намного превышает статистику потока, рассеянного из воксела. Выполнив сканирование в плоскости х, у, мы получаем массив альбедных сумм, по которым можно выполнить математическую реконструкцию распределения плотности в сечении.
В данной работе исследована возможность в условиях одностороннего доступа организовать процесс сканирования, обеспечивающий необходимое количество интегральных измерений (альбедных сумм) и возможность реконструкции распределения плотности в сечении с применением алгоритма обратного проецирования. Схема одностороннего сканирования показана на рис. 1.
Функция R (5, ф) представляет собой радоновский образ объекта и является распределением рассеянного комптоновского излучения в зависимости от положения щели коллиматора s при заданном угле ф. Из рис. 2 видно, что для каждого слоя объекта значение R(s, ф) пропорционально области объекта, в которой происходит рассеяние (LjHo). Поэтому при моделировании данных, получаемых при сканировании, удобно представить данные, как зависимость от длины L, равной длине хорды для круга R1 на расстоянии от центра s, за вычетом длины хорды для круга R2, при сканировании на отверстии.
R(s, <Р) = Н0 ■ L(s, (р),
где L{s,(p) = Ll{s,(p)-Ln{s,(p),Li - длина хорды для круга R1 на расстоянии от центра 5 (рис. 2); L2 - длина хорды для круга R2 на расстоянии от центра (S0^) - s) (рис. 2), где S0(ф) -расстояние от центра отверстия до линии соответствующей S^) = 0.
Рис. 1. Схема сканирования с получением альбедных сумм: 1, 2 - соответственно, первичный и вторичный коллиматоры; 3 - линия рассеяния; 4 - детектор; 5 - тест-объект
Рис. 2. Получение данных при сканировании
п
Из рис. 2 видно, что: 5 = г-совС/?); /? = — ср — а\ х = г-со8(а); у = г-$т(а).
Исходя из этого получается:
,-,.„.(г-<,+в))-г.*К„+в>--Г. =^п«+,.со5(„),
Бо (р) = х0-япр + у0-ад5р.
Выражение верно для любой точки (х, у), поэтому связь между координатами получаемых данных и координатами изображения описывается формулой:
£(<р) = х • ыгир + у • соър .
Полученные таким образом с помощью комптоновски рассеянного излучения данные о внутренней структуре объекта аналогичны получаемым при регистрации излучения, прошедшего сквозь объект. Как и в трансмиссионной томографии, необходимо решать обратную задачу - по массиву интегральных данных математически реконструировать детальное распределение плотности. Исходя из теоремы Котельникова, такая реконструкция требует, чтобы массив
экспериментальных данных имел дискретность (т. е. шаг перемещения) по выбранному направлению сканирования не более чем V от требуемого пространственного разрешения по этому же направлению. Поэтому в конструкцию механизма сканирования заложена возможность перемещения системы источник-детектор относительно объекта контроля с шагом 0,1мм и измерение в каждом положении интегрального сигнала с детектора.
Для реализации алгоритма обратного проецирования вычислялась величина \Вр\(х, j) по данным f, где f = R(si,pi) для М равномерно распределенных в пространстве параллельных лучей в каждом ракурсе. Пусть А обозначает угол между направлениями ракурсов (так что А = к/М), a d - шаг между параллельными лучами. Тогда
[Вр](х, v) = )R(s,<p)d<p,
0
{Вр\(х,у) = ’\R(x-sin(p + y-cos(p,(p)dcp .
0
Для вычислений правую часть выражения аппроксимируем суммой Римана для данного интеграла
М-1
A^R(x- sin(mA) + у • cos(mA),тА),
т=1
а затем проведем оценку для каждого значения т величины R (х ■ sin(»/A) + у ■ cos(tnA), тА) по известным значениям R(nd, тА) (-N < п < N) путем интерполирования.
При линейной интерполяции п выбирают так, чтобы
nd < х • sin (tilА) + у • cos (тА) < (п + 1 )с/. и вычисляют R (х • sin(/H/\) + у ■ cos(mА), тА) по формуле:
(и + 1W - х • sin (тиД) - v-cosfwA) л/ , x-sin(»7A) + v-cosfmA)-nd л//
---- ---------1-- — ------11 • R (nd, тА) +------ --- — ----- --- ----R((n +1 )d, тА).
d d
После вычисления значения \Вр\ (r, ф) для каждого элемента изображения производится
мультипликативная нормировка.
Таким образом, получаем следующий алгоритм действий:
а) Установка начальных параметров х, y0, ю.
б) Сбор данных по указанной ранее схеме.
в) Аппроксимация интеграла \Вр\ (r, ф) суммой Римана (по собранным данным).
г) Линейная интерполяция
sin(mA)
д) Вычисление суммы Римана для каждого элемента изображения (r, ф).
е) Мультипликативная нормировка полученного изображения.
Для проверки предлагаемого метода был разработан эксперимент, в котором первичный пучок коллимирован в виде ленты толщиной 0,5 мм и шириной 50 мм. Коллиматоры детектора имеют щели с раскрытием 0,5 мм и возможность поворота относительно оси сцинтиллятора для настройки на разную глубину в изделии. На рис. 3 показан тест-объект, его синограмма и томограмма.
в
Рис. 3. Тест-объект и результаты эксперимента: а) тест-объект; б) синограмма; в) томограмма
Полученный результат показывает возможность применения данного алгоритма при создании томографических систем на обратно рассеянном комптоновском излучении.
Исследованы также методы улучшения качества реконструированного изображения с использованием различных интерполяционных функций и модификация алгоритма обратного проецирования на основе выделения контуров.
Для интерполяции сигнала использованы: метод интерполяции по ближайшему значению. Интерполяционной функцией в этом методе является линейная функция, с помощью которой выполняется интерполяция значений многочленом нулевого порядка.
При интерполяции по ближайшему значению вычисляют Я(г со8(тП - ф), тП) по величинам Я(пй, тП), где й - шаг между параллельными лучами и п выбирают таким образом, чтобы величина | пй - г со8(тП - ф)| принимала наименьшее значение. Затем складываются лучевые суммы для лучей по одному из каждого ракурса, которые являются ближайшими к точке (г, ф), и результат умножается на П. Это демонстрируют результаты моделирования.
Метод линейной интерполяции. При линейной интерполяции п выбирают так, чтобы п<3 < г со$(тО - ф)<(п + 1)б/, и вычисляется Я(г со$(тО - ф), тП) по формуле
(Я + 1У-,-СОЕ(МАч»)д(тЛ ^ + ,-С05(МА — „)-П_1Щп + Ц(# тД) б/ б/
Линейная интерполяция является несколько более сложной: вместо лучевых сумм одного луча складывают линейную интерполяцию лучевых сумм двух лучей, которые находятся по обеим сторонам от точки (г, ф).
Лучшие результаты получаются, если использовать более сложные интерполяционные функции: Гаусса, 8тс(х), трехмерного 5-сплайна. На рис. 4 приведены примеры реконструкции
по методу выделения контуров с использованием интерполяционных функций Гаусса (а) и 8Шс(х) (б).
Рис. 4. Реконструкция по методу выделения контуров для интерполяции по Гауссу (а) и по sinc (х) (б)
Как видно, полученное изображение адекватно передает характерные особенности внутренней структуры объекта. Целесообразно использовать полученный алгоритм лишь как основу более сложного алгоритма, позволяющего получать более точные изображения. Представленные схема сканирования и модифицированный алгоритм реконструируют один слой объекта. Для получения полной томографической картины необходимо его применить к данным, получаемым из каждого слоя [3, 4].
Таким образом, предложенная схема сканирования для комптоновского томографа позволяет применять алгоритм реконструкции методом обратного проецирования в комптонов-ской томографии, что позволит повысить разрешающую способность метода.
СПИСОК ЛИТЕРАТУРЫ
1. Сизиков В.С. Математические методы обработки результатов измерений. - СПб.: Политехника, 2001. - 240 с.: ил.
2. Неразрушающий контроль и диагностика/ под ред. В.В. Клюева. - М.: Машиностроение, 1995. - 488 с.
3. Хермен Г. Восстановление изображений по проекциям: Основы реконструктивной томографии / Пер. с англ. - М.: Мир, 1983. - 349 с.
4. Кривовяз С.В. Схема сбора данных и алгоритм реконструкции по обратным проекциям для комптоновской томографии // Современные техника и технологии 2006: Сборник материалов конференции. - Томск, 2006. - С. 45-49.
Поступила 04.02.2013 г.
