卫星图像的分割，第1部分

这是一系列脚本中的第一个例子
描述典型的数字
图像处理任务。在第一部分中，我们将熟悉
图像处理的基本方法和算法
于卫星图像的分割任务。
分割是映射
图像像素到不同的类。
在最简单的情况下，您可以分离
条件"前景"中的对象
从条件"背景"中的对象。
在这种情况下，分割的结果
可以是相同的二进制掩码矩阵
大小作为原始图像，其中
元素1对应
到感兴趣的对象，并且元素0对应于背景。
涵盖的主题:
*数字图像的导入和可视化
*线性和非线性滤波
*形态学操作
*开发自定义功能

下载必要的库，导入数据

下一行代码应该
如果应用程序中缺少某些内容，请注释掉。:

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

连接下载的库:

using Images, ImageShow, ImageContrastAdjustment, ImageBinarization, ImageMorphology

加载和可视化原件
具有功能的彩色图像 load.
处理后的图像
是卫星图像
城市地图的一小块区域,
其中包含两者
发展地区和森林带。
我们的任务是尝试分离
来自城市的森林。

I = load("$(@__DIR__)/map_small.jpg")       # загрузка исходного цветного изображения

预处理

作为一项规则，预处理任务
包括图像分析以突出显示
对象的特征（它们与背景的差异）,
对比度和/或亮度变化,
过滤以去除噪声或平滑纹理,
等等。

让我们考虑哪些预处理算法
在我们的情况下是合适的。鉴于图像
是在颜色，你可以尝试考虑分开
通道（红色，绿色和蓝色）。要做到这一点
，使用功能 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,:,:]) ]

为了得到二进制掩码,
我们需要与一个
灰度图像（所谓的强度矩阵）。
我们可以将彩色图像转换为
由函数的强度矩阵 Gray 或
考虑其中一个通道。正如我们所看到的,
在蓝色通道中，图像中的森林区域
强度不同于
的城区。分割的任务是分离
像素的较亮区域来自较暗区域。

[RGB.(CV[3,:,:]) simshow(ones(h,20)) Gray.(I)]

现在让我们使用线性平滑滤波器
使紧密间隔的强度均匀
蓝色通道中的像素。让我们比较两个线性
滤波器-高斯滤波器
一个9乘9像素的滑动窗口，以及
作为一个简单的平均滤波器,窗口测量7乘7.
让我们将它们应用于矩阵 BLUE 使用函数
imfilter，我们将转移不同的内核。

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

使用所用内核的选定参数,
结果并没有根本不同，我们相当
满意。作为主要图像
进一步的操作，我们采取一个简单的结果
平均滤波器，矩阵 avgIMG.

二值化

二值化是获得矩阵的过程
包含二进制值，即1或0。
在图像中，这对应于白色和
黑色像素。在分割中，白色像素
会对应
到其中一个类，而黑色像素到另一个。

让我们比较两种二值化方法--Otsu方法,
可通过功能 binarize，和一个简单的
强度与阈值的比较。

binOTSU = binarize(avgIMG, Otsu());         # бинаризация методом Отсу
binTHRESH = avgIMG .> 0.25;                 # бинаризация по порогу 0.25
BW_AVG = binTHRESH;
[ simshow(binOTSU) simshow(zeros(h,20)) simshow(binTHRESH) ]

结果也几乎相同,
之后我们将继续使用二进制图像
与阈值的简单比较。
尽管我们得到了一个二进制文件
矩阵，一个面具，上面有很多噪音：区域
有缺口，"狭缝"，有很多
小物体。制作面具
更平滑，无噪音，我们将使用
二元形态学算法。

形态学操作

可以进行形态学操作
在二进制图像和
强度矩阵。但要了解基本的
算法，我们将自己限制为二进制
形态。

基本算法，如扩张,
打开，关闭，或侵蚀，使用
特殊内核-结构元素。

我们可以创建参数化函数
典型的"菱形"或"方形"型结构元件,
但是，我们有兴趣创建一个圆形结构
元素。因此，我们二值化二维
阈值为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) ]

现在让我们执行形态学闭合操作
与功能 closing
具有自定义结构元素。这将允许
你要"关闭"白色之间的小孔和缝隙
像素，以及平滑物体（blob）的形状。
我们认为物体是白色的"组合"区域
由黑色像素包围的像素，或边框
的图像。
关闭操作后，我们将使用
的功能 area_open 我们将删除大小的对象
不超过500像素。

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

现在我们反转我们的二进制掩码
并移除尺寸小于5000像素的对象。
我们也将平滑倒面具
随着关闭操作。:

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

分割结果

我们要确保由此产生的结果
二值掩码对原始图像进行分段。
为此，需要对图像进行组合。

对于我们选择的可视化，我们将添加30%
二进制掩模的强度到红色和
原始彩色图像的绿色通道。
因此，我们得到一个半透明
黄色面具。:

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 ]

在视觉上，我们可以评估结果为好。
我们正在目睹森林面积的成功分配。
该算法不是通用的，因为
许多人
算法参数匹配相似
测试图像。然而，示例的目的是
演示开发原型的过程
中央研究所的算法。

非线性滤波

让我们考虑另一种方法
图像预处理。在这种情况下，我们将使用
图像的直方图的拉伸
在灰度，以及一个非线性
纹理过滤器。

首先，让我们将该函数应用于单色图像
adjust_histogram 用的方法
LinearStretching. 我们得到了一个非常
对比图像与大量的
亮和暗像素。

IMG = adjust_histogram(Gray.(I), LinearStretching(dst_minval = -1, dst_maxval = 1.8));
[ Gray.(I) simshow(ones(h,20)) IMG ]

现在让我们使用所谓的范围过滤器。
这是一个非线性滤波器，
确定滑动核心内部阵列的最大强度和最小强度之间的差值
.

为此，让我们编写一个简单的自定义函数
range，这是基于内置功能
extrema，它返回两个值。

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

让我们将一个自定义函数应用于窗口滑动
使用函数跨越图像 mapwindow.
通过Otsu方法对非线性滤波的结果进行二值化。

rangeIMG = mapwindow(range, IMG, (7,7));        # функция `range` в окне 7х7
BW_RNG = binarize(rangeIMG, Otsu());            # бинаризация методом Отсу
[ rangeIMG simshow(zeros(h,20)) BW_RNG ]

在城市发展领域
，像素强度变化更强烈且
空间上比森林带区域更频繁。
因此，纹理滤镜可用于
图像中平均强度的单独区域
标准后略有变化
线性滤波器。

自定义形态函数

让我们结合二值化后的操作
成一个功能。 myMorphology 并应用它
到二值图像。:

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

myMorphology (generic function with 1 method)

可视化生成的蒙版:

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 ]

尝试比较linear的结果
和非线性滤波编程。

自定义分割功能

最后，我们将组装所有顺序
预处理，非线性滤波
和形态学算法转换成另一个函数。 mySegRange,
它接受彩色图像作为输入:

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

mySegRange (generic function with 1 method)

并将其应用于该区域的更大快照。:

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)

这种操作的实际意义是能力
对选定的数据进行分析和计算
自动模式下的对象。
了解地图上像素的线性尺寸（以米为单位），您可以,
例如,
确定对象在地图上的区域。
在我们的例子中，我们可以计算面积（以km^2为单位）
被森林带占据，城市发展，以及
作为森林占总面积的百分比
正在考虑的地点。:

pixsquare = (800/74)^2;                   # площадь одного пикселя растра карты (км^2)
forest_area_km2 = round((sum(myMASK) * pixsquare) / 1e6, sigdigits = 4) # площадь леса

17.73

city_area_km2 =  round((sum(.!myMASK) * pixsquare) / 1e6, sigdigits = 4) # площадь города

41.08

forest_percentage = round(100*sum(myMASK)/length(I2), sigdigits = 4)    # процент леса

30.15

结论

利用卫星图像处理的例子，我们认识了
用各种方法进行图像预处理,
二值化和形态学操作
细分任务，并且还学会了如何快速开发
自定义数字图像处理功能。