CC BY
24
спективно применение ГРП с такими добавками воздуха при стробоскопической регистрации быстро-протекающих процессов с частотой до 100 кГц.
СПИСОК ЛИТЕРАТУРЫ
1. Гордеев В . А . , Даманский Е . А . Время' постановления искровых камер / Препринт ТИ, 1967,- С.6.
2. Мартынов В. И., Ры льдов В . В . Послеразрядная деионизация в искровых камерах с малым зазором//ПТЭ-I975.-T.3.- С.72-75.
3. F i s е h е J . and Zorn. Spark Chamber Symposium // RSI.- 1961.- V.32.- P.480.
4. Mijamoto S . Discharge Chmwer. Some Consideration of Its Characteristics // Nuovo Cimento, 1963,- V. 27 - №6,- P. 1325.
5. Беспалов В.И., Зайцев А.К., Кононов М.Ю., Кулешов В.К. Собственная нерезкость газоразрядных преобразователей, работающих с высокоэнергетическим тормозным излучением. // Дефектоскопия. - 1988. — № 1. -С. 71-78.
6. Самарский A.A., Гулин А . В . Численные методы.-М: Наука, 1989.-432 с.
7. Лоз а некий Э.Д., Фирсов О.В. Теория искры. -М: Атомиздат, 1975. -272 с.
8. A.c. 1208500 (СССР). Способ преобразования изображений в радиационном интроскопе с импульсным излучением / М.Ю. Кононов. В.Н. Ланшаков, В.К. Кулешов //Опубл. в Б.И., 1981.
9. Гороховский Ю.Н. Спектральные исследования фотографического процесса. -М: Изд-во физ.-ма т.лит.. 1960.-392 с.
10. Гороховский Ю.Н., Баранова В.П. Свойства черно-белых фотографических пленок. Сентрометрический справочник. -М: Наука, 1970. -388 с.
11. Фризер X . Фотографическая регистрация информации. Пер. с нем. / Под ред. К.В. Вендоровского. -М: Мир, 1978. -670 с.
12. Кулешов A.C., Кулешов В.К. Собственная нерезкость и частотноконтрастные характеристики газоразрядных преобразователей, работающих в лавинном режиме // Дефектоскопия. -1985.- №11.- С. 40-44.
13. Заневский Ю.В. Проволочные детекторы элементарных частиц. - М: Атомиздат, 1978.- 234 с.
УДК 620.179
Б. И. КАПРАНОВ, В. А. ШАВЕРИН, В. В. ВАРГА, Ю. В. АЛХИМОВ, О. А. СИДУЛЕНКО, В. Я. МАКЛАШЕВСКИЙ, В. Н. ФИЛИНОВ
АЛГОРИТМЫ И ПРОГРАММНОЕ ОБЕСПЕЧЕНИЕ КОМПТОНОВСКОГО ТОМОГРАФА "ТОМ8СА1Ч-200"
Создан комплекс алгоритмов и программного обеспечения для предварительной обработки сигналов с блока сцинтилляционных детекторов, формирования и обработки комптоновских изображений.
1. Комптоновская томография
Алгоритмы синтеза томограмм
Простейший метод визуализации - 3-координатное механическое сканирование рассеивающего объема (РО) по объекту контроля и регистрация в каждой точке интенсивности потока квантов в детекторе.
Применение линейки детекторов позволяет, осуществляя двухкоординатное сканирование, получить набор данных о каждом сечении объекта контроля. То есть, при прохождении сканера по одной координате У получаем двухмерный массив данных или массив сечений рассеяния. При сдвиге сканера по другой координате г и прохождении вновь по ¥ получаем следующее сечение объекта контроля и так далее.
Большую производительность формирования изображения можно получить в системах с двухко-ординатными позиционно-чувствительными детекторами-матрицами. В качестве элементов двухкоор-динатного пространственного отбора используются либо пинхол, либо микроканальные коллиматоры, применяемые гамма-камерами и осуществляющие параллельный перенос гамма-изображения на чувствительную поверхность детектора. Распределение интенсивности обратнорассеянного излучения на детекторе пропорционально распределению активности наведенных в объекте источников, а следовательно, и распределению плотности в реконструируемом слое.
Однако в силу ослабления излучения на изображение будет наложена помеха, определяемая распределением плотности в слоях, через которые проходит прямое и рассеянное излучение. Чем сложнее распределение плотности объекта, тем менее адекватна реконструкция.
В [1] описана система для выявления расслоений, трещин, инородных включений и других неод-нородностей, приводящих к изменению плотности материала при контроле слоистых конструкций из углепластиков. Источником излучения служит рентгеновская трубка на 250 кВ, а рассеянное излучение с помощью щелевого коллиматора переносится на позиционно-чувствительный детектор. Установка позволяет контролировать детали из углепластика толщиной до 70 мм при высокой чувствительности к дефектам.
Пример практически реализованной системы визуализации с пинхолом, построенной на таком принципе, приведен в [2]. Объект контроля засвечивается ленточным пучком, создаваемым щелевым коллиматором. Комптоновское изображение формируется на плоском детекторе (рентгеновская плёнка или двухкоординатный детектор) с помощью пинхольной апертуры. И в том и в другом случае формируемые изображения не свободны от помех, связанных с ослаблением и многократным рассеянием.
Одним из основных параметров комптоновской томографии, как и трансмиссионной, является пространственная разрешающая способность.
При разработке систем обычно рассматривают этот вопрос в приближении параллельных пучков первичного и рассеянного излучения [3-5], что не позволяет судить о реальной разрешающей способности. Более подробно вопрос о разрешающей способности комптоновских сканеров рассмотрен в [6] на примере анализа системы "ComScan".
Анализ геометрии с реальными, т.е. расходящимися пучками выполнен в НИИ интроскопии ТПУ при разработке системы "TomScan"[7, 8]. Рассмотрена геометрия контроля, представленная на рис.1. Отличие от обычно используемой модели состоит в том, что рассеивающий объем имеет форму не ромба, а более сложной фигуры. Его размер в направлении сканирования (DF) намного превышает диаметр зондирующего пучка. Поэтому переход через границу раздела двух сред сопровождается не скачкообразным переходом уровня сигнала S или его' логарифма с линии плотности среды 1 на линию
плотности среды 2, а происходит плавно, по линии с ограниченной крутизной. В связи с этим встаёт задача: по плавному переходу восстановить исходное положение границы раздела.
В [9] предложен алгоритм, основанный на проведении в точке максимальной крутизны логарифма и восстановлении линии раздела через эту точку. Предложенный алгоритм обеспечивает точность восстановления положения границы ±0,25 мм. Более детально этот вопрос изучен в [10]. Приведён пример сканирования тест-объекта (рис.2,а) рассеивающим объёмом (РО) по рис.1 с высотой РО DF= 8 мм. Тест-объект представляет собой многослойную структуру из слоёв органопластика разной плотности и воздуха. Суммарная толщина - 40 мм. Энергия первичного излучения - 200 кэВ. Массив значений сигнала с детектора, получаемый в процессе сканирования РО с шагом 0,1 мм, представлен графически на рис.2,б, а восстановленный массив плотности - на рис.2,е.
Проанализированы различные критерии для определения пространственного положения границы раздела двух сред. Для повышения точности предложен
8N .(х)
алгоритм, основанный на фиксации положения максимума функции £iV((x)+—^—.
Для определения значения плотностей двух сред, разделённых границей г, могут быть также использованы различные алгоритмы. Предложены два варианта. В [8] плотности находятся по точкам пересечения перпендикуляра с продолжениями линий плотностей. В работе используется алгоритм, основанный на решении интегрального уравнения, представляющего сигнал на детекторе Ns(x) в виде
интеграла свёртки распределения плотности в объекте р(х) с апертурной функцией рассеивающего объёма (АФРО) (функцией распределения чувствительности РО по оси X). В этом случае, имея массив значений Ns(x) и зная АФРО, можно восстановить исходное распределение плотности р(х). Погрешность восстановления р(х) в очень сильной степени зависит от того, насколько точно определена АФРО и насколько она деформируется (за счёт ослабления первичного и рассеянного излучения и
Рис. I. Геометрия контроля с реальными расходящимися пучками
многократного рассеяния) при движении РО в объекте вдоль оси X. В [11] приводится достигнутое разрешение по плотности ±10%.
При разработке комплекса программного обеспечения в настоящей работе в качестве отправного было принято требование о необходимости обеспечить пространственное разрешение порядка 0,1-0,3 мм при использовании первичного и рассеянного пучков с поперечными размерами ~1 мм. Эта задача может быть решена, если шаг перемещения РО составляет 0,1 мм. В этом случае при проходе РО по одной линии сканирования получается массив чисел (256 элементов), каждое из которых представляет собой количество квантов, рассеянных в детектор во всех точках РО. Здесь возможны две ситуации: размеры рассеивающего объёма малы по сравнению с размерами неоднородностей плотности. В этом случае в пределах всего РО плотность р можно считать постоянной. Элемент объема контроли-
1,5 1,0
0,5
Р(х) г/см-3
1мм
шшшш
//ЩУ/ЩЩщ
Ww/vv
Щ////У/ШЛ
25
1 щ
'А
1мм
ЛЛ
Рис.2. Тест-объект («); массив значений сигнала с детектора (шаг О, I мм) (Г>); массив плотности распределения (»)
руемого объекта, в пределах которого плотность практически не изменяется, называют элобом. То есть рассмотренная ситуация соответствует случаю, когда размеры РО не превышают размеры элоба. Для этого случая разработан рекуррентный алгоритм вычисления распределения плотности, представленный в [12].
Рассмотрим схему сканирования (рис.3). Узкий пучок гамма-излучения от точечного моноэнергетического источника S с энергией Е или от рентгеновской трубки с эффективной энергией Е (поз.1)
фокусируется коллиматором (поз.2) и нормально падает на контролируемый объект. Рассеянное в элементе объёма материала (поз.3) на угол 6, излучение, проходя через коллиматор (поз.5), попадает на детектор (поз.-/).
Коллиматор представляет собой цилиндрические каналы квадг ратного сечения с площадью поперечного сечения S, и 5,, соответственно и длинами каналов /?,. и hit.
Осуществляя пошаговое двухкоординатное сканирование с шагом по горизонтали (т.е. вдоль контролируемого сечения), равным апертуре коллиматора источника, и по вертикали с шагом, равным апертуре коллиматора детектора, деленной на sin8i(, мы условно разбиваем контролируемое сечение на элементы объёма (элобы) с усреднённой по элобу плотностью p/t; где ; - номер шага сканирования по вертикали или номер строки матрицы р,*; к - номер шага сканирования по горизонтали или номер столбца в матрице pik. Рассеяние квантов в этих злобах вызывает появление сигналов на детекторе.
Сканируя по горизонтали первый приповерхностный ряд элобов, для чего пошагово перемещают систему источник-детектор параллельно контролируемому сечению, получают на детекторе сигналы, однозначно связанные с р,* этих элобов. Альбедо-выход из рассеивающего объёма с координатами (2/с)
0 (1 - exp(-p|/,<7(,Up + V2|ip)))-^-N.,, (1)
mfPi
у/
vfffiWVW7W\
/ v /
// ?v
/1
Щ
/
Рис.3. Геометрия поточечного сканирования
« JVn da.
.0
До2
о " Р|/Дцр . »-ер/
где М\к - количество квантов, рассеянных в элобе р,4 и регистрируемых детектором в единицу времени; Цр,Цр - массовые коэффициенты ослабления излучения с энергиями Е0,ЕХ, Ех - энергия рассеянного пучка:
8. Заказ 867
1 + а(1 ~С08вх) '
(2)
Здесь а =
£р(МэВ) 0,511
Решение трансцендентного уравнения (1) относительно р|А осуществляется численным методом. Переместим систему источник-детектор на шаг по вертикали и снова осуществим горизонтальное сканирование. Сигнал, регистрируемый детектором при сканировании второго ряда элобов, связан с ослаблением в злобах второго ряда:
Ы-
2 к
а оЮ
Рг к
йр^Рц-и )))х
■ехр(-о(црр2А. 1-ехр(-д(ц°р1А. +^рЛр2к))
(3)
йрр2А- + Цр ^2р|А.+1 р|А. + л/2р2А. М
После подстановки вычисленных с помощью (1) плотностей р,* в (3) получают трансцендентное уравнение относительно р 1к, которое решается аналогично (1). Также восстанавливаются плотности всех нижележащих слоев. Сигнал, соответствующий рассеянию в элементе с индексами г, /с, имеет вид:
Р->. ......._ „
а
МрРд- +Мр>/2р,А-+,-_|
Г с<=/, + 1
-^- '-ехр(-«().1р Ерд.+л/2ц; X Р/,))-
ИрР«7с+Ир^2рЛ. /=2 /=2
С=к + 1~ 2
(4)
Решая уравнение (4) на каждом шаге сканирования по глубине, получают распределение плотности в поперечном сечении с пространственным разрешением, определяемым размерами пучков. Такой рекуррентный алгоритм автоматически корректирует затухание первичного и рассеянного излучений. Определённые трудности представляет учёт расходимости пучков и многократного рассеяния. При практической реализации алгоритма расходимость учитывается "в среднем". Для этого в выражении (7) в геометрический фактор Г вводится поправка g, равная отношению площади сечения реального РО к Г0, вычисленного для не расходящихся пучков при тех же размерах коллиматоров.
Недостатком отмеченного алгоритма является также влияние неровности поверхности, так как при низком качестве стартовых итераций погрешность вычисления с увеличением глубины накапливается. Для уменьшения этой погрешности на поверхности объекта устанавливают пластину из однородного материала с близкими атомным номером и плотностью и с толщиной, немного превышающей размер РО по глубине. Исходное положение перед сканированием - РО находится целиком в эталонной пластинке.
Данный алгоритм поточечного сканирования обеспечивает пространственное разрешение 2,4 мм по глубине и I мм в продольном направлении, поэтому его используют в режиме "БЫСТРЫЙ ПРОСМОТР" для обнаружения крупных локальных или протяжённых дефектов и различных закладок.
Для реконструкции распределения плотности с высоким пространственным разрешением используют алгоритм, основанный на сканировании. Шаг перемещения РО 0,1 мм как в поперечном, так и в продольном направлении. И здесь решается обратная задача - нахождение детального распределения плотности в сечении р (х, у) по известным апертурным функциям РО в поперечном и продольном направлениях. Рассмотрим рассеивающий объём, обозначенный на рис.1 как ИЕРС. Сканирование состоит в смещении блока источник-детектор (5,/)) относительно изделия толщиной Н таким образом, что РО проходит через переднюю границу, движется внутри изделия и выходит через нижнюю границу. Общий характер зависимости от координаты РО (х) интенсивности сигнала, снимаемого с детектора, приведён на рис.5, кривая Щх).
Рис.4. Геометрия линейного сканирования
Участкам перехода РО через границы соответствует "зона входа" и "зона выхода". Объект считается однородным, если его распределение плотности р(х)=соп51 (рис.6).
Главная цель обработки массива данных 7У(л) состоит в восстановлении по нему неизвестного распределения плотности р(х). Для РО конечных размеров характер изменения Щх) на участках перехода через границу неоднородностей (зона входа и выхода) определяется, кроме всего прочего, ослаблением первичного и рассеянного излучения внутри РО.
Понятие РО конечных размеров предполагает, что пренебречь этим ослаблением нельзя. Это обозначает, что апертурная функция РО (функция, описывающая изменение сигнала в детекторе при переходе через РО границы раздела) деформируется по мере его движения в слое. В связи с этим первым этапом обработки массива данных Щх), является коррекция ослабления в пределах РО. На участке движения РО внутри слоя ход функции Щх) достаточно точно для однородного материала описывается экспоненциальной функцией Щх). Амплитуда А определяется энергией первичного излучения Е0 и
1Ч(х)
ч,^К(х)=А(Е„р)е
"тх)
-1(Е„р)х
Рис.5. Зависимость /\'(л-), Ы'(х), Ф(л)
Рис.6. Взаимодействие апертуры РО и слоя
плотностью материала р, а постоянная экспоненты т - коэффициентами ослабления первичного цр и
рассеянного ц* потоков.
Переход АФРО через границу неоднородности осуществляется в направлении стрелки А (рис.6). Результат физического взаимодействия такой апертуры с объектом состоит в выходе квантов рассеянного излучения из области геометрического перекрытия апертуры и объекта. Удельный выход квантов из каждой точки пропорционален электронной плотности объекта в этой точке. Тогда общее число квантов, вышедшее из области перекрытия пропорционально интегралу по площади перекрытия объема БОРЕ и плотности р{х,у). То есть можно записать:
Ф'(х7 ) = 2 ¡к(хТ -х)р(х)с1х .
(5)
Здесь 1г(хт-х) будет представлять собой апертурную функцию РО (АФРО) или ядро интеграль ного уравнения. Из (5) видно, что Ф'(х) равно удвоенной свертке функций к{х), р(х):
О при хт < О, х^р при 0 <хтЕ, (й-х^Щу при Е<хт< Д О при хт > Б.
к{хт)-
(6)
При нормировке на Ф,пах =р тм\Ь{х)с^х получим
и{хт) = ] - х) РМ Лх = и(х)• Т"
•■'0 0 •'о
Выражение (7) представляет собой интегральное уравнение 1-го рода с интегралом типа свертки. Для решения такого интегрального уравнения могут быть использованы два алгоритма: решение полиномиальное (РП) и с помощью преобразования Фурье [13, 14].
Алгоритм полиномиального решения предназначен для восстановления по выборке экспериментальных данных и решения интегрального уравнения (10) в виде алгебраического полинома оптимальной степени. Он заключается в построении полиномов различных степеней, каждый из которых минимизирует функционал эмпирического риска для своей степени и выборке из них полинома, для которого риск принимает наименьшее значение. В алгоритме РП искомое приближение представляется линейной комбинацией полиномов Чебышева (¿¡{х)
р(х) = Ч(х) =
Ё^еДх).
./•=0
(8)
При фиксированной степени полинома к коэффициенты разложения а, определяются из условия минимизации величины функционала эмпирического риска
и,-1
О ./'=0
7*
I /=|
(9)
где с?,- - дисперсии замеров II,.
Оценка качества приближенного решения (т.е. величина функционала эмпирического риска, достигнутая на решении нормальной системы) производится по той же случайной выборке, по которой это приближение строилось. И для любой случайной выборки определяется выражением:
Л
М
(/с + 1)(1п/-1п(/с + 1)+1)-1пг|
(10)
где
- вероятность, с которой справедлива эта оценка; а* = (ф'ф) ]Ф'1/ - решение нормальной системы уравнений.
Изменяя степень полинома, которым описывается приближенное решение интегрального уравнения (10), получают различные значения оценки (10). Та степень, при которой оценка принимает наименьшее значение, является оптимальной степенью алгебраического полинома, приближающего решение интегрального уравнения. Сам этот полином определяется набором коэффициентов а, являющихся решением нормальной системы, которая соответствует оптимальной степени полинома.
Более оптимальным в смысле аппаратурной реализации является алгоритм с использованием преобразования Фурье [11,16]. Смысл этого алгоритма состоит в следующем.
Применяя преобразования Фурье к (7) и используя лемму Бореля о свертке, получим
Щ со) =/г(ю)р(ю). (П)
Отсюда р(ю) = . Применяя обратное преобразование Фурье, получим
Р(хг) = 1 т7техр(-уюхГУсо.
2л о п{ю)
(12)
Это решение неустойчиво, так как левая часть уравнения (7) известна приближенно, с погрешностью, определяемой в основном статистикой. В отмеченных работах время измерения всегда выбиралось таким, чтобы обеспечить статическую погрешность не хуже ±0,1% (106 квантов за время измерения). Так как и(хГ)= С/?(хг)+ У(хг), где У{хт) - помеха, то
Функция У(хТ) носит случайный характер, поэтому отношение ^^ может не иметь обратного
преобразования Фурье (ОПФ) из-за влияния высоких частот со. Но даже если и имеет ОПФ
п( со)
1У{хт), то уклонение функции ТГ{хт) от нуля может быть сколь угодно большим. Следовательно, в качестве приближенного решения (5) с приближенной левой частью нельзя брать точное решение этого уравнения. Точного решения может не существовать, а если оно и существует, то не обладает свойством устойчивости к малым отклонениям левой части С/(хг). Для повышения устойчивости надо
"подавить" влияние высоких частот помехи У(хт). Для этого можно функцию умножить на ре-
гуляризующий множитель /(со,а), зависящий от параметра а. Одним из типов регуляризирующего множителя может быть множитель типа
где ¿(со)=/г(со)*/г(-со)=|/2(со)|2, М(ю)=ю2Г, г - произвольное положительное число. Умножая решение на /(со,а), получим
/ ч и (а) ч и{ со) Л(о>) А(—ш)
/?(со) /г(со) |/7((о)
/г(со) 1г(- со)
|2 +асо2
или
р(со) =
|/г(со)
I2 +асо2
Применяя к (13) ОПФ, получим окончательное решение
Ц(со) /?(- со)
р(хГ) = —- /
2я о |/г(со) | +асо* Параметр а находится из критерия невязки, т.е. величины
ехр(-уйох-Г) ¿/со,
5 = ^рЫр(хг,а )~и(хт))гс1хт.
Параметр а, найденный из условия (17), считается оптимальным, и решение, вычисленное при этом а, является оптимальным решением.
Алгоритм решения задачи был разработан с учетом толщины сканируемого слоя 50 мм и шага сканирования Ах =0,1 мм, что дает сетку х( =/Ах, где 1-0,...,511. Левая часть (13) представлялась набором экспериментальных данных Щх) с отмеченной выше погрешностью. Интеграл при вычислении свертки в невязке представлялся разложением по формуле трапеции. Интегральные прямые и ОПФ были заменены их дискретными аналогами. Дискретные ПФ были реализованы БПФ, что существенно сократило время обработки данных. Восстановленная функция с помощью БПФ представлена на рис.2,е. Алгоритм реализован подпрограммой МАКБ5 [11]. Структурная схема алгоритма приведена на рис.7.
С целью оптимизации шага сканирования были получены массивы данных N для /=0,...,1023 (для шага Ах =0,05 мм) и для /=0,...,255 (для шага Ах =0,2 мм). Обработка этих данных показала, что при /=1023 существенного повышения пространственного разрешения и точности восстановления р(х) при заданной
погрешности результатов измерения N в 0,1% не происходит. С другой стороны время обработки и используемый подпрограммами объем памяти возрастают в среднем в четыре раза. Увеличение шага до 0,2 мм заметно ухудшает пространственное разрешение. Поэтому в качестве оптимального шага сканирования принят шаг Ах =0,1 мм и массив данных ./V из 512 элементов.
Алгоритмы РП и БПФ восстанавливают распределение плотности по /-й
(14)
(15)
(16)
(17)
Рис.7. Структурная схема алгоритма ПФ
линии сканирования в поперечном направлении. Полный набор экспериментальных данных содержит 256 таких массивов (номер массива - j) с шагом в продольном направлении также 0,1 мм.
После набора полного массива данных 256x256 элементов и восстановления плотности в каждом из v'-x поперечных массивов производится транспонирование с организацией / продольных массивов. Транспонирование осуществляется подпрограммой TRANS.
Затем с помощью тех же алгоритмов (РП или БПФ) производится восстановление распределения плотности в продольном направлении p(j>) с учетом АФРО в направлении оси Y. В результате такого двухкоординатного восстановления получается двухкоординатное распределение плотности p(x,j') в поперечном сечении объекта. Полученный двумерный массив фильтруется одной из имеющихся подпрограмм фильтрования и нормируется на диапазон целых чисел от 0 до 255. Этот массив представляет собой комптоновскую томограмму поперечного сечения. Томограмма выводится на экран монитора в черно-белой палитре с 256 полутонами.
2. Предварительная обработка информации со сцинтилляционных детекторов
2.1. Интерполяция зависимостей
Сигнал со сцинтилляционного детектора в виде импульсов поступает на счетчик. Идет счет импульсов за время Т, затем результат преобразуется в двоичный код и поступает в память компьютера.
В режиме "БЫСТРОГО ПРОСМОТРА" получаем значения количества рассеянных квантов за время измерения Тв позициях X объекта контроля, отличающихся на заданный интервал, содержащий несколько элементарных шагов по 0,1 мм (например 4 шага). Получаем массив Щх) при 7=сош1, Х=(Х0, X,. .... Х„_и Х„), где Х,ГХ,^=4. Т.е. в режиме "БЫСТРОГО ПРОСМОТРА" нам известен только каждый четвертый элемент массива. Остальные значения нужно восстановить. Рассмотрим кратко основные методы интерполяции.
Канонический полином: пусть функция Ы(х) задана таблицей значений, полученной из эксперимента или путем вычисления в последовательности значении аргумента х0, хь..., л"„. Выбранные значения аргумента х называются узлами таблицы.
В качестве аппроксимирующей функции ф(х) полиномом Р„(х) степени п в каноническом виде выбирается ср(а)=/,„(л')=с0+С|Х+г2а'2+...+с),х". Свободными параметрами интерполяции г, являются коэффициенты полинома. Интерполяция полиномами обладает такими преимуществами, как простота вычислений их значений, дифференцирования и интегрирования. Коэффициенты полинома г, определяются из условий Лагранжа Рп(х^=Ы-„ 0 < /' < п. Интерполяционный полином Лагранжа записывается в следующем виде: " " х-х,
= - ■ Практическое применение полино-
;=о /=о X,- — Ху ./*<
ма Лагранжа оправдано только в случае, когда интерполяционная функция вычисляется в сравнительно небольшом количестве точек х.
Полином для данной выше таблицы в виде, предложенном Ньютоном, выглядит так: Р„(х)=А0+А 1(х-хй)+А1(х-х0)(х-х])+...+А„(х-х0)(х-х,)...(а-х„_|). Коэффициенты поли-
нома определяются из условий Лагранжа Р„(х,)-Ыь 0 < г < п. Полагаем х=х0, тогда в формуле все слагаемые, кроме А0, обращаются в нуль, следовательно, А0 = Л^ [15]. Существуют также и другие разновидности интерполяций и аппроксимирующих функций [16].
Разновидности аппроксимирующих функций, таких, как параболическая регрессия, экспоненциальная и т.д., нам не подходят, так как статистическую природу сигнала с детекторов нельзя предсказать заранее. Поэтому в работе применялась линейная интерполяция.
С Конец
Рнс.8. Структурная схема интерполяции функции
В данной работе использован алгоритм полиномиального аппроксимирования с использованием полинома пятого порядка. Этот порядок оптимален, так как при уменьшении степени полинома мы не получаем точного восстановления функции, а при увеличении степени на единицу время, затрачиваемое машиной на интерполяцию, увеличивается в два раза.
На рис.8, приведена блок-схема алгоритма интерполяции.
Программа сначала запрашивает имя файла, в котором хранятся 100 значений табличной функции. Значения считываются в массив А1[1]. Затем запрос имени файла (Fl) для записи восстановленной функции из массива А2[1]. Далее берутся четыре элемента из массива AI [I] (kl:=l; nl:=4; for i:=kl to nl do). Для них вычисляются коэффициенты полинома В0 и В, (xx:=kl;). Потом идет запись пяти элементов (В0+В|*хх; хх:=хх+0.2; и так 5 раз) в массив А2[1]. Значения kl и nl увеличиваются на единицу (kl:=kl + l; nl:=nl + l;). Если n 1 < 100, то процесс повторяется. Если nl>100, то файлы чтения и записи закрываются. Выход. В результате получается массив из 512 значений, каждое из которых соответствует сигналу в стандартных позициях сканера с шагом 0,1 мм (N).
2.2. Нормировка
После восстановления функции N,{x) возникает проблема визуального представления графика функции. Так как разрешающая способность монитора компьютера ограничена (по вертикали 480 пикселов), то весь диапазон полученных значений N(x) отобразить одновременно невозможно. Кроме того, яркость свечения пикселов в результирующем изображении (или двухмерный оцифрованный массив чисел типа Byte) изменяется в пределах 0..255, т.е. 256 градаций серого цвета. Поэтому нормировку функции N(x) (или одномерного массива) необходимо произвести на 255. Т. е. диапазон значений нормированной функции будет от 0 до 255. На рис.9 представлена блок-схема алгоритма нормировки. Массив данных - z[I], где 1=1. .512.
Программа считывает данные из файла в массив. Затем находится максимальный элемент массива. Число 255 делится на этот элемент, получаем коэффициент. Далее умножаем каждый элемент массива на коэффициент, и полученные данные записываем в выходной файл.
Операция нормирования реализуется как отдельной программой NORM для формирования нормированного файла данных, так и подпрограммой NORMIROVKA в программе вывода графиков на экран компьютера GRAF.
2.3. Фильтрация
При сканировании функция N(x) содержит искажения, обусловленные статистической природой потока квантов и шумами. Поэтому полученные данные являются смесью сигнала с шумом и кроме того имеются искажения, вносимые интерполяцией.
Рис.9. Блок-схема алгоритма нормировки
В обоих случаях есть необходимость сглаживания получаемой функции или другими словами фильтрация. Фильтрацию осуществляет подпрограмма БИЛИ, в основу которой заложен алгоритм линейной фильтрации с переменным набором точек.
Структурная схема алгоритма программы фильтрования приведена на рис.10. Программа позволяет задать количество циклов фильтрования и число точек фильтрации, что позволяет выявить оптимальное сглаживание функции. Входными параметрами программы также являются имя файла для фильтрации - исходный файл; имя отфильтрованного файла - выходной файл.
Рис.10. Структурная схема алгоритма программы фильтрования
Пусть мы задали один цикл фильтрации по пяти точкам. Тогда программа открывает исходный файл и считывает данные в массив а 1 [I]. Затем берутся первые пять элементов из а1[1] (1:=1; а 1 [I] до а1[1+4]) и находится среднее для этих значений. Среднее значение присваивается 1+2 элементу в массиве а2[1+2]. Далее переменная I увеличивается на единицу и операция повторяется пока не достигнут конец массива данных. После этого массив значений а2[1] записывается в выходной файл.
Алгоритмы вычитания подставки или коррекции фона в массивах данных при отсутствии объекта контроля в работе не приведены, поскольку они просты и реализованы стандартными алгоритмами.
2.4. Визуализация процессов предварительной обработки
Для визуализации представленных операций и их результатов разработана программа GRAF. Структурная блок-схема программы представлена на рис.11. Программа GRAF позволяет одновременно выводить на дисплей до трех графиков функций распределения плотностей. Каждый график имеет свой цвет для лучшего восприятия и более точной оценки. Вывод трех графиков одновременно
Рис.1 I. Структурная блок-схема алгоритма программы GRAF
обусловлен тем, что необходимо обеспечить возможность сравнения полученных на этапах обработки функций. В программе GRAF используется разрешение монитора 640x480, что позволяет отобразить нужный диапазон интенсивностей, равный 256, и количество элементов в одномерном массиве, равное 512. Программа просмотра графиков также оснащена возможностью определения номера элемента в массиве данных и соответствующей ему интенсивности (нажатие клавиши "пробел"). Выход из программы - по нажатию "Enter".
3. Обработка двумерных массивов
В результате предварительной обработки получаем 2-мерный массив 256x256 элементов. Он уже, представляет собой оцифрованное изображение поперечного сечения в 256 градациях серого цвета, рис.12 (контрольный образец: углепластик р-1,3б г/см\ толщина 1 мм, 2 отверстия 1x1 мм, 2 отверстия 2x2мм, А=0,1 мм). В томографии для анализа изображений, определения различных характеристик дефектов и др. операций дальнейшие преобразования проводят над изображением. В настоящей работе использованы основные операции над изображениями, применяемые в томографии [17].
По аналогии с операциями над одномерными массивами существуют операции над 2-мерными массивами (оцифрованным изображением).
3.1. Основные операции над изображениями, применяемые в томографии
На рис.13 показаны различные ядра свертки. Как видно из рисунка, большинство матриц ядер имеют размерность 3x3 и нечетное число рядов и колонок. Этот формат коэффициентов свертки в пределах ядра используется в промышленности в качестве стандарта. Большой размер ядра повышает гибкость процесса свертки.
К сожалению, некоторые детали простого расчета весовой суммы свертки усложняют его реализацию. Первая и самая основная трудность связана с кромками изображения. Т. е. ядро свертки выходит краем за границы изображения. Такое искажение происходит на верхней, правой, левой и на нижней границах изображения. Для устранения такого явления можно использовать несколько методов. Наиболее исчерпывающими решениями являются: 1) данные на краях изображения не учитываются или 2) данные изображения могут копироваться, для того чтобы "синтезировать" дополнительные данные границы.
Вторая трудность связана с динамичностью диапазона новых значений преобразуемых элементов. Для большинства ядер новые значения преобразуемых элементов находятся в пределах действительных значений от 0 до 255. К сожалению, ядро размытия дает значения, выпадающие из диапазона действительных значений элементов изображения. По этой причине в алгоритм свертки вводится масштабирование.
Низкочастотные пространственные фильтры оставляют низкочастотные компоненты изображения нетронутыми и ослабляют высокочастотные компоненты. На рис.13 показаны три различных низкочастотных ядра. Отметим, что сумма значений ядра для всех низкочастотных фильтров равна 1.
Низкочастотная фильтрация используется для улучшения резкости изображения. Если отфильтрованное изображение с помощью низкочастотного фильтра вычитается из исходного изображения, то в результате будет относительное увеличение содержания высокочастотных компонент без увеличения шума в изображении. Субъективно результирующее изображение кажется резче, чем исходное. Это может быть использовано для освещения областей изображения, которые затемнены.
Низкочастотные пространственные фильтры
1/9 1/9 1/9
1/9 1/9 1/9
1/9 1/9 1/9 LP1
1/10 1/10 1/10
1/10 1/10 1/5 1/10 1/10 1/10 LP2 £= 1
1/16 1/8 1/16
1/8 1/4 1/8
1/16 1/8 1/16 LP3
Высокочастотные пространственные фильтры
НР1
о
о
5 -1 НР2 £=1
1 -2 1
-2 5 -2
1 -2 1 НРЗ
Ядро размытия 1
1 1
Е=25
Рис. i 3. Различные ядра свертки
Высокочастотные фильтры выделяют высокочастотные компоненты изображения, оставляя содержание низкочастотных компонент нетронутым. При использовании высокочастотной фильтрации возможно усиление края изображения. На рис.13 показаны 3 ядра высокочастотных фильтров.
Усредненное фильтрование использует значение элементов, содержащихся в области примыкания, для определения нового значения рассматриваемого элемента. Однако новое значение элемента не вычисляется алгоритмически из элементов в области примыкания. Вместо этого оно располагает элементы в области примыкания в возрастающем порядке и отбирает средние значения как новые значения рассматриваемых элементов. Алгоритм усредняющего фильтра показан на рис.14.
Результатом усредненного фильтрования является то, что любой случайный шум, содержащийся
Исходное изображение
Колонка, строка
Медианное или среднее значение
27 22 29
42 60j 29
27 25 29
Значение пикселов, отсортированные в возрастающем порядке
-> |22|25|27|27|29|29|29|42]б0"1
*
Значение 29 заменяет 60 в выходном Окружение точки 3x3 Если 60 было значением точки шума, создавшей
Обрабатывается центральная точка со значением 60
белое пятно на исходном изображении, пятно будет отфильтровано
Рис. 14. Работа усредняющего фильтра Рис. 15. Томограмма после фильтрования
в изображении, будет эффективно устранен. Это происходит потому, что любое случайное резкое изменение в интенсивности элемента в переделах области примыкания будет сортироваться, т. е. оно будет помещено либо на вершину, либо на нижнюю часть отсортированных значений области примыкания и не будет учитываться, так как для нового значения элемента всегда отбирается усредненное значение. Отбор случайных значений элементов показан на рис.14.
Результаты программы фильтрования представлены на рис.15.
4. Визуализация
4.1. Алгоритм и подпрограмма просмотра оцифрованного массива
данных на этапах обработки
На стадии обработки изображения возникает проблема визуализации промежуточных оцифрованных изображений. Для решения этой проблемы был разработан алгоритм и подпрограмма просмотра массивов данных.
Данные массива имеют тип "Byte" - это символьный тип. Значением символьного типа является множество всех символов ПК. Каждому символу присваивается целое число в диапазоне 0...255. Это число служит кодом внутреннего представления символа. Конечная цель работы - создание изображения в стандартном графическом формате, BMP, PCX и т. д. Т. е. в ASCII коде, а это и есть код внутреннего представления символа. Итак, конечное изображение в стандартном графическом формате обрабатывается и визуализируется любым графическим редактором.
Промежуточные изображения являются цифровыми. Модели компьютеров фирмы IBM позволяют выдать на монитор стандартными путями в системе DOS всего 64 градации серого цвета, а необходимо 256. Отсюда видно, что на каждую стандартную градацию серого приходится четыре значения нашего массива данных. Т. е. каждое четвертое наше значение имеет серый оттенок, а все остальные имеют значения стандартной цветовой палитры. Значит стоит задача перепрограммирования цветной стандартной палитры в серую.
В данной работе использовалась видеоплата Trident VGA, видеопамять 500 Kb. Загружаемый драйвер для работы монитора в режиме 640x480x256 - 8514ai.exe.
Для видеоплат IBM 8514 и VGA в языке программирования Turbo Pascal процедура SetRGBPalette(ColorNum, Red, Green, Blue: integer) позволяет модифицировать палитру, где
- ColorNuin - номер цвета (от нуля до 255);
- Red, Green, Blue - интенсивности лучей, соответственно красный, зеленый, синий.
Генерация полутоновой палитры, которая будет восприниматься как равнораспределенная, усложняется двумя обстоятельствами: нелинейностью мониторов для компьютера и нелинейностью системы зрительного восприятия [28].
На основе таблицы значений полутоновой палитры для 64 уровней градаций серого был разработан алгоритм перегрузки цветовой палитры в полутоновую, на 256 уровней. Как говорилось выше, каждый четвертый уровень соответствует действительно серому цвету, а трем последующим уровням нужно задать значения, лежащие на пределе. Т. е. отличающиеся от серого, но незаметные человеческому глазу. Любой цвет - это смешивание трех лучей (Red, Green, Blue). Чтобы получить оттенок серого цвета, надо смешать RGB в одинаковых пропорциях, так и делалось в работе для каждого четвертого значения. Для других значений RGB смешивались с небольшим отклонением от равности пропорций.
С помощью многочисленных экспериментов, результаты которых оценивались визуально, был разработан следующий алгоритм перегрузки цветовой палитры в серую:
- i := 2; - переменной i присваивается значение 2;
- while i < 254 do begin - цикл while выполняется пока i меньше 254;
- SetRGBPalette( i-2, i-4, i, i ); - перегрузка с завалом влево луча Red;
- SetRGBPalette( i-1, i, i, i );-перегрузка без завала;
- SetRGBPalette( i, i+4, i, i ); - перегрузка с завалом вправо луча Red;
- SetRGBPalette( i+1, i, i+4, i+4 ); - перегрузка с завалом вправо лучей Green и Blue;
- i := i + 4; - увеличение счетчика i на 4;
- end; - конец цикла;
- SetRGBPalette(0, 0, 0, 0); - восстанавливаем черный цвет.
4.2. Алгоритм и подпрограмма преобразования оцифрованного изображения в графический формат "BMP"
Сами по себе оцифрованные изображения имеют ограниченную ценность. Они становятся полезными, когда находятся в такой форме, что ими могут манипулировать другие прикладные программы.
Для этого необходимо привести данные изображения в "стандартную форму", которую могут интерпретировать и использовать другие прикладные программы [18]. Если изображение может быть введено в рисующую программу или программу электронной сверстки, его можно затем ввести в презентации, просмотры слайдов, отчеты, публикации и т.д.
BMP-файлы имеют формат, представленный на рис.16.
Файл содержит битовое изображение, начинается со структуры BITMAPFILEHEADER. Эта структура описывает тип файла и его размер, а также смещение области битов изображения.
Сразу после структуры BITMAPFILEHEADER в файле следует структура BITMAPINFO, которая содержит описание изображения и таблицу цветов. Описание изображения (размеры изображения, метод компрессии, размер таблицы цветов и т. д.) находится в структуре BITMAPINFOHEADER. В некоторых случаях (не всегда) в файле может присутствовать таблица цветов (как массив структур RGBQUAD), присутствующих в изображении.
Биты изображения обычно располагаются сразу после таблицы цветов.
BITMAPFILEHEADER
BITMAPINFOHEADER
Массив RGBQUAD (таблица цветов)
Биты изображения
BITMAPINFO
Рис. 16. Формат BMP файла
4.3. Алгоритм преобразования в BMP ф
1.Запрос подпрограммой имени BMP файла (Имя.Ьтр).
2.Запись в BMP файл BITMAPINFO (см. выше рис. 16).
З.Запрос имени файла с оцифрованным изображением.
4.Перевод оцифрованного изображения в ASCII код и запись в BMP файл.
5.Выход из подпрограммы.
Преобразование оцифрованного изображения в графический формат BMP осуществляется с помощью подпрограммы 1078ВМР, блок схема которой приведена на рис.17.
Полученный BMP файл обрабатывается любым современным графическим редактором, модифицируется, выбирается область интереса и т. п. По комптоновскому изображению дается заключение о наличии включений в объект контроля.
5. Пакет программного обеспечения комптоновского томографа
"TomScan-200"
Пакет представляет из себя набор программ, использующий механизм всплывающего меню, реализованный в графическом режиме адаптера Trident VGA с разрешением 640x480x256.
Применение графических возможностей ПК с адаптером VGA и выше обусловило возможность создания пакета программ для формирования изображений, наблюдения изображений в чёрно-белой палитре (256 полутонов серого), проводить операции с изображением. При этом можно формировать видео-базу на винчестере, дискетах, магнитной ленте, цифровых дисках.
Пакет программного обеспечения предусматривает возможность:
- математического моделирования процесса сканирования объекта рассеивающим объёмом произвольной формы или линейными альбедными суммами с заданием извне параметров контролируемого изделия, источника излучения, коллимационной системы;
- управления механизмом сканирования с двухкоординатным перемещением объекта относительно детекторного блока с шагом 0,1 мм;
- предварительной обработки информации с детекторного блока;
- формирования массива 256x256 элементов с амплитудой каждого элемента А^,, изменяющейся
от 0 до 255 (меню "Обработка");
- перенормировки и транспонирования массивов;
- просмотра оцифрованных файлов на этапах обработки (меню "Просмотр");
- установки параметров системы пространственного отбора и рентгеновской установки (меню "Установка");
- установки параметров обработки оцифрованных изображений (подменю "Параметры Обработки");
- пуска системы сканирования (меню "Пуск");
1 - обращения к справочной службе (меню "Help").
** Меню "Установка" включает в себя следующие окна и операции:
"Сканер" - окно параметров сканера (начальные, конечные координаты сканирования; напряжение на трубке; шаг сканирования по X и У);
"Выход" - выход из оболочки в операционную систему, завершение работы.
** Меню "Пуск" содержит команду пуска сканирующего устройства.
** Меню "Просмотр" позволяет визуализировать оцифрованное изображение на этапах обработки, выбор файла осуществляется с любого доступного компьютеру устройства (подпрограмма VIW).
** Меню "Обработка" включает в себя следующие окна и операции:
"Параметры обработки" - установка параметров для операций нормирования, транспонирования, фильтрования, умножения (соответственно подпрограммы NORM, TRANS, FILTR, UMNOJ)\
"Нормировка" - осуществляет нормирование оцифрованного изображения на 255 (255 - количество градаций оттенков серого);
"Транспонирование" - обеспечивает процессы вращения оцифрованного изображения;
о р м а т
Рис.17. Блок-схема программы 1078ВМР
"Фильтрование" - сглаживание шумов, помех изображения;
"Умножение" - домножение каждого элемента двухмерного изображения на заданный коэффициент.
При выборе какого-либо из пунктов оболочки, внизу экрана высвечивается надпись, поясняющая действие выбранной операции.
Пакет программного обеспечения является коммуникабельным, легко усовершенствуется подключением измененных или дополнительных подпрограмм или драйверов, прост в обращении с ним.
Структурная схема оболочки программы MENU приведена на рис.18.
Программа MENU использует подпрограммы, которые описаны выше. В основу создания программы заложен принцип объектно-ориентированного программирования. При завершении выполнения любого из пунктов меню управление передается в главное меню.
Рис. 18. Структурная схема программы MENU
Все подпрограммы написаны на языке Turbo Pascal и сгенерированы в виде ЕХЕ-файлов. Работа над программой MENU еще не закончена, так как на данном этапе созданы новые подпрограммы -это преобразование оцифрованного изображения в BMP формат, хранение данных изображения в ASCII коде и ряд разновидностей пространственных фильтров. В перспективе планируется переход на платформу Windows, изменение языка программирования, усовершенствование алгоритма программы с учетом требований пользователей.
СПИСОК ЛИТЕРАТУРЫ
1. Cappellini V., Guzzardy R., LicitzaG. A new 3-D object recognition tecnique using Compton tomography / Presented of the conference an advances in image processing and pattern recognition,Piza,Italy, 10-12 Dec., 1985.
2. Parish R.W., Cason D . W . J .// Private Communication. NDT International,-1977,-P. 181.
3. Преображенский Н.Г., Толиина С.П., Филинов В.Н. О некоторых возможностях комптоновской томографии: Тезисы 2-го Всесоюзн. сим. по вычислительной томографии - 1985,- С. 133-134.
4. Преображенский Н.Г., Тол ни на С.П. Спектротомография в пространстве импульсов на основе компто-новского рассеяния // Сб.: Линейные и нелинейные задачи вычислительной томографии - Н.: ВЦ СО АН СССР, 1985,-С.132-141.
5. T о л п и и а С . П . Алгоритмы комптоновской томографии в дефектоскопии. УДК 620.179.15.
6. Holt R.S., Cooper M.J., Jacson D.F. X and Y ray imaging tecniques // Nucl.lnstr.and Melh. in Phys.Research.- 1984,- Vol.221.- P.98-104,
7. Капранов Б . И . , Ч а н и н Г . С . , Маклашевский В . Я . , Филинов В.Н. Цифровая обработка изображений в томографии на комптоиовском обратном рассеянии // 14-я Российская научно-техническая конференция.-М„ I996.-C.352.
8. Варга В . В . , Капранов Б . И . и др. Реконструкция изображений в 3-D томографии на комптоиовском обратном рассеянии // 14-я Российская научно-техническая конференция,- М., 1996,- С.333.
9. Das G . , Pad hi Н.С. Directional Compton profile studies in KC1 //J. Phys.C.Solid State Phys.- 1987,- Vol.20- P.5253-5260.
10. Bridge B. A theoretical feasibility study of the use of Compton backscatter gamma-ray tomography (CBGT) for undewater offshore NDT // Brit.J. N.D.T.- 1985,- Vol.27.- P.357-363.
11. Strecker H, Scatter imaging of aluminium eastings using an X-ray fan beam and pinhole camera // Subst.to Mat.Eval.-1987.
12. Капранов Б. И., Челяди н А. М . , Гор бань Ю. И., Шаверин В. А. Способ измерения распределения плотности. А. С. № 167999, 1991 г.
13. Макс Ж. Методы и техника обработки сигналов при физических измерениях. Т. I. Основные принципы и классические методы / Пер. с франц. под ред. Н.Г. Волкова - М.: Мир, 1983.- С. 18-29.
14. Гольденберг Л.М. и др. Цифровая обработка сигналов//Учебное пособие для высших учебных заведений. 2-е изд.- М.: Радио и связь, 1990,- С. 123-138.
15. Мудров А . Е . Численные методы для ПЭВМ на языках БЕЙСИК, ФОРТРАН и ПАСКАЛЬ,-Томск.: Раско, 1992. -270 с.
16. Бруевич А.Н., Евтянов С.И. Аппроксимация нелинейных характеристик и спектры при гармоническом воз-действии.- М.: Советское радио, 1965.- 343 с.
17. Линяли К. Практическая обработка изображений / Пер. с англ.- М.: Мир, 1996,- 510 с.
18. Фролов А.В.. Фролов Г.В. Библиотека системного прораммиста. Т. 14. Графический интерфейс в MS WINDOWS.- М.: ДИАЛОГ-МИФИ, 1994,- 288 с.
