Engee 文档
Notebook

卫星图像分割，第 1 部分

这是一系列脚本中的第一个示例 描述典型数字图像处理任务 图像。在第一部分中，我们将了解 基本方法和图像处理算法 分割卫星图像的基本方法和图像处理算法。 分割是将图像中的像素 映射成不同的类别。 在最简单的情况下，可以将物体从传统的 "前景 "中分离出来。 传统 "前景 "中的物体 (前景）中的物体与传统 "背景 "中的物体分开。 在这种情况下，分割的结果 可以是一个二元掩码--一个与原始图像大小相同的矩阵。 矩阵的尺寸与原始图像相同，其中 元素 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 像素的高斯滤波器和 简单的平均滤波器。 让我们使用函数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 进行形态闭合操作。这样我们就可以 封闭 "白色像素之间的小孔和间隙，以及平滑物体（blob）的形状。 像素之间的小洞和间隙，以及平滑物体（blobs）的形状。 我们认为物体是由白色像素 "连接 "起来的区域。 被黑色像素包围的白色像素区域，或图像的边界。 的边界。 关闭操作后，我们将使用 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

分割结果

我们要确保生成的 二进制掩码能分割原始图像。 为此，我们需要合并图像。

对于我们选择的可视化效果，我们将把二进制掩膜的每个 的强度。 绿色通道。 这样，我们将得到一个半透明的 黄色遮罩：

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

从视觉效果上看，效果还不错。 我们可以看到成功地选择了森林区域。 该算法并非通用算法，因为 许多参数 算法是根据类似的 测试图像。不过，这个例子的目的是 演示开发 DSC 算法原型的过程。 DSC 算法原型的开发过程。

非线性过滤

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

首先，让我们对单色图像应用 函数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

这些操作的实际意义在于 在自动模式下对选定对象进行分析和计算。 在自动模式下对选定对象进行分析和计算。 了解地图上像素的线性尺寸（以米为单位）、 例如 确定地图上物体的面积。 在我们的例子中，我们可以计算出森林带、城市区域和其他区域的面积（单位：平方公里）、 的面积（以平方公里为单位），以及 森林占总面积的百分比 的百分比：

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

结论

以卫星图像处理为例，我们学习了 不同的图像预处理方法、 二值化和形态学操作的不同方法。 分割任务的不同方法，并学会了如何快速开发 定制数字图像处理功能。