Сегментация спутникового снимка, часть 1
ussmoСегментация спутникового снимка, часть 1
Это первый пример из цикла скриптов,
описывающих типичные задачи цифровой обработки
изображений. В первой части мы познакомимся с
базовыми подходами и алгоритмами обработки изображений
для задачи сегментации спутникового снимка.
Сегментация - это сопоставление пикселей
изображения различным классам.
В самом простом случае можно отделять
объекты на условном "переднем плане"
(foreground) от объектов условного "фона" (background).
В этом случае результатом сегментации
может служить бинарная маска - матрица таких же
размеров, что и исходное изображение, где
элементы 1 соответствуют интересующим нас
объектам, а элементы 0 - фону.
Рассматриваемые темы:
- импорт и визуализация цифровых изображений
- линейная и нелинейная фильтрация
- морфологические операции
- разработка пользовательских функций
Загрузка необходимых библиотек, импорт данных
Следующую строчку кода стоит разкомментировать,
если чего-то из применяемого не хватает:
# Pkg.add(["Images", "ImageShow", "ImageContrastAdjustment", "ImageBinarization", "ImageMorphology"])
Подключение загруженных библиотек:
using Images, ImageShow, ImageContrastAdjustment, ImageBinarization, ImageMorphology
Загрузка и визуализация исходного
цветного изображения функцией load.
В качестве обрабатываемого изображения
выступает спутниковый снимок
небольшого участка карты города,
на котором присутствуют как
области застройки, так и лесополоса.
Наша задача - попробовать отделить
лес от города.
I = load("$(@__DIR__)/map_small.jpg") # загрузка исходного цветного изображения
Предобработка
Как правило, в задачи предобработки
входят анализ изображения для выделения
характерных признаков объектов (отличия их от фона),
изменения контраста и/или яркости,
фильтрация для удаления шумов или сглаживания текстуры,
и проч.
Рассмотрим, какие из алгоритмов предобработки
подойдут в нашем случае. Учитывая, что изображение
цветное, можно попробовать рассмотреть отдельные
каналы (красный, зелёный и синий). Для этого
воспользуемся функцией channelview:
h = size(I, 1); # высота изображения, пикс
CV = channelview(I); # разложение изображения на каналы
[ RGB.(CV[1,:,:], 0.0, 0.0) RGB.(0.0, CV[2,:,:], 0.0) RGB.(0.0, 0.0, CV[3,:,:]) ]
Для того, чтобы получить бинарную маску,
нам нужно работать с изображением в
градациях серого (так называемой матрицей интенсивности).
Мы можем перевести цветное изображение в
матрицу интенсивности функцией Gray, или же
рассмотреть один из каналов. Как мы видим,
в синем канале область леса на изображении
сильнее отличается по интенсивности от области
города. И задача сегментации сводится к отделению
более ярких областей пикселей от более тёмных.
[RGB.(CV[3,:,:]) simshow(ones(h,20)) Gray.(I)]
Теперь воспользуемся линейным сглаживающим фильтром,
чтобы выровнять интенсивность близко расположенных
пикселей в синем канале. Сравним два линейных
фильтра - фильтр Гаусса со
скользящим окном размером 9 на 9 пикселей, а также
простой осредняющий фильтр с окном размером 7 на 7.
Применим их к матрице BLUE при помощи функции
imfilter, которой мы передадим разные кернелы.
BLUE = CV[3,:,:]; # интенсивность синего канала RGB
gaussIMG = imfilter(BLUE, Kernel.gaussian(2)); # линейный фильтр Гаусса
kernel_average = ones(7,7) .* 1/(7^2); # кернел среднего арифметического
avgIMG = imfilter(BLUE, kernel_average); # линейный фильтр среднего
[simshow(gaussIMG) simshow(ones(h,20)) simshow(avgIMG)]
При выбранных параметрах используемых кернелов,
результат принципиально не отличается и нас вполне
устраивает. В качестве основного изображения для
дальнейших операций возьмём результат простого
осредняющего фильтра - матрицу avgIMG.
Бинаризация
Бинаризация - это процесс получения матрицы,
содержащий бинарные значения, то есть 1 или 0.
На изображении это соответствует белым и
чёрным пикселям. В сегментации белые пиксели
будут соответствовать
одному из классов, а черные - другому.
Сравним два метода бинаризации - метод Отсу,
доступный посредством функции binarize, и простое
сопоставление интенсивности пороговому значению.
binOTSU = binarize(avgIMG, Otsu()); # бинаризация методом Отсу
binTHRESH = avgIMG .> 0.25; # бинаризация по порогу 0.25
BW_AVG = binTHRESH;
[ simshow(binOTSU) simshow(zeros(h,20)) simshow(binTHRESH) ]
Результаты так же почти идентичны,
и мы продолжим с бинарным изображением после
простого сравнения с порогом.
Несмотря на то, что мы получили бинарную
матрицу, маску, на ней весьма много шума: области
имеют разрывы, "прорези", присутствует много
мелких объектов. Чтобы сделать маску более
"гладкой" и свободной от шума, мы воспользуемся
алгоритмами бинарной морфологии.
Морфологические операции
Морфологические операции можно осуществлять
как над бинарными изображениями, так и над
матрицами интенсивности. Но для понимания базовых
алгоритмов будем ограничиваться бинарной
морфологией.
В основных алгоритмах, таких как дилатация,
раскрытие, закрытие или эрозия, используются
специальные кернелы - структурные элементы.
Нам доступны параметризованные функции для создания
типичных структурных элементов типа "ромб" или "квадрат",
впрочем, нам интересно создать круглый структурный
элемент. Поэтому мы бинаризуем двумерное
распределение Гаусса по порогу в 0.0025.
kernel_d = strel_diamond((13, 13), r=6); # структурный элемент "ромб"
kernel_b = strel_box((13, 13), r=4); # структурный элемент "квадрат"
se = Kernel.gaussian(3) .> 0.0025; # структурный элемент "пользовательский"
[ simshow(kernel_d) simshow(zeros(13,13)) simshow(kernel_b) simshow(zeros(13,13)) simshow(se) ]
Теперь осуществим операцию морфологического закрытия
функцией closing
с пользовательским структурным элементом. Это позволит
"закрыть" небольшие отверстия и щели между белыми
пикселями, а также сгладить форму объектов (блобов).
Объектами мы считаем "объединённые" области белых
пикселей, окружённых черными пикселями, или границей
изображения.
После операции закрытия мы воспользуемся
функцией area_open и удалим объекты, размер
которых не превышает 500 пикселей.
closeBW = closing(BW_AVG,se); # морфологическое закрытие
openBW = area_opening(closeBW; min_area=500) .> 0; # удаление "малых" объектов
[ simshow(closeBW) simshow(zeros(h,20)) simshow(openBW) ]
Теперь инвертируем нашу бинарную маску,
и удалим объекты размером менее 5000 пикселей.
Инвертированную маску также сгладим операцией
закрытия:
fillBW = imfill(.!openBW, (0,5000)); # удаление "крупных" объектов
MASK_AVG = closing(fillBW,se); # морфологическое закрытие (сглаживание)
[ simshow(fillBW) simshow(ones(h,20)) simshow(MASK_AVG) ]
Результат сегментации
Мы хотим убедиться, что полученная
бинарная маска сегментирует исходное изображение.
Для этого изображения нужно совместить.
Для выбранной нами визуализации добавим по 30%
интенсивности бинарной маски в красный и
зелёный каналы исходного цветного изображения.
Таким образом, получим полупрозрачную
маску жёлтого цвета:
sv1_AVG = StackedView(CV[1,:,:] + (MASK_AVG./3), CV[2,:,:] +
(MASK_AVG./3), CV[3,:,:]); # наложение "жёлтой маски"
forest_AVG = colorview(RGB, sv1_AVG);
sv2_AVG = StackedView(CV[1,:,:] + (.!MASK_AVG./3), CV[2,:,:] +
(.!MASK_AVG./3), CV[3,:,:]); # наложение "жёлтой маски"
city_AVG = colorview(RGB, sv2_AVG);
[ forest_AVG simshow(ones(h,20)) city_AVG ]
Визуально можем оценить результат, как неплохой.
Мы наблюдаем успешное выделение лесного массива.
Алгоритм не является универсальным, так как
многие параметры
алгоритмов подбирались к похожим
тестовым изображениям. Впрочем, целью примера является
демонстрация процесса разработки прототипов
алгоритмов ЦОИ.
Нелинейная фильтрация
Рассмотрим ещё один подход к предобработке
изображения. В этом случае мы будем использовать
растяжение гистограммы изображения
в градациях серого, а также нелинейный
текстурный фильтр.
Для начала применим к монохромному изображению
функцию adjust_histogram с методом
LinearStretching. Мы получили очень
контрастное изображение с большим числом
ярких и тёмных пикселей.
IMG = adjust_histogram(Gray.(I), LinearStretching(dst_minval = -1, dst_maxval = 1.8));
[ Gray.(I) simshow(ones(h,20)) IMG ]
Теперь воспользуемся так называемым фильтром диапазона.
Это нелинейный фильтр, который внутри скользящего
кернела определяет разницу между максимальным и минимальным
значением интенсивности массива.
Для этого напишем простую пользовательскую функцию
range, в основе которой лежит встроенная функция
extrema, возвращающая два значения.
function range(array) # функция нелинейной фильтрации
(minval, maxval) = extrema(array); # мин. и макс. значения массива
return maxval - minval; # диапазон (макс. - мин.)
end
Применим пользовательскую функцию к скользящему
по изображению окну при помощи функции mapwindow.
Результат нелинейной фильтрации бинаризуем методом Отсу.
rangeIMG = mapwindow(range, IMG, (7,7)); # функция `range` в окне 7х7
BW_RNG = binarize(rangeIMG, Otsu()); # бинаризация методом Отсу
[ rangeIMG simshow(zeros(h,20)) BW_RNG ]
В области городской застройки
интенсивность пикселей меняется сильнее и
пространственно чаще, чем в области лесополосы.
Поэтому текстурные фильтры могут применяться для
разделения областей на изображениях, средняя интенсивность
которых слабо различается после стандартных
линейных фильтров.
Пользовательская функция морфологии
Соберём последующие за бинаризацией операции
в одну функцию myMorphology и применим её
к бинарному изображению:
function myMorphology(BW)
se = Kernel.gaussian(3) .> 0.0025;
closeBW = closing(BW,se);
openBW = area_opening(closeBW; min_area=500) .> 0;
fillBW = imfill(.!openBW, (0,5000));
MASK = closing(fillBW,se);
return MASK
end
Визуализируем результирующую маску:
MASK_RNG = myMorphology(BW_RNG);
sv1_RNG = StackedView(CV[1,:,:] + (MASK_RNG./3), CV[2,:,:] +
(MASK_RNG./3), CV[3,:,:]);
forest_RNG = colorview(RGB, sv1_RNG);
sv2_RNG = StackedView(CV[1,:,:] + (.!MASK_RNG./3), CV[2,:,:] +
(.!MASK_RNG./3), CV[3,:,:]);
city_RNG = colorview(RGB, sv2_RNG);
[ forest_RNG simshow(ones(h,20)) city_RNG ]
Попробуйте сравнить результаты обработки при линейной
и нелинейной фильтрации программным образом.
Пользовательская функция сегментации
Наконец, соберём все последовательные
алгоритмы предобработки, нелинейной фильтрации
и морфологии в ещё одну функцию mySegRange,
принимающей на вход цветное изображение:
function mySegRange(I)
CV = channelview(I);
IMG = adjust_histogram(Gray.(I), LinearStretching(dst_minval = -1, dst_maxval = 1.8));
rangeIMG = mapwindow(range, IMG, (7,7));
BW_RNG = binarize(rangeIMG, Otsu());
MASK_RNG = myMorphology(BW_RNG);
return MASK_RNG
end
И применим её к более крупному снимку местности:
I2 = load("$(@__DIR__)/map_big.jpg");
CV2 = channelview(I2);
myMASK = mySegRange(I2);
sview = StackedView(CV2[1,:,:] + (.!myMASK./3), CV2[2,:,:] + (.!myMASK./3), CV2[3,:,:]);
colorview(RGB, sview)
Практический смысл подобных операций - это возможность
осуществлять анализ и вычисления над выделенными
объектами в автоматическом режиме.
Зная линейные размеры пикселя на карте в метрах, можно,
к примеру,
определить площадь объектов на карте.
В нашем случае мы можем подсчитать площадь (в км^2),
занимаемую лесополосой, городской застройкой, а также
процентное отношение леса к общей площади
рассматриваемого участка:
pixsquare = (800/74)^2; # площадь одного пикселя растра карты (км^2)
forest_area_km2 = round((sum(myMASK) * pixsquare) / 1e6, sigdigits = 4) # площадь леса
city_area_km2 = round((sum(.!myMASK) * pixsquare) / 1e6, sigdigits = 4) # площадь города
forest_percentage = round(100*sum(myMASK)/length(I2), sigdigits = 4) # процент леса
Заключение
На примере обработки спутникового снимка мы познакомились
с различными методами предобработки изображений,
бинаризации и морфологических операций для задачи
сегментации, а также узнали, как быстро разрабатывать
пользовательские функции цифровой обработки изображений.