Документация Engee
Notebook

Сегментация спутникового снимка, часть 2

Второй пример из цикла по цифровой обработке изображений. В первой части мы рассматривали задачу сегментации спутникового снимка и базовые подходы к её решению. Во второй части мы рассмотрим дополнительные методы, доступные начинающему разработчику, в частности:

  • цветовая сегментация

  • фильтр Собеля

  • SRG-алгоритм

Загрузка необходимых библиотек, импорт изображения {#Загрузка-необходимых-библиотек,-импорт-данных}

Следующую строчку кода стоит разкомментировать, если чего-то из применяемого не хватает:

In [ ]: 
# Pkg.add(["Images", "ImageShow", "ImageContrastAdjustment", "ImageBinarization", "ImageMorphology", "ImageFiltering", "ImageSegmentation"])

Подключаем необходимые библиотеки:

In [ ]: 
using Images, ImageShow, ImageContrastAdjustment, ImageBinarization, ImageMorphology, ImageFiltering, ImageSegmentation

Загрузка исходного цветного изображения (спутникового снимка):

In [ ]: 
I = load("$(@__DIR__)/map_small.jpg")
Out[0]:
No description has been provided for this image

Цветовая сегментация

В первой части мы практически не обращались к цветовой информации, но очевидно, что лесной массив в основном зелёный, в то время как город содержит много серого и белого. А значит, мы можем попробовать выделить области, близкие по цвету к наблюдаемому на изображении оттенку зелёного. Мы будем бинаризовать отдельные каналы изображения, а доступ к ним мы получим, как и ранее, при помощи функции channelview:

In [ ]: 
(h, w) = size(I);
CV = channelview(I);
[ RGB.(CV[1,:,:], 0.0, 0.0) RGB.(0.0, CV[2,:,:], 0.0) RGB.(0.0, 0.0, CV[3,:,:]) ]
Out[0]:
No description has been provided for this image

Типичный "зелёный пиксель" находится, как можно предположить, в центре изображения. Выделим "центральный" пиксель, зная ширину и высоту изображения, и посмотрим на значения его интенсивности в цветовых каналах:

In [ ]: 
midh = Int(round(h/2));
midw = Int(round(w/2));
midpixel = CV[:,midh, midw]
Out[0]:
3-element Array{N0f8,1} with eltype N0f8:
 0.22N0f8
 0.267N0f8
 0.118N0f8

Теперь, используя стандартные логические операции сравнения, бинаризуем каждый из каналов. В красном канале логической единице будут соответствовать пиксели с интенсивностью от 0.15 до 0.25, в зелёном - от 0.2 до 0.3, а в синем - от 0.05 до 0.15.

In [ ]: 
BIN_RED = (CV[1,:,:] .> 0.15) .& (CV[1,:,:] .< 0.25) ;
BIN_GREEN = (CV[2,:,:] .> 0.2) .& (CV[2,:,:] .< 0.3);
BIN_BLUE = (CV[3,:,:] .> 0.05) .& (CV[3,:,:] .< 0.15);
simshow([BIN_RED BIN_GREEN BIN_BLUE])
Out[0]:
No description has been provided for this image

А теперь объединим три бинарные маски в одну операцией логического "И" - теперь белый пиксель будет только там, где в исходном цветном изображении все три проверяемых условия (диапазона) выполняются:

In [ ]: 
BIN = BIN_RED .& BIN_GREEN .& BIN_BLUE;
simshow(BIN)
Out[0]:
No description has been provided for this image

Сделать из этого пёстрого набора пикселей маску нам вновь поможет морфология. В этот раз мы возьмём небольшой структурный элемент (7х7 пикселей) в форме ромба:

In [ ]: 
se = strel_diamond((7,7))
Out[0]:
7×7 ImageMorphology.StructuringElements.SEDiamondArray{2, 2, UnitRange{Int64}, 0} with indices -3:3×-3:3:
 0  0  0  1  0  0  0
 0  0  1  1  1  0  0
 0  1  1  1  1  1  0
 1  1  1  1  1  1  1
 0  1  1  1  1  1  0
 0  0  1  1  1  0  0
 0  0  0  1  0  0  0

И оценим результат операции морфологического закрытия:

In [ ]: 
closeBW = closing(BIN,se);
simshow(closeBW)
Out[0]:
No description has been provided for this image

Результирующий вид маски мы получим, удалив небольшие "блобы", а также осуществив операцию закрытия ещё раз, теперь уже для сглаживания границ основных "крупных" областей объединённых белых пикселей:

In [ ]: 
openBW = area_opening(closeBW; min_area=500) .> 0;
se2 = Kernel.gaussian(3) .> 0.0025;
MASK_colorseg = closing(openBW,se2);
simshow(MASK_colorseg)
Out[0]:
No description has been provided for this image

Наложим инвертированную маску на исходное изображение знакомым способом:

In [ ]: 
sv_colorseg = StackedView(CV[1,:,:] + (.!MASK_colorseg./3), CV[2,:,:] + 
                            (.!MASK_colorseg./3), CV[3,:,:]);
view_colorseg = colorview(RGB, sv_colorseg)
Out[0]:
No description has been provided for this image

Фильтр Собеля

Фильтр Собеля — это оператор, используемый в обработке изображений для обнаружения границ (краёв объектов) на основе вычисления градиента яркости в горизонтальном и вертикальном направлениях. Он применяет дискретный аналог производной, подчёркивая области резкого изменения интенсивности пикселей.

Временно забываем про цвет и будем работать с картой яркости, то есть с градациями серого:

In [ ]: 
imgray = Gray.(I)
Out[0]:
No description has been provided for this image

Оператор Собеля использует два ядра (небольшие матрицы).

  • По оси X (вертикальные границы):

    [ -1  0  1 ]
[ -2  0  2 ]
[ -1  0  1 ]

  • По оси Y (горизонтальные границы):

    [ -1 -2 -1 ]
[  0  0  0 ]
[  1  2  1 ]

    Принцип работы:

    • Каждое ядро свертывается (convolution) с изображением, выделяя перепады яркости в своём направлении.

    • Результаты применяются к исходному изображению для получения:

      • Gx (градиент по X),

      • Gy (градиент по Y).

    • Общий градиент определяется как корень из суммы квадратов Gx и Gy

In [ ]: 
# Ядра Собеля для осей X и Y
sobel_x = Kernel.sobel()[1]  # Горизонтальный градиент (вертикальные границы)
sobel_y = Kernel.sobel()[2]  # Вертикальный градиент (горизонтальные границы)

# Применение свертки
gradient_x = imfilter(imgray, sobel_x)
gradient_y = imfilter(imgray, sobel_y)

# Общий градиент (объединение X и Y)
gradient_magnitude = sqrt.(gradient_x.^2 + gradient_y.^2);
imsobel = gradient_magnitude ./ maximum(gradient_magnitude)
Out[0]:
No description has been provided for this image

Бинаризуем результат фильтрации методом Отсу без дополнительных аргументов:

In [ ]: 
BW = binarize(imgray, Otsu());
simshow(BW)
Out[0]:
No description has been provided for this image

Добавим немного морфологической магии для получения результирующей маски:

In [ ]: 
se = Kernel.gaussian(3) .> 0.0035;
closeBW = closing(BW,se);
noobj = area_opening(closeBW; min_area=1000) .> 0;
se2 = Kernel.gaussian(5) .> 0.0027;
smooth = closing(noobj,se2);
smooth_2 = opening(smooth,se2);
MASK_sobel = area_opening(.!smooth_2; min_area=500) .> 0;
simshow(MASK_sobel)
Out[0]:
No description has been provided for this image

Визуализируем результат:

In [ ]: 
sv_sobel = StackedView(CV[1,:,:] + (.!MASK_sobel./3), CV[2,:,:] + 
                            (.!MASK_sobel./3), CV[3,:,:]);
view_sobel = colorview(RGB, sv_sobel)
Out[0]:
No description has been provided for this image
In [ ]: 
sv_sobel = StackedView(CV[1,:,:] + (MASK_sobel./3), CV[2,:,:] + 
                            (MASK_sobel./3), CV[3,:,:]);
view_sobel = colorview(imsobel, sv_sobel)

Фильтр Собеля использовался нами как ещё один текстурный фильтр. Мы воспользовались различием в равномерности изменения яркости областей леса и города.

SRG-алгоритм

Seeded Region Growing (SRG) — это классический алгоритм сегментации изображений, основанный на идее объединения пикселей или регионов в группы по схожим характеристикам, начиная с заранее заданных "зародышевых точек" (seed points).

Основные принципы работы:

  1. Выбор начальных точек (seeds)

    • Пользователь или автоматический метод выбирает одну или несколько стартовых точек (пикселей), которые принадлежат интересующему объекту.

  2. Определение критерия схожести

    • Обычно используется разница интенсивностей, цвета или текстуры между соседними пикселями.

    • Например, пиксель добавляется в регион, если его интенсивность отличается от среднего значения региона не более чем на порог T.

  3. Итеративное расширение регионов

    • Алгоритм проверяет соседей уже включенных пикселей и добавляет их в регион, если они удовлетворяют критерию.

    • Процесс продолжается, пока нельзя добавить новые пиксели.

В качестве начальных выберем точки с координатами [200, 50] и [200,300]. Оценим визуально их цвета - от этого сильно зависит успешность работы алгоритма:

In [ ]: 
[ I[200,50] I[200,300] ]
Out[0]:
No description has been provided for this image

Устанавливаем координаты начальных точек, и получаем сегменты (результат сегментации) функцией seeded_region_growing, также передав ей исходное цветное изображение. Средние значения цвета внутри сегмента можно узнать функцией segments_mean:

In [ ]: 
seeds = [(CartesianIndex(200,50),1), (CartesianIndex(200,300),2)]
segments = seeded_region_growing(I, seeds);
sm = segment_mean(segments);
[ sm[1] sm[2] ]
Out[0]:
No description has been provided for this image
In [ ]: 
simshow(map(i->segment_mean(segments,i), labels_map(segments)))
Out[0]:
No description has been provided for this image

Бинарную маску мы получим из матрицы лейблов. Нас интересуют пиксели с лейблом 1:

In [ ]: 
lmap = labels_map(segments);
BW_smart = lmap .== 1;
simshow(BW_smart)
Out[0]:
No description has been provided for this image

Ну и немного простой морфологии:

In [ ]: 
se = Kernel.gaussian(2) .> 0.004;
closeBW = closing(BW_smart,se);
MASK_smart = area_opening(.!closeBW; min_area=500) .> 0;
simshow(MASK_smart)
Out[0]:
No description has been provided for this image

Наложим на исходное изображение:

In [ ]: 
sv_smart = StackedView(CV[1,:,:] + (.!MASK_smart./3), CV[2,:,:] + 
                            (.!MASK_smart./3), CV[3,:,:]);
view_smart = colorview(RGB, sv_smart)
Out[0]: