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

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

цветовая сегментация
фильтр Собеля
SRG-алгоритм

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

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

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

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

using Images, ImageShow, ImageContrastAdjustment, ImageBinarization, ImageMorphology, ImageFiltering, ImageSegmentation

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

I = load("$(@__DIR__)/map_small.jpg")

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

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

(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,:,:]) ]

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

midh = Int(round(h/2));
midw = Int(round(w/2));
midpixel = CV[:,midh, midw]

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.

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])

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

BIN = BIN_RED .& BIN_GREEN .& BIN_BLUE;
simshow(BIN)

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

se = strel_diamond((7,7))

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

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

closeBW = closing(BIN,se);
simshow(closeBW)

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

openBW = area_opening(closeBW; min_area=500) .> 0;
se2 = Kernel.gaussian(3) .> 0.0025;
MASK_colorseg = closing(openBW,se2);
simshow(MASK_colorseg)

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

sv_colorseg = StackedView(CV[1,:,:] + (.!MASK_colorseg./3), CV[2,:,:] + 
                            (.!MASK_colorseg./3), CV[3,:,:]);
view_colorseg = colorview(RGB, sv_colorseg)

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

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

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

imgray = Gray.(I)

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

По оси 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

# Ядра Собеля для осей 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)

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

BW = binarize(imgray, Otsu());
simshow(BW)

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

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)

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

sv_sobel = StackedView(CV[1,:,:] + (.!MASK_sobel./3), CV[2,:,:] + 
                            (.!MASK_sobel./3), CV[3,:,:]);
view_sobel = colorview(RGB, sv_sobel)