AN APPROACH TO THE PROBLEM OF EDGE ENHANCEMENT ON MULTISPECTRAL IMAGES Ivan G. Kazantsev Institute of Computational Mathematics and Mathematical Geophysics SB RAS, Novosibirsk, Russia Abstract In this work, we investigate an approach for combination of the images in several spectral intervals, obtained as the result of image processing in remote sensing. The technique consists in the reduction of two-dimensional algorithms to a set of one-dimensional edge detection procedures. A new approach of the amplitude and phase data information stored in spectral channels is proposed. The approach investigated is applicable to processing of multispectral satellite imagery. Keywords: Multispectral images, edge enhancement ОБ ОДНОМ ПОДХОДЕ К ЗАДАЧЕ УЛУЧШЕНИЯ ГРАНИЦ НА МНОГОСПЕКТРАЛЬНЫХ ИЗОБРАЖЕНИЯХ Казанцев И. Г. Институт вычислительной математики и математической геофизики СО РАН, Новосибирск В работе рассматривается комплексирование изображений в нескольких спектральных диапазо- нах, получаемых при дистанционных исследованиях. Метод основан на сведении двумерных алгорит- мов к набору одномерных процедур выделения границ. Предлагается также новый подход к использо- ванию амплитудных и фазовых данных различных каналов. Исследуемые методы применимы в обра- ботке многоспектральной информации дистанционного спутникового зондирования земной поверхно- сти. Ключевые слова: многоспектральные изображения, улучшение границ. Введение. В работе рассматривается проблема использования многоспектральной ин- формации для улучшения контраста границ на примере задачи восстановления изображений по рентгеновским проекциям, регистрируемым детекторами высокого спектрального разре- шения. Предполагается, что границы присутствуют на многоспектральных визуальных дан- ных в виде перепадов яркостей, наблюдаемых одновременно для всех каналов в исследуемых областях изображения. Модель формирования изображений известна и представлена некото- рым интегральным преобразованием, обращение которого может быть получено итерацион- ными методами [1]. Для известного распределения линейного коэффициента ослабления f ( x, y, z ) среды, источника A и детектора B, рассмотрим модель формирования данных в виде преобразования Радона: 𝐵 𝑝(𝐴, 𝐵) = ∫ 𝑓(𝑥, 𝑦, 𝑧) 𝑑𝑙. 𝐴 Здесь dl – элемент длины на линии AB. Задача состоит в реконструкции f по данным прони- кающего излучения p (A, B), регистрируемым большим множеством пар «источник - детек- тор», организованных в проекции. Итерационные методы рассматривают задачу восстановле- ния изображений как дискретную линейную систему Ax=b, где проекционные данные b (дис- кретный аналог p) являются взвешенными суммами определенных элементов (пикселов) изоб- ражения, лексикографически представленного в виде вектора x. Для различных видов дистан- ционного зондирования определяется своя системная матрица как для проникающего, так и для отраженного, а также других типов взаимодействия излучения и исследуемого объекта. Метод. Одним из активно используемых на практике методов является алгоритм Кач- мажа, известный также как ART (Algebraic reconstruction techniques) [2], использующий по- строчную обработку проекционной матрицы и потому экономный. В этой работе рассматри- вается итерационный алгоритм ART, модифицированный на основе методов условной мини- мизации полной вариации (TV – total variation), с целью улучшения качества реконструкции, уменьшения артефактов и уточнения границ на реконструкции при малом числе направлений просвечивания. Предполагается, что объекты имеют кусочно-постоянные коэффициенты ослабления. Это проявляется в скачкообразном изменении значений восстанавливаемой функ- ции от объекта (фона) к объекту, что наблюдается во многих задачах томографии и обработки изображений. Предполагается также, что эти скачки имеют место в одних и тех же областях изображений для всех спектральных каналов (рис. 1). Пусть Am – вектор-строка (m = 1, . . . , M) матрицы A длиной N; b – вектор-столбец правой части длиной M; λ - параметр релаксации; nIter – число итераций ART. В псевдокодах оче- редная итерация метода ART имеет вид (𝑘) (𝑘+1) (𝑘) < 𝐴𝑚 , 𝑥𝐴𝑅𝑇 > 𝑇 𝑥𝐴𝑅𝑇 ← 𝑥𝐴𝑅𝑇 + 𝜆 𝐴 < 𝐴𝑚 , 𝐴𝑇𝑚 > 𝑚 98 Рис. 1. Дискретная модель многоспектральной томографии. Слева – геометрическая модель прохож- дения луча AB через цифровое изображение размером 5 × 8 пикселов. Пикселы, пересекаемые лучом, обозначены вдоль луча как отсчеты (1 2 3 4 5 6 7) и составляют дискретную ступенчатую последова- тельность. Справа – профили значений коэффициента ослабления вдоль луча AB для различных энер- гий проникающего излучения в каналах Ch1, Ch2, … . Модуль градиента двумерного изображения x вычисляется в пикселе (s, t) по формуле |∇𝑥𝑠,𝑡 | = √(𝑥𝑠,𝑡 − 𝑥𝑠−1,𝑡 )2 + (𝑥𝑠,𝑡 − 𝑥𝑠,𝑡−1 )2 Целевой функцией выбрана L1-норма градиента, называемая полной вариацией (total varia- tion – TV) изображения x: ||𝑥||𝑇𝑉 = ∑ |∇𝑥𝑠,𝑡 | = ∑ √(𝑥𝑠,𝑡 − 𝑥𝑠−1,𝑡 )2 + (𝑥𝑠,𝑡 − 𝑥𝑠,𝑡−1 )2 𝑠,𝑡 𝑠,𝑡 Задача условной оптимизации формулируется так: найти изображение x 𝑥 = 𝑎𝑟𝑔𝑚𝑖𝑛 ||𝑥||𝑇𝑉 с ограничениями: 𝐴𝑥 = 𝑏, 𝑥 ≥ 0 Полученный алгоритм называется ART-TV и в псевдокоде имеет вид 1: function ART-TV (A, b, M, N, λ, α, nIter, nTV ) 2: x (0) ←0 3: for k ← 1 to nIter do 4: x (k) ← x (k−1) 5: x ART ← x (k) 6: for m ← 1 to M do 7: xART ← xART +λ/ ATm ART 8: end for 9: x POS ← max(0, x ART ) 10: β ←|| x POS − x ART || 2 11: x TV ← x POS 12: for q ← 1 to nTV do TV 13: dx TV ←∇ x ||x||TV 14: dx TV ← dxTV / || dx TV || 2 15: x TV ← x TV − α · β · dx TV 16: end for 17: x (k) ← max(0, x TV ) 18: end for 19: return x(nIter) 20: end function 99 INPUT: матрица A размером 𝑀 × 𝑁; вектор-столбец правой части b длины M; λ - пара- метр релаксации ART; α - параметр минимизации TV; nIter - число итераций ART; nTV – число итераций цикла минимизации полной вариации. OUTPUT: восстановленное изображение x(nIter). Проведены вычислительные эксперименты с реконструкциями по 20 проекциям мето- дами ART и ART-TV (Рис. 2, Рис. 3 и Рис. 4). Рис. 2. Дискретная модель. Слева – фантом 128х128х128. Справа – изображение центрального горизонтального слоя. Рис. 3. Реконструкция фантома вычислена по 20 проекциям каждая с 384 детекторами, в результате пяти ART итераций и пяти итераций при минимизации TV. Слева – минимизация полной вариации вычислена на всем двумерном изображении. Нормализированная среднеквадратичная ошибка равна 0.07. Справа – профили центральных столбцов фантома и реконструкций. Рис. 4. Реконструкции фантома вычислена по 20 проекциям каждая с 384 детекторами, в результате пяти ART итераций и пяти итераций при минимизации TV. Слева – минимизация осуществлена на одномерных массивах пикселов, посещаемых лучом просвечивания. Нормализированная среднеквад- ратичная ошибка равна 0.05. Справа – сравниваются профили центральных столбцов фантома и реконструкций. Использовалась минимизация полной вариации в двух вариантах – по всему двумерному носителю изображения (известный подход) и по наборам пикселов, лежащих на лучах просве- чивания (предлагаемый в этой работе). Во втором подходе все формулы вычисления вариаций 100 модифицированы на одномерный случай. Идея состоит в снижении размерности полной вари- ации и управляемой сегментации границ и скачков на многоспектральных реконструкциях. На данный момент изображения обрабатываются раздельно для каждого канала. Литература по обнаружению скачков в одномерных сигналах пополняется постоянно и появляется много но- вых эффективных методов [3]. Если имеется K спектральных каналов регистрации, а текущий луч пересекает поле изображения в L пикселах, то можно исследовать матрицы одномерных профилей (Рис. 1) размером 𝐾 × 𝐿 для совместного анализа границ на многоспектральной ин- формации, поскольку границы областей объективно существуют одновременно в определен- ном пикселе на всех каналах. В задачах дистанционного зондирования данные уже представ- лены многоспектральными изображениями и проблема улучшения и выделения границ может рассматриваться для произвольного набора прямых линий, набрасываемых на поле снимка случайно или согласно выбранному правилу. Возможно обобщение данного подхода на про- извольные криволинейные траектории обхода областей изображений синхронно для всех спектральных каналов с последующим применением оптимизационных методов [4], [5]. Рассмотрим множество изображений 𝑓𝑘 , 𝑘 = 1, … , 𝐾 в K спектральных диапазонах. Вы- числим для каждого изображения его амплитуду и фазу: (𝐴𝑘 , 𝜑𝑘 ) . Они также представляют собой некоторые изображения. Составим новые изображения {𝐵𝑘1 ,𝑘2 = (𝐴𝑘1 , 𝜑𝑘2 )} для каждых 1 ≤ 𝑘1 , 𝑘2 ≤ 𝐾 с амплитудными 𝐴𝑘1 и фазовыми изображениями 𝜑𝑘2 соответственно. Тогда можно рассматривать большую блочную 𝐾 × 𝐾 матрицу, диагональные блоки которой явля- ются исходными изображениями, а блоки вне диагонали - перекрестно синтезированные изоб- ражения, с целью классификации и исследования взаимной информации в спектральных ка- налах. Это предмет исследований в будущем. Заключение. В статье применен хорошо адаптируемый к различным геометриям про- свечивания и теоретически точный алгебраический итерационный метод ART-TV, в котором на каждой итерации осуществляется минимизация полной вариации градиента реконструиру- емого изображения и его подобластей. Рассматриваемые подходы могут быть применены в задаче улучшения границ на многоспектральных данных дистанционного зондирования вве- дением одномерной обработки на линейных подобластях изображений. Работа выполнена при финансовой поддержке РФФИ (грант № 16-07-00066) и Про- граммы I.33П Президиума РАН (проект № 0315-2015-0012). ЛИТЕРАТУРА [1] Гонсалес Р., Вудс Р. Цифровая обработка изображений. М.: Техносфера, 2005. 1072 с. [2] 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. [3] Condat L. A Direct Algorithm for 1-D Total Variation Denoising // IEEE Signal Processing Letters. – 2013. – V. 20. – P. 1054–1057. [4] Rudin L.I, Osher S., Fatemi E. Nonlinear total variation based noise removal algorithms // Physica D. – 1992. – V. 60. – P. 259–268. [5] Chambolle A. An Algorithm for Total Variation Minimization and Applications // Journal of Mathe- matical Imaging and Vision. – 2004. – V. 20. – P. 89–97. 101