Документация 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]:
No description has been provided for this image

SRG-алгоритм, хоть и является "умным", не относится к техниками машинного обучения. Это классический алгоритм, в котором нет этапов обучения и автоматического выделения признаков непосредственно из данных. Впрочем, его можно объединять с техниками машинного обучения, например, для автоматического выбора начальных точек.

Заключение

Во второй части мы рассмотрели дополнительные методы обработки изображений для задачи сегментации изображения, такие как пороговая цветовая сегментация, фильтр Собеля для текстурного анализа, а также SRG-алгоритм.

Далее нас ждёт машинное обучение для задач цифровой обработки изображений.