CC BY
28
УДК 519.652
КОМПЛЕКСИРОВАНИЕ МНОГОСПЕКТРАЛЬНОЙ ИНФОРМАЦИИ В ЗАДАЧАХ УТОЧНЕНИЯ ГРАНИЦ НА ИЗОБРАЖЕНИЯХ В ТОМОГРАФИИ
Иван Гаврилович Казанцев
Институт вычислительной математики и математической геофизики СО РАН, 630090, Россия, г. Новосибирск, пр. Академика Лаврентьева, 6, доктор физико-математических наук, доцент, старший научный сотрудник лаборатории обработки изображений, тел. (383)330-73-32, e-mail: kig@ooi.sscc.ru
В работе предложены подходы комплексирования изображений в нескольких спектральных диапазонах, получаемых при реконструкции изображений в томографии. Метод основан на сведении двумерной полной вариации, многократно вычисляемой в итерациях оптимизационных алгоритмов томографии, к набору одномерных процедур выделения границ. Предлагается также новый подход к использованию амплитудных и фазовых данных различных каналов. Исследуемые методы применимы в обработке многоспектральной информации дистанционного спутникового зондирования земной поверхности.
Ключевые слова: многоспектральные изображения, томография.
COMBINATION OF MULTISPECTRAL INFORMATION IN THE PROBLEM OF DIGITAL IMAGE EDGES ENHANCEMENT IN TOMOGRAPHY
Ivan G. Kazantsev
Institute of Computational Mathematics and Mathematical Geophysics SB RAS, 630090, Russia, Novosibirsk, 6 Akademik Lavrentiev Prospect, D. Sc., Senior scientific researcher of Image Processing Laboratory, tel. (383)330-73-32, e-mail: kig@ooi.sscc.ru
In this work, an approach for combination of the images in several spectral intervals, obtained as the result of image reconstruction in tomography is suggested. The technique is based on reduction of two-dimensional total variation, multiply used in iterations of optimization techniques of tomography, to a set of one-dimensional edge detection procedures. A new approach of the amplitude and phase data stored in spectral channels is proposed. The approach investigated can be applied to processing of multispectral satellite imagery.
Key words: multispectral images, tomography.
В работе рассматривается задача восстановления изображений по рентгеновским проекциям, регистрируемым детекторами высокого спектрального разрешения. Идеализированная интегральная модель рентгеновской компьютерной томографии хорошо известна. Для известного распределения линейного коэффициента ослабления f (x, y, z) среды, источника A и детектора B малых размеров, модель формирования данных имеет вид лучевого преобразования Радона:
р(А,В) = SBAf(x,y,z) dl.
Здесь dl - элемент длины на линии AB. Задача состоит в реконструкции f по данным p (A, B), регистрируемым большим множеством пар «источник - детектор», организованных в проекции. Итерационные методы рассматривают задачу восстановления изображений как дискретную линейную систему Ax = b, где проекционные данные b (дискретный аналог p) являются взвешенными суммами определенных элементов (пикселов) изображения, лексикографически представленного в виде вектора x. Одним из активно используемых на практике методов является алгоритм Качмажа, известный также как ART (Algebraic reconstruction techniques) [1], использующий построчную обработку проекционной матрицы и потому экономный. В этой работе рассматривается итерационный алгоритм ART, модифицированный на основе методов условной минимизации полной вариации (TV - total variation), с целью улучшения качества реконструкции, уменьшения артефактов и уточнения границ на реконструкции при малом числе направлений просвечивания. Предполагается, что объекты имеют кусочно-постоянные коэффициенты ослабления. Это проявляется в скачкообразном изменении значений восстанавливаемой функции от объекта (фона) к объекту, что наблюдается во многих задачах томографии и обработки изображений. Предполагается также, что эти скачки имеют место в одних и тех же областях изображений для всех спектральных каналов (рис. 1).
а) б)
Рис. 1. Дискретная модель многоспектральной томографии:
а) геометрическая модель прохождения луча AB через цифровое изображение размером 5x8 пикселов. Пикселы, пересекаемые лучом, обозначены вдоль луча как отсчеты (1 2 3 4 5 6 7) и составляют дискретную ступенчатую последовательность; б) профили значений линейного коэффициента ослабления вдоль луча AB для различных энергий проникающего излучения в каналах Ch 1, ...
Пусть Am - вектор-строка (m = 1, . . . , M) матрицы A длиной N; b - вектор-столбец правой части длиной M; X - параметр релаксации; nIter - число итераций ART. В псевдокодах очередная итерация метода ART имеет вид
(fc + 1) „СЮ ■ 1 <ArrvXART>AT
^ Л DT* ^ ^ Л DT* + ^ . .Т Лгг
ART ART ' 'L <л АТ >
m-
Модуль градиента двумерного изображения x вычисляется в пикселе (s, t) по формуле
|Vxs,t| = V (xs,t — xs-1,t)2 + (xs,t — xs,t-l)2 •
Целевой функцией выбрана Li - норма градиента, называемая полной вариацией (total variation - TV) изображения x:
|х|
I TV = ^ |Vxs,t| = ^ J(xs,t — xs-iit)2 + (xSit — xSit-i) s,t s,t
Задача условной оптимизации формулируется так: найти изображение х
х = агдтт ||х||ту
с ограничениями:
Ах = Ь,х > 0
Полученный алгоритм называется ART-TV и в псевдокоде имеет вид
1 2
3
4
5
6
7
8
9
10 11 12
13
14
15
16
17
18
19
20
function ART-TV (A, b, M, N, 1, a, niter, nTV) * 0 ^0
for k ^ 1 to niter do
* (k) ^ x (k-1)
* ART x (k)
for m ^ 1 to M do
*ART ^ *ART +X<Am, *art>/<A m » ^ m^ n
end for
* POS ^ max(0, x ART )
V <
xPOS - x ART II 2
x TV — x POS
for q — 1 to nTV do
• TV - Fx \\X\\tv
dxTV —i
dx TV — dxTV / \\ dx TV \
x TV — x TV - a ■ в ■ dx
TV
TV
end for
x
(k)
——
max(0, x TV)
end for return x(nIter) end function
ART
TV
INPUT: матрица A размером MxW; вектор-столбец правой части b длины M; X - параметр релаксации ART; a - параметр минимизации TV; niter - число итераций ART; nTV - число итераций цикла минимизации полной вариации. OUTPUT: восстановленное изображение x(nIter\
2
2
а)
б)
Рис. 2. Дискретная модель: а) фантом 128 х 128 х 128; б) изображение центрального горизонтального слоя
в) г)
Рис. 3. Две реконструкции фантома вычислены по 20 проекциям каждая с 384 детекторами, в результате пяти ART итераций и пяти итераций
при минимизации TV:
а) минимизация полной вариации, вычисленной на всем двумерном изображении; нормализированная среднеквадратичная ошибка равна 0.07; б) профили центральных столбцов фантома и реконструкций; в) минимизация осуществлена на одномерных массивах пикселов, посещаемых лучом просвечивания; нормализированная среднеквадратичная ошибка равна 0.05; г) сравниваются профили центральных столбцов фантома и реконструкций
Проведены вычислительные эксперименты с реконструкциями по 20 проекциям методами ART и ART-TV (рис. 2 и 3). Использовалась минимизация полной вариации в двух вариантах - по всему двумерному носителю изображения (известный подход) и по наборам пикселов, лежащих на лучах просвечивания (предлагаемый в этой работе). Во втором подходе все формулы вычисления вариаций модифицированы на одномерный случай. Идея состоит в снижении размерности полной вариации и управляемой сегментации границ и скачков на многоспектральных реконструкциях. На данный момент изображения обрабатываются раздельно для каждого канала. Литература по обнаружению скачков в одномерных сигналах пополняется постоянно и появляется много новых эффективных методов [2]. Если имеется K спектральных каналов регистрации, а текущий луч пересекает поле изображения в L пикселах, то можно исследовать матрицы одномерных профилей (см. рис. 1) размером К X L для совместного анализа границ на многоспектральной информации, поскольку границы областей объективно существуют одновременно в определенном пикселе на всех каналах.
Рассмотрим множество изображений fk,k = l,...,K в K спектральных диапазонах. Вычислим для каждого изображения его амплитуду и фазу: (Ak,^k). Они также представляют собой некоторые изображения. Составим новые изображения {Вк к = (Aki, ^к2)} для каждых 1 < kt,k2 < К с амплитудными Aki и фазовыми изображениями фк2 соответственно. Тогда можно рассматривать большую блочную К X К матрицу перекрестно синтезированных изображений с целью классификации и исследования взаимной информации в спектральных каналах. Это предмет исследований в будущем.
В статье применен хорошо адаптируемый к различным геометриям просвечивания и теоретически точный алгебраический итерационный метод ART-TV, в котором на каждой итерации осуществляется минимизация полной вариации градиента реконструируемого изображения и его подобластей. Рассматриваемые подходы могут быть применены в задаче улучшения границ на многоспектральных данных дистанционного зондирования введением одномерной обработки на линейных подобластях изображений.
Работа выполнена частично при финансовой поддержке Российского фонда фундаментальных исследований (проект № 16-07-00066) и Программы 1.33ППрезидиума РАН (проект № 0315-2015-0012).
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Sidky E. Y., Pan X. Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization // Physics in Medicine and Biology. - 2008. - V. 53. -P.4777-4807.
2. Condat L. A Direct Algorithm for 1-D Total Variation Denoising // IEEE Signal Processing Letters. - 2013. - V. 20. - P. 1054-1057.
© И. Г. Казанцев, 2017
