Сегментация спутникового снимка, часть 1
Сегментация спутникового снимка, часть 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) # процент леса
Заключение¶
На примере обработки спутникового снимка мы познакомились с различными методами предобработки изображений, бинаризации и морфологических операций для задачи сегментации, а также узнали, как быстро разрабатывать пользовательские функции цифровой обработки изображений.