卫星图像分割,第 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 像素的高斯滤波器和
简单的平均滤波器。
让我们使用函数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。 在图像中,这对应于白色和 黑色像素。在分割过程中,白色像素 将对应于 其中一类,而黑色像素则对应另一类。
让我们比较两种二值化方法--大津方法、
大津的二值化方法(可通过函数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)的形状。
像素之间的小洞和间隙,以及平滑物体(blobs)的形状。
我们认为物体是由白色像素 "连接 "起来的区域。
被黑色像素包围的白色像素区域,或图像的边界。
的边界。
关闭操作后,我们将使用
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) ]
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 ]
从视觉效果上看,效果还不错。 我们可以看到成功地选择了森林区域。 该算法并非通用算法,因为 许多参数 算法是根据类似的 测试图像。不过,这个例子的目的是 演示开发 DSC 算法原型的过程。 DSC 算法原型的开发过程。
非线性过滤¶
让我们考虑另一种方法来预处理 的另一种方法。在这种情况下,我们将使用 图像直方图拉伸 以及非线性纹理过滤器。 纹理过滤器。
首先,让我们对单色图像应用
函数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
窗口。
非线性滤波的结果将通过大津方法进行二值化处理。
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
将生成的掩码可视化:
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 ]
尝试比较用线性 和非线性滤波处理的结果。
自定义分割功能¶
最后,让我们把所有的顺序
预处理、非线性过滤。
和形态学处理的所有功能,再整合成一个函数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
并将其应用到更大的地形图像中:
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)
这些操作的实际意义在于 在自动模式下对选定对象进行分析和计算。 在自动模式下对选定对象进行分析和计算。 了解地图上像素的线性尺寸(以米为单位)、 例如 确定地图上物体的面积。 在我们的例子中,我们可以计算出森林带、城市区域和其他区域的面积(单位:平方公里)、 的面积(以平方公里为单位),以及 森林占总面积的百分比 的百分比:
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) # процент леса
结论¶
以卫星图像处理为例,我们学习了 不同的图像预处理方法、 二值化和形态学操作的不同方法。 分割任务的不同方法,并学会了如何快速开发 定制数字图像处理功能。