前言

本系列文章为学校课程《数字图像处理》布置的一些小project作业

除了给出每个问题的解答和代码，我也会附上相关知识点，以方便后续复习

前置概念

形态学处理

概念

形态学处理是指将数学形态学作为⼯具从图像中提取对于表达和描述区域形状有⽤的图像分量。

形态学处理表现为⼀种邻域运算形式；⼀种特殊定义的邻域称之为“结构单元”（Structure Element），在每个像素位置上它与⼆值图像对应的区域进⾏特定的逻辑运算，逻辑运算的结果为输出图像的相应像素。

形态学运算的效果取决于结构单元的⼤⼩、内容以及逻辑运算的性质

基本处理定义

二值形态学处理

二值图像中，黑色表示1，白色表示0

平移

此处的结构单元x在平移运算中通常为1个点

扩张

扩张使得图像扩大

腐蚀

腐蚀会使得图像缩小

如果有完美匹配则将结构中心处设为1

性质

扩张和腐蚀具有以下性质

方形结构单元越大，膨胀和腐蚀的效果越强

形态学变换

结构开(open)变换

结构开变换等同于先腐蚀然后再扩张

公式：$A \circ B = (A \ominus B) \oplus B$

目的：使轮廓平滑，抑制A物体边界的⼩离散点或尖峰，在研究物体的形态分布时常⽤。⽤来消除⼩物体、在纤细点处分离物体、平滑较⼤物体的边界的同时并不明显改变其⾯积。

下图能够非常直观的说明平滑作用

结构闭(close)变换

结构闭变换等同于先扩张然后再腐蚀

公式：$A \bullet B = (A \ominus B) \oplus B$

目的：也是⽤于图像光滑。但与开变换相反，闭变换⽤来填充物体内细⼩空洞、连接邻近物体、平滑其边界的同时并不明显改变其⾯积。

下图能够非常直观的说明填补空洞的作用

图像分割

概述

图像分割即将数字图像分割成多个部分或区域的过程。这些区域通常是根据图像中的像素特性（如颜色、亮度、纹理）来定义的，目的是使分割后的区域在某种意义上具有内部一致性，同时与其他区域有明显的区别。

图像分割算法可以看作是基于亮度值的不连续性和相似性

基于阈值的分割

基础

基于阈值的图像分割是一种将图像转换为二值图像的方法，即图像中的每个像素点要么是前景（通常表示为白色或1），要么是背景（通常表示为黑色或0）。这种方法特别适用于对象和背景在亮度或颜色上有显著差异的图像。

在基于阈值的分割中，关键是确定一个或多个阈值，以便将像素分为不同的类别。不同类型的阈值分割方法侧重于如何选择这些阈值。

基本全局阈值

当目标和背景像素的灰度分布非常不同时，可以对整个图像使用单个（全局）阈值

下面是一种全局阈值常用的迭代算法：

为全局阈值T选择一个初始估计值
使用T分割图像，这将产生两组像素：由灰度值大于T的所有像素组成的$G_1$，由所有小于等于T的像素组成的$G_2$
对$G_1$和$G_2$中的像素分别计算平均灰度值$m_1$和$m_2$
在$m_1$和$m_2$之间计算一个新的阈值：
$T= \frac{1}{2} (m_1+m_2)$
重复步骤2到步骤4，直到连续迭代中的两个T值间的差小于某个预定义的值$\Delta T$为止

当与目标和背景相关的直方图模式之间存在一个非常清晰的波谷时，上述算法很有效

基本自适应阈值（可变阈值）

可变阈值，也称为局部或自适应阈值，其基本思想是为图像的不同区域使用不同的阈值。这种方法认识到图像的各个部分可能由于光照、阴影、噪声或其他因素的影响，而具有不同的亮度和对比度特性。因此，相比于全局阈值，它可以更加精细地处理图像的局部特征。

自适应阈值相对于全局阈值的优势体现在：

光照不均匀：

自适应阈值能够更好地处理图像中光照变化的区域，因为它可以为每个局部区域计算一个专门的阈值。
局部对比度变化：

在一些区域可能存在的高对比度物体旁边的低对比度物体可以通过自适应阈值更容易被分割出来。
噪声和纹理：

全局阈值可能会将噪声错误地识别为前景或背景。而自适应阈值可以根据局部区域的纹理和噪声水平调整阈值，从而降低噪声对分割结果的影响。

一个较为简单的自适应阈值算法如下：

初始化：选择一个窗口大小windowSize，这将定义用于计算局部阈值的邻域大小。选择一个调整常数C，它将用于从局部平均值中减去，以计算每个像素的阈值。
遍历图像的每个像素：对于图像中的每个像素p(i, j)，执行以下步骤：

a. 确定局部邻域：计算出当前像素p(i, j)的局部邻域。这个邻域是以(i, j)为中心，大小为windowSize x windowSize的像素块。

b. 提取局部区域：从原始图像中提取出步骤2a定义的局部邻域。如果邻域超出图像边界，则只考虑落在图像内的部分。
计算局部统计量：对于步骤2b中提取的每个局部邻域，计算局部统计量，本算法中即计算其所有像素的平均灰度值m。
计算局部阈值：使用局部统计量和调整常数C来计算当前像素的局部阈值T_local。阈值可以用以下公式计算：T_local = m * (1 - C)。
应用阈值：对于每个像素p(i, j)，根据其局部阈值T_local来决定其在二值图像中的值。如果p(i, j)的灰度值大于T_local，则在二值图像中该像素设置为1（或255，取决于所需的二值图像格式），否则设置为0。
输出二值图像：重复步骤2至步骤5，直到图像中的每个像素都被处理。最终得到的二值图像是根据局部阈值分割后的结果。

HW1 二值形态学处理

作业要求：

一个二进制数组代表了一张黑白图像的一部分，如下图所示。假设该片段周围的所有像素都是黑色背景，请在这部分图像上执行下列操作。

(a) 使用如下结构元素进行扩张（Dilation）操作，原点用圆圈标出。

(b) 使用如下结构元素进行腐蚀（Erosion）操作，原点用圆圈标出。

(c) 使用如下结构元素进1行扩张（Dilation）操作，原点用圆圈标出。

(d) 使用如下结构元素进行腐蚀（Erosion）操作，原点用圆圈标出。

(e) 使用上述结构元素进行开操作（Opening）。

(f) 使用上述结构元素进行闭操作（Closing）。

这里我们将数组视为二值图像来展示

首先我们定义题目中待处理的的图像数组以及结构单元

binary_image = [
    0 0 0 0 0 0 0;
    0 0 1 1 0 0 0;
    0 0 0 1 0 1 0;
    0 0 0 1 1 0 0;
    0 0 1 1 1 1 0;
    0 1 0 1 0 1 0;
    0 0 0 0 0 0 0
];
se1 = [0 1 1 1 0];
se2 = [1 0; 0 1];

扩张操作使用matlab自带的函数imdilate()

dilation1 = imdilate(binary_image, se1);
imwrite(dilation1, 'dilation1.png');
dilation2 = imdilate(binary_image, se2);
imwrite(dilation2, 'dilation2.png');

腐蚀操作使用函数imerode()

erosion1 = imerode(binary_image, se1);
imwrite(erosion1, 'erosion1.png');
erosion2 = imerode(binary_image, se2);
imwrite(erosion2, 'erosion2.png');

开操作使用函数imopen()

opening1 = imopen(binary_image, se1);
imwrite(opening1, 'opening1.png');
opening2 = imopen(binary_image, se2);
imwrite(opening2, 'opening2.png');

闭操作使用函数imclose()

closing1 = imclose(binary_image, se1);
imwrite(closing1, 'closing1.png');
closing2 = imclose(binary_image, se2);
imwrite(closing2, 'closing2.png');

最后将形态学处理的结果与原图进行对比

figure;
subplot(3,3,1), imshow(binary_image), title('原图');
subplot(3,3,2), imshow(dilation1), title('使用SE1扩张');
subplot(3,3,3), imshow(erosion1), title('使用SE1腐蚀');
subplot(3,3,4), imshow(dilation2), title('使用SE2扩张');
subplot(3,3,5), imshow(erosion2), title('使用SE2腐蚀');
subplot(3,3,6), imshow(opening1), title('使用SE1进行开操作');
subplot(3,3,7), imshow(opening2), title('使用SE2进行开操作');
subplot(3,3,8), imshow(closing1), title('使用SE1进行闭操作');
subplot(3,3,9), imshow(closing2), title('使用SE2进行闭操作');

得到的结果如下图

我们可以看到，由于SE1结构单元的元素呈水平排列，所以经过扩张操作，原图的水平方向上白色区域变宽

同时SE1腐蚀导致原图水平方向上变窄且细节消失

SE2由于是呈对角线，所以扩张操作增加了对角方向的白色区域，而腐蚀减少了对角方向的白色像素

SE1的开操作将上半部分小的白色趋于滤去，而SE2的开操作相对而言保留了更多不同的区域

而SE1结构单元的闭操作填补了原图的细小空洞，而SE2闭操作后与原图一致，这是由于对角方向上基本上都是连贯的，不需要填补

HW2 图像分割

作业要求：

图1展示了blobz1.png和blobz2.png两幅图像。两者的区别在于blobz1的照明几乎是均匀的，而blobz2的照明非常不均匀。这个问题的目标是构建一个基于全局灰度阈值的算法，用于分割每个图像。

blobz1.png如下

blobz2.png如下

首先我们需要观察一下这两张图片的直方图，用于确认什么分割方法更适合

img1 = imread('blobz1.png');
img2 = imread('blobz2.png');

% 显示原始图像及其直方图
figure;
subplot(2,2,1); imshow(img1); title('blobz1.png');
subplot(2,2,2); imshow(img2); title('blobz2.png');
subplot(2,2,3); imhist(img1); title('blobz1 直方图');
subplot(2,2,4); imhist(img2); title('blobz2 直方图');
saveas(gcf, 'histograms.png');

得到下图结果

可以发现，由于图一光照较为均匀，直方图中呈现出了两个较为明显的波峰，适合直接使用全局阈值分割

而图二的波形由于光照不均匀是多峰的，所以需要使用自适应阈值分割来处理

首先是全局阈值算法，我们使用简单的迭代算法：

为全局阈值T选择一个初始估计值
使用T分割图像，这将产生两组像素：由灰度值大于T的所有像素组成的$G_1$，由所有小于等于T的像素组成的$G_2$
对$G_1$和$G_2$中的像素分别计算平均灰度值$m_1$和$m_2$
在$m_1$和$m_2$之间计算一个新的阈值：
$T= \frac{1}{2} (m_1+m_2)$
重复步骤2到步骤4，直到连续迭代中的两个T值间的差小于某个预定义的值$\Delta T$为止

首先设置初始参数

1
2
3

T1 = double(max(img1(:)) + min(img1(:))) / 2;
deltaT1 = 0.5; % 阈值变化的停止条件
prevT1 = -1; % 初始化为一个负数

这里T的初始值我们设置为使用最大值和最小值计算的平均值

停止条件设置为0.5

接下来是循环求取阈值

while abs(T1 - prevT1) >= deltaT1
    G1 = img1(img1 > T1);
    G2 = img1(img1 <= T1);
    m1 = mean(G1(:));
    m2 = mean(G2(:));
    prevT1 = T1;
    T1 = 0.5 * (m1 + m2);
end

求得阈值后，我们直接使用判断语句得到结果的图像矩阵

1 2	binaryImage1 = img1 > T1; imwrite(binaryImage1, 'result/global_img1.png');

我们可以得到如下和原图的对比

全局迭代阈值算法对于图一的分割效果非常好

但是到了图二，全局阈值的分割效果就很不好了，下图为图二的全局阈值处理结果

可以看到，由于原图的左下角光照较差，大部分像素点的像素值较低，处于直方图左侧，低于算出的阈值，被作为前景处理

接下来我们使用自适应算法：

初始化：选择一个窗口大小windowSize，这将定义用于计算局部阈值的邻域大小。选择一个调整常数C，它将用于从局部平均值中减去，以计算每个像素的阈值。
遍历图像的每个像素：对于图像中的每个像素p(i, j)，执行以下步骤：

a. 确定局部邻域：计算出当前像素p(i, j)的局部邻域。这个邻域是以(i, j)为中心，大小为windowSize x windowSize的像素块。

b. 提取局部区域：从原始图像中提取出步骤2a定义的局部邻域。如果邻域超出图像边界，则只考虑落在图像内的部分。
计算局部统计量：对于步骤2b中提取的每个局部邻域，计算局部统计量，本算法中即计算其所有像素的平均灰度值m。
计算局部阈值：使用局部统计量和调整常数C来计算当前像素的局部阈值T_local。阈值可以用以下公式计算：T_local = m * (1 - C)。
应用阈值：对于每个像素p(i, j)，根据其局部阈值T_local来决定其在二值图像中的值。如果p(i, j)的灰度值大于T_local，则在二值图像中该像素设置为1（或255，取决于所需的二值图像格式），否则设置为0。
输出二值图像：重复步骤2至步骤5，直到图像中的每个像素都被处理。最终得到的二值图像是根据局部阈值分割后的结果。

这里我们将自适应算法写为一个函数，窗口大小以及C作为参数传入

function binaryImage = adaptiveThresholding(img, windowSize, C)
    halfWindowSize = floor(windowSize / 2);
    binaryImage = zeros(size(img));
    for i = 1:size(img, 1)
        for j = 1:size(img, 2)
			...
        end
    end
    binaryImage = uint8(binaryImage * 255);
end

由于每个点都需要计算相应阈值，我们需要循环遍历

rowMin = max(1, i - halfWindowSize);
rowMax = min(size(img, 1), i + halfWindowSize);
colMin = max(1, j - halfWindowSize);
colMax = min(size(img, 2), j + halfWindowSize);
localRegion = img(rowMin:rowMax, colMin:colMax);
localMean = mean(localRegion(:));
binaryImage(i, j) = img(i, j) > localMean * (1 - C);