Satellite image segmentation, part 1¶
This is the first example in a series of scripts describing typical digital image processing tasks images. In the first part we will get acquainted with basic approaches and image processing algorithms for the task of segmenting a satellite image. Segmentation is the mapping of pixels of an image into different classes. In the simplest case it is possible to separate objects in the conventional "foreground" (foreground) from objects in the conventional "background". In this case, the result of segmentation can be a binary mask - a matrix of the same size as the original image. dimensions as the original image, where elements 1 correspond to the objects of interest objects of interest, and elements 0 correspond to the background. Topics covered:
- import and visualisation of digital images
- linear and non-linear filtering
- morphological operations
- development of user-defined functions
Loading necessary libraries, importing data¶
The next line of code should be uncommented, if something is missing:
# Pkg.add(["Images", "ImageShow", "ImageContrastAdjustment", "ImageBinarization", "ImageMorphology"])
Connecting loaded libraries:
using Images, ImageShow, ImageContrastAdjustment, ImageBinarization, ImageMorphology
Loading and visualising the source
colour image with the function load.
The image to be processed
is a satellite image
of a small section of the city map,
on which there are both
building areas and a forest belt.
Our goal is to try to separate
the forest from the city.
I = load("$(@__DIR__)/map_small.jpg") # загрузка исходного цветного изображения
Pre-processing¶
Typically, preprocessing tasks include analysing the image to identify Characteristic features of objects (distinguishing them from the background), changing contrast and/or brightness, filtering to remove noise or smooth textures, etc.
Let's see which of the preprocessing algorithms
will work in our case. Given that the image
is coloured, we can try to consider separate
channels (red, green and blue). For this purpose
let's use the function 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,:,:]) ]
In order to get the binary mask,
we need to work with the image in
greyscale image (called intensity matrix).
We can convert a colour image into
intensity matrix with the function Gray, or we can
consider one of the channels. As we can see,
in the blue channel, the forest area in the image
is more intense than the area
of the city. So the segmentation task is all about separating out
the brighter areas of pixels from the darker areas.
[RGB.(CV[3,:,:]) simshow(ones(h,20)) Gray.(I)]
Now we'll use a linear smoothing filter,
to equalise the intensity of the closely spaced
pixels in the blue channel. Let's compare two linear
filters, a Gaussian filter with
sliding window of 9 by 9 pixels, and a
a simple averaging filter with a 7 by 7 window.
Let's apply them to the matrix BLUE using the function
imfilterto which we pass different kernels.
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)]
With the chosen parameters of the used kernels,
the result is not fundamentally different and we are quite
satisfactory. As the main image for
for further operations we will take the result of a simple
averaging filter - the matrix avgIMG.
Binarisation¶
Binarisation is the process of obtaining a matrix, containing binary values, i.e. 1 or 0. In an image, this corresponds to white and black pixels. In segmentation, the white pixels will correspond to to one of the classes and black pixels to the other.
Let's compare two binarisation methods - Otsu's method,
available through the function binarize, and a simple
intensity matching to a threshold value.
binOTSU = binarize(avgIMG, Otsu()); # бинаризация методом Отсу
binTHRESH = avgIMG .> 0.25; # бинаризация по порогу 0.25
BW_AVG = binTHRESH;
[ simshow(binOTSU) simshow(zeros(h,20)) simshow(binTHRESH) ]
The results are also almost identical, and we'll continue with the binary image after a a simple comparison with the threshold. Even though we've got a binary matrix, a mask, there's a lot of noise on it: the areas have discontinuities, "gaps", there are many small objects. To make the mask more "smooth and free of noise, we're going to use binary morphology algorithms.
Morphological operations¶
Morphological operations can be performed both on binary images and on intensity matrices. But to understand the basic algorithms, we'll limit ourselves to binary morphology.
Basic algorithms such as dilatation, opening, closing, or erosion, they use special kernels - structural elements.
Parameterised functions are available to create typical rhombus or square structural elements, but we're interested in creating a circular structural element. element. So we binarise the two-dimensional Gaussian distribution at a threshold of 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) ]
Now let's perform a morphological closure operation
function closing
with a custom structural element. This will allow us to
to "close" small holes and gaps between white
pixels, as well as smooth the shape of objects (blobs).
We consider objects to be "joined" areas of white
pixels surrounded by black pixels, or the boundary of the
of the image.
After the close operation, we'll use
the area_open function and remove objects
that are less than 500 pixels.
closeBW = closing(BW_AVG,se); # морфологическое закрытие
openBW = area_opening(closeBW; min_area=500) .> 0; # удаление "малых" объектов
[ simshow(closeBW) simshow(zeros(h,20)) simshow(openBW) ]
Now invert our binary mask, and remove objects smaller than 5000 pixels. The inverted mask will also be smoothed with the operation closure:
fillBW = imfill(.!openBW, (0,5000)); # удаление "крупных" объектов
MASK_AVG = closing(fillBW,se); # морфологическое закрытие (сглаживание)
[ simshow(fillBW) simshow(ones(h,20)) simshow(MASK_AVG) ]
Segmentation result¶
We want to make sure that the resulting binary mask segments the original image. To do this, the images need to be combined.
For our chosen visualisation, we will add 30% each intensity of the binary mask to the red and green channels of the original colour image. Thus, we will get a semi-transparent mask of yellow colour:
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 ]
Visually we can evaluate the result as not bad. We observe successful selection of the forest area. The algorithm is not a universal algorithm, because many parameters algorithms were adjusted to similar test images. However, the purpose of the example is to to demonstrate the process of developing prototypes of DSC algorithms.
Nonlinear filtering¶
Let's consider another approach to preprocessing the of the image. In this case, we'll use image histogram stretching in grayscale, as well as a non-linear texture filter.
To begin with, let's apply to a monochrome image
function adjust_histogram with the method
LinearStretching. We got a very
contrasty image with a large number
of bright and dark pixels.
IMG = adjust_histogram(Gray.(I), LinearStretching(dst_minval = -1, dst_maxval = 1.8));
[ Gray.(I) simshow(ones(h,20)) IMG ]
Now let's use the so-called range filter. It's a non-linear filter that within a sliding kernel determines the difference between the maximum and minimum intensity value of the array.
For this purpose, let's write a simple user-defined function
rangewhich is based on the built-in function
extrema, which returns two values.
function range(array) # функция нелинейной фильтрации
(minval, maxval) = extrema(array); # мин. и макс. значения массива
return maxval - minval; # диапазон (макс. - мин.)
end
Let's apply the custom function to the sliding window
window using the function mapwindow.
The result of nonlinear filtering is binarised by Otsu method.
rangeIMG = mapwindow(range, IMG, (7,7)); # функция `range` в окне 7х7
BW_RNG = binarize(rangeIMG, Otsu()); # бинаризация методом Отсу
[ rangeIMG simshow(zeros(h,20)) BW_RNG ]
In the urban area. pixel intensity changes more strongly and spatially more frequently than in the forest belt area. Therefore, texture filters can be used to separating regions in images whose average intensity of which there is little difference after standard linear filters.
Custom morphology function¶
Let's collect the operations following binarisation
into one function myMorphology and apply it
to the binary image:
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
Visualise the resulting mask:
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 ]
Try to compare the results of processing with linear and non-linear filtering by software.
Custom segmentation function¶
Finally, let's gather all the sequential
preprocessing, nonlinear filtering.
and morphology into one more function mySegRange,
which takes a colour image as input:
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
And apply it to a larger image of the terrain:
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)
The practical meaning of such operations is the ability to perform analyses and calculations on selected objects in automatic mode. Knowing the linear dimensions of a pixel on a map in metres, for example, determine the area of objects on the map. In our case, we can calculate the area (in km^2), of the forest belt, the urban area, as well as the percentage of the forest to the total area of the area under consideration:
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) # процент леса
Conclusion¶
On the example of satellite image processing we have learnt different methods of image preprocessing, binarisation and morphological operations for the task of segmentation task, and learnt how to quickly develop custom digital image processing functions.