Satellite image segmentation, part 2¶
The second example from the series on digital image processing. In Part 1 we considered the problem of satellite image segmentation and basic approaches to its solution. In the second part we will consider additional methods available to the novice developer, in particular:
colour segmentation
Sobel filter
SRG algorithm
Loading required libraries, importing image {#Loading-required-libraries,-importing-data}¶
The next line of code should be uncommented if something is missing from the applied code:
# Pkg.add(["Images", "ImageShow", "ImageContrastAdjustment", "ImageBinarization", "ImageMorphology", "ImageFiltering", "ImageSegmentation"])
Connect the necessary libraries:
using Images, ImageShow, ImageContrastAdjustment, ImageBinarization, ImageMorphology, ImageFiltering, ImageSegmentation
Loading the original colour image (satellite image):
I = load("$(@__DIR__)/map_small.jpg")
Colour segmentation¶
In the first part we hardly addressed colour information, but it is obvious that the forest area is mostly green, while the city contains a lot of grey and white. Which means we can try to extract regions that are close in colour to the shade of green observed in the image. We will binarise individual channels of the image, and we will access them, as before, using the function 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,:,:]) ]
A typical "green pixel" is located, as we can assume, in the centre of the image. Let's select the "central" pixel, knowing the width and height of the image, and look at its intensity values in the colour channels:
midh = Int(round(h/2));
midw = Int(round(w/2));
midpixel = CV[:,midh, midw]
Now, using standard logical comparison operations, we binarise each of the channels. In the red channel, pixels with intensities from 0.15 to 0.25 will correspond to a logical one, in the green channel - from 0.2 to 0.3, and in the blue channel - from 0.05 to 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])
And now let's combine three binary masks into one by the operation of logical "AND" - now the white pixel will be only where in the original colour image all three checked conditions (ranges) are fulfilled:
BIN = BIN_RED .& BIN_GREEN .& BIN_BLUE;
simshow(BIN)
Once again morphology will help us to make a mask out of this motley set of pixels. This time we will take a small structural element (7x7 pixels) in the shape of a rhombus:
se = strel_diamond((7,7))
And evaluate the result of the morphological closure operation:
closeBW = closing(BIN,se);
simshow(closeBW)
We obtain the resulting mask by removing small "blobs" and performing the closure operation again, now to smooth the boundaries of the main "large" areas of the merged white pixels:
openBW = area_opening(closeBW; min_area=500) .> 0;
se2 = Kernel.gaussian(3) .> 0.0025;
MASK_colorseg = closing(openBW,se2);
simshow(MASK_colorseg)
Let's apply the inverted mask to the original image in the familiar way:
sv_colorseg = StackedView(CV[1,:,:] + (.!MASK_colorseg./3), CV[2,:,:] +
(.!MASK_colorseg./3), CV[3,:,:]);
view_colorseg = colorview(RGB, sv_colorseg)
Sobel filter¶
The Sobel filter is an operator used in image processing to detect boundaries (edges of objects) by calculating the luminance gradient in the horizontal and vertical directions. It applies a discrete analogue of the derivative, highlighting areas of sharp changes in pixel intensity.
We will temporarily forget about colour and work with the luminance map, i.e. grayscale:
imgray = Gray.(I)
The Sobel operator uses two kernels (small matrices).
- The X axis (vertical bounds):
[ -1 0 1 ]
[ -2 0 2 ]
[ -1 0 1 ]
- Y axis (horizontal boundaries):
[ -1 -2 -1 ]
[ 0 0 0 ]
[ 1 2 1 ]
Principle of Operation:
Each kernel convolves (convolution) with the image, highlighting brightness variations in its direction.
The results are applied to the original image to obtain:
Gx (X gradient),
Gy (Y gradient).
The total gradient is defined as the root of the sum of the squares of Gx and 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)
Binarise the filtering result by Otsu method without additional arguments:
BW = binarize(imgray, Otsu());
simshow(BW)
Let's add some morphological magic to get the resulting mask:
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)
Visualise the result:
sv_sobel = StackedView(CV[1,:,:] + (.!MASK_sobel./3), CV[2,:,:] +
(.!MASK_sobel./3), CV[3,:,:]);
view_sobel = colorview(RGB, sv_sobel)
sv_sobel = StackedView(CV[1,:,:] + (MASK_sobel./3), CV[2,:,:] +
(MASK_sobel./3), CV[3,:,:]);
view_sobel = colorview(imsobel, sv_sobel)
We used the Sobel filter as another texture filter. We took advantage of the difference in brightness uniformity between the forest and city areas.
SRG algorithm¶
Seeded Region Growing (SRG) is a classical image segmentation algorithm based on the idea of grouping pixels or regions into groups based on similar characteristics, starting from predefined "seed points".
Basic principles of operation:
**Select seed points* 1.
- The user or an automatic method selects one or more starting points (pixels) that belong to the object of interest.
Similarity criterion definition
Typically used is the difference in intensity, colour or texture between neighbouring pixels.
For example, a pixel is added to a region if its intensity differs from the region average by no more than a threshold
T.
Iterative region expansion
The algorithm checks the neighbours of already included pixels and adds them to the region if they satisfy the criterion.
The process continues until no new pixels can be added.
Let's choose points with coordinates [200, 50] and [200,300] as the initial ones. Let's evaluate their colours visually - the success of the algorithm depends on it:
[ I[200,50] I[200,300] ]
Set the coordinates of the initial points, and get segments (segmentation result) by the function seeded_region_growing, also passing the initial colour image to it. The average colour values within a segment can be obtained by the function segments_mean:
seeds = [(CartesianIndex(200,50),1), (CartesianIndex(200,300),2)]
segments = seeded_region_growing(I, seeds);
sm = segment_mean(segments);
[ sm[1] sm[2] ]
simshow(map(i->segment_mean(segments,i), labels_map(segments)))
We get the binary mask from the label matrix. We are interested in pixels with the label 1:
lmap = labels_map(segments);
BW_smart = lmap .== 1;
simshow(BW_smart)
And some simple morphology:
se = Kernel.gaussian(2) .> 0.004;
closeBW = closing(BW_smart,se);
MASK_smart = area_opening(.!closeBW; min_area=500) .> 0;
simshow(MASK_smart)
Let's overlay the original image:
sv_smart = StackedView(CV[1,:,:] + (.!MASK_smart./3), CV[2,:,:] +
(.!MASK_smart./3), CV[3,:,:]);
view_smart = colorview(RGB, sv_smart)
SRG algorithm, although it is "smart", does not belong to machine learning techniques. It is a classical algorithm that has no learning steps and no automatic feature extraction directly from the data. However, it can be combined with machine learning techniques, for example, to automatically select starting points.
Conclusion¶
In the second part, we reviewed additional image processing techniques for the image segmentation task, such as thresholded colour segmentation, Sobel filter for texture analysis, and SRG algorithm.
Next, we are looking forward to machine learning for digital image processing tasks.