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

Автор
avatar-ussmoussmo
Notebook

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

Это первый пример из цикла скриптов, описывающих типичные задачи цифровой обработки изображений. В первой части мы познакомимся с базовыми подходами и алгоритмами обработки изображений для задачи сегментации спутникового снимка. Сегментация - это сопоставление пикселей изображения различным классам. В самом простом случае можно отделять объекты на условном "переднем плане" (foreground) от объектов условного "фона" (background). В этом случае результатом сегментации может служить бинарная маска - матрица таких же размеров, что и исходное изображение, где элементы 1 соответствуют интересующим нас объектам, а элементы 0 - фону. Рассматриваемые темы:

  • импорт и визуализация цифровых изображений
  • линейная и нелинейная фильтрация
  • морфологические операции
  • разработка пользовательских функций

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

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

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

Подключение загруженных библиотек:

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

Загрузка и визуализация исходного цветного изображения функцией load. В качестве обрабатываемого изображения выступает спутниковый снимок небольшого участка карты города, на котором присутствуют как области застройки, так и лесополоса. Наша задача - попробовать отделить лес от города.

In [ ]:
I = load("$(@__DIR__)/map_small.jpg")       # загрузка исходного цветного изображения
Out[0]:
No description has been provided for this image

Предобработка

Как правило, в задачи предобработки входят анализ изображения для выделения характерных признаков объектов (отличия их от фона), изменения контраста и/или яркости, фильтрация для удаления шумов или сглаживания текстуры, и проч.

Рассмотрим, какие из алгоритмов предобработки подойдут в нашем случае. Учитывая, что изображение цветное, можно попробовать рассмотреть отдельные каналы (красный, зелёный и синий). Для этого воспользуемся функцией channelview:

In [ ]:
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,:,:]) ]
Out[0]:
No description has been provided for this image

Для того, чтобы получить бинарную маску, нам нужно работать с изображением в градациях серого (так называемой матрицей интенсивности). Мы можем перевести цветное изображение в матрицу интенсивности функцией Gray, или же рассмотреть один из каналов. Как мы видим, в синем канале область леса на изображении сильнее отличается по интенсивности от области города. И задача сегментации сводится к отделению более ярких областей пикселей от более тёмных.

In [ ]:
[RGB.(CV[3,:,:]) simshow(ones(h,20)) Gray.(I)]
Out[0]:
No description has been provided for this image

Теперь воспользуемся линейным сглаживающим фильтром, чтобы выровнять интенсивность близко расположенных пикселей в синем канале. Сравним два линейных фильтра - фильтр Гаусса со скользящим окном размером 9 на 9 пикселей, а также простой осредняющий фильтр с окном размером 7 на 7. Применим их к матрице BLUE при помощи функции imfilter, которой мы передадим разные кернелы.

In [ ]:
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)]
Out[0]:
No description has been provided for this image

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

Бинаризация

Бинаризация - это процесс получения матрицы, содержащий бинарные значения, то есть 1 или 0. На изображении это соответствует белым и чёрным пикселям. В сегментации белые пиксели будут соответствовать одному из классов, а черные - другому.

Сравним два метода бинаризации - метод Отсу, доступный посредством функции binarize, и простое сопоставление интенсивности пороговому значению.

In [ ]:
binOTSU = binarize(avgIMG, Otsu());         # бинаризация методом Отсу
binTHRESH = avgIMG .> 0.25;                 # бинаризация по порогу 0.25
BW_AVG = binTHRESH;
[ simshow(binOTSU) simshow(zeros(h,20)) simshow(binTHRESH) ]
Out[0]:
No description has been provided for this image

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

Морфологические операции

Морфологические операции можно осуществлять как над бинарными изображениями, так и над матрицами интенсивности. Но для понимания базовых алгоритмов будем ограничиваться бинарной морфологией.

В основных алгоритмах, таких как дилатация, раскрытие, закрытие или эрозия, используются специальные кернелы - структурные элементы.

Нам доступны параметризованные функции для создания типичных структурных элементов типа "ромб" или "квадрат", впрочем, нам интересно создать круглый структурный элемент. Поэтому мы бинаризуем двумерное распределение Гаусса по порогу в 0.0025.

In [ ]:
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) ]
Out[0]:
No description has been provided for this image

Теперь осуществим операцию морфологического закрытия функцией closing с пользовательским структурным элементом. Это позволит "закрыть" небольшие отверстия и щели между белыми пикселями, а также сгладить форму объектов (блобов). Объектами мы считаем "объединённые" области белых пикселей, окружённых черными пикселями, или границей изображения. После операции закрытия мы воспользуемся функцией area_open и удалим объекты, размер которых не превышает 500 пикселей.

In [ ]:
closeBW = closing(BW_AVG,se);                           # морфологическое закрытие
openBW = area_opening(closeBW; min_area=500) .> 0;      # удаление "малых" объектов
[ simshow(closeBW) simshow(zeros(h,20)) simshow(openBW) ]
Out[0]:
No description has been provided for this image

Теперь инвертируем нашу бинарную маску, и удалим объекты размером менее 5000 пикселей. Инвертированную маску также сгладим операцией закрытия:

In [ ]:
fillBW = imfill(.!openBW, (0,5000));            # удаление "крупных" объектов
MASK_AVG = closing(fillBW,se);                  # морфологическое закрытие (сглаживание)
[ simshow(fillBW) simshow(ones(h,20)) simshow(MASK_AVG) ]
Out[0]:
No description has been provided for this image

Результат сегментации

Мы хотим убедиться, что полученная бинарная маска сегментирует исходное изображение. Для этого изображения нужно совместить.

Для выбранной нами визуализации добавим по 30% интенсивности бинарной маски в красный и зелёный каналы исходного цветного изображения. Таким образом, получим полупрозрачную маску жёлтого цвета:

In [ ]:
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 ]
Out[0]:
No description has been provided for this image

Визуально можем оценить результат, как неплохой. Мы наблюдаем успешное выделение лесного массива. Алгоритм не является универсальным, так как многие параметры алгоритмов подбирались к похожим тестовым изображениям. Впрочем, целью примера является демонстрация процесса разработки прототипов алгоритмов ЦОИ.

Нелинейная фильтрация

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

Для начала применим к монохромному изображению функцию adjust_histogram с методом LinearStretching. Мы получили очень контрастное изображение с большим числом ярких и тёмных пикселей.

In [ ]:
IMG = adjust_histogram(Gray.(I), LinearStretching(dst_minval = -1, dst_maxval = 1.8));
[ Gray.(I) simshow(ones(h,20)) IMG ]
Out[0]:
No description has been provided for this image

Теперь воспользуемся так называемым фильтром диапазона. Это нелинейный фильтр, который внутри скользящего кернела определяет разницу между максимальным и минимальным значением интенсивности массива.

Для этого напишем простую пользовательскую функцию range, в основе которой лежит встроенная функция extrema, возвращающая два значения.

In [ ]:
function range(array)                           # функция нелинейной фильтрации
    (minval, maxval) = extrema(array);          # мин. и макс. значения массива
    return maxval - minval;                     # диапазон (макс. - мин.)
end

Применим пользовательскую функцию к скользящему по изображению окну при помощи функции mapwindow. Результат нелинейной фильтрации бинаризуем методом Отсу.

In [ ]:
rangeIMG = mapwindow(range, IMG, (7,7));        # функция `range` в окне 7х7
BW_RNG = binarize(rangeIMG, Otsu());            # бинаризация методом Отсу
[ rangeIMG simshow(zeros(h,20)) BW_RNG ]
Out[0]:
No description has been provided for this image

В области городской застройки интенсивность пикселей меняется сильнее и пространственно чаще, чем в области лесополосы. Поэтому текстурные фильтры могут применяться для разделения областей на изображениях, средняя интенсивность которых слабо различается после стандартных линейных фильтров.

Пользовательская функция морфологии

Соберём последующие за бинаризацией операции в одну функцию myMorphology и применим её к бинарному изображению:

In [ ]:
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
Out[0]:
myMorphology (generic function with 1 method)

Визуализируем результирующую маску:

In [ ]:
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 ]
Out[0]:
No description has been provided for this image

Попробуйте сравнить результаты обработки при линейной и нелинейной фильтрации программным образом.

Пользовательская функция сегментации

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

In [ ]:
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
Out[0]:
mySegRange (generic function with 1 method)

И применим её к более крупному снимку местности:

In [ ]:
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)
Out[0]:
No description has been provided for this image

Практический смысл подобных операций - это возможность осуществлять анализ и вычисления над выделенными объектами в автоматическом режиме. Зная линейные размеры пикселя на карте в метрах, можно, к примеру, определить площадь объектов на карте. В нашем случае мы можем подсчитать площадь (в км^2), занимаемую лесополосой, городской застройкой, а также процентное отношение леса к общей площади рассматриваемого участка:

In [ ]:
pixsquare = (800/74)^2;                   # площадь одного пикселя растра карты (км^2)
forest_area_km2 = round((sum(myMASK) * pixsquare) / 1e6, sigdigits = 4) # площадь леса
Out[0]:
17.73
In [ ]:
city_area_km2 =  round((sum(.!myMASK) * pixsquare) / 1e6, sigdigits = 4) # площадь города
Out[0]:
41.08
In [ ]:
forest_percentage = round(100*sum(myMASK)/length(I2), sigdigits = 4)    # процент леса
Out[0]:
30.15

Заключение

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