【图像处理笔记】图像分割之基于灰度不连续性的分割方法
形态学处理相同，图像分割操作的输入是图像，输出是从图像中提取出来的属性。本章的大多数分割算法都基于图像灰度值的两个基本性质之一：不连续性和相似性。第一类方法根据灰度的突变（如边缘）将图像分割为多个区域；第二类方法根据一组预定义的准则把一幅图像分割为多个区域。本节的重点是以灰度局部剧烈变化检测为基础的分割方法。我们感兴趣的三类图像特征是孤立的点、线和边缘。边缘像素是图像中灰度突变的那些像素，而边缘（或边缘线段）是相连边缘像素的集合。边缘检测子是用来检测边缘像素的局部图像处理工具。线可视为一条细边缘线段，其两侧的背景灰度要么远大于线像素的灰度，要么远小于线像素的灰度。最后，孤立点可视为一个背景（前景）像素围绕的前景（背景）像素。
1 基础知识
1.1 分割
图像分割需满足以下条件：
a)分割必须是完全的，即图像为所有子区域的并集，每个像素都必须在一个区域内
b)所有子区域都是一个连通集，即一个区域中的点以某些预定义的方式连接（如这些点必须是8连通的）
c)各个区域必须是不相交的
d)分割后的区域中的像素必须满足某些性质，如某一区域中的所有像素都有相同的灰度
e)两个邻接区域像素满足的性质是不同的
分割中的基本问题就是把一幅图像分成满足前述条件的多个区域。通常单色图像的分割算法依据的是灰度值的两个性质之一：不连续性和相似性。
在基于灰度值不连续性的第一类算法中，我们假设这些区域的边界彼此完全不同，并且与背景不同，以便能够根据灰度的局部不连续性来检测边界，主要方法是基于边缘的分割。下图(a)显示了一个恒定灰度区域叠加到一个恒定灰度暗色背景上的一幅图像，这两个区域组成了整幅图像。图(b)是根据灰度不连续性计算内部区域的边界的结果。边界内侧和外侧都是0，因为在这些区域中灰度不存在不连续性。图(c)是为了分割图像，我们对边界上或边界内的像素分配一个灰度级（譬如白色），而对边界外部的所有像素分配另一个灰度（譬如黑色）。分割后的区域满足上面的条件。
下面三幅图说明了基于区域的分割。图(d)类似于图(a)，但内部区域的灰度形成了一幅纹理模式。图(e)显示了计算图中灰度不连续的结果。灰度中的大量寄生变化使得我们难以识别原图像中的唯一边界，因为很多非零灰度变化连接到了边界，所以基于边缘的分割不是一种合适的方法。然而，我们注意到，外部区域是恒定的，因此在解决这个分割问题时，只需一个能够区分纹理区域和恒定区域之间的不同的谓词逻辑。像素值的标准差是实现这一目的一种测度，因为纹理区域的标准差非零，而其他区域的标准差为零。图(f)显示了将原图像分成多幅大小为8x8的子区域后的结果。若子区域标准差为正，将这个子区域标记为白色。结果是，该区域边缘的周围成“块状”外观，因为8x8大小的一组方形都被标记为相同的灰度。
1.2 导数
一阶导数和二阶导数尤其适用于检测灰度的局部突变，数学函数的导数可根据有限差分来定义。
首先，根据泰勒级数，我们得到一维函数f(x)在任意点x处的一阶导数的近似：
一般来说，用于表示导数的泰勒级数的项越多，近似就越精确。包含更多的项意味着在近似中使用更多的点，产生更小的误差。但是对于一阶导数，我们只是用线性项。当dx=1时，上式为
当dx=-1时，
我们可以用前项差分得到灰度差：
反向差分
中心差分
对于相同数量的点，中心差分的误差较小。因此，导数通常表示为中心差分。
基于中心差分的二阶导数等于
为了得到三阶中心导数，我们需要在x的两边多加一个点。也就是说，我们需要f(x+2)和f(x-2)的泰勒级数展开，分别令dx=2和dx=-2得到，
类似地，四阶有限差分为
下图显示一个简化的剖面线，从左往右分别遇到了点、线和边缘时的特性，其中斜坡横跨4个像素，噪声点是1个像素，线是3个像素，台阶边缘的过渡发生在相邻像素之间。为简单起见，灰度级的数量限制为8。
我们可以发现：
（1）一阶导数通常产生粗边缘；（一阶导数在整个斜坡上不为零，而二阶导数在斜坡开始和结尾处不为零）
（2）二阶导数对精细细节（如细线、孤立点和噪声）有更强的响应；（二阶导数在增强急剧变化方面更为激进，增强细节（包括噪声）方面，二阶要好）
（3）二阶导数在灰度斜坡和台阶过渡处会产生双边缘响应；（二阶导数进入边缘和离开边缘符号相反，“双边缘”效应可用来定位边缘）
（4）二阶导数的符号可用于确定边缘的过渡是从亮到暗还是暗到亮。（亮到暗：负；暗到亮：正）
计算图像中每个像素位置的一阶导数和二阶导数的方法是空间卷积。通过计算核系数与核覆盖区域的灰度值的乘积之和得到滤波器在核的中心点的响应。
2 孤立点的检测
根据前一节得出的结论，我们知道点检测应以二阶导数为基础：
式中，偏导数用二阶有限差分计算得到：
上述表达式可用拉普拉斯核实现。当滤波器在这一点的响应的绝对值超过一个规定的阈值，则我们说在核的中心位置(x,y)检测到了一个点。在输出图像中，将这样的点标注为1，将其他点标注为0，就得到了一幅二值图像。换句话说，输出表示如下：
其中g(x,y)是输出图像，T是一个非负阈值。
示例：图像中孤立点的检测
下面是喷气发动机涡轮叶片的x射线图像。图像右上象限的叶片中有一个由单个黑色像素表示的通孔。利用拉普拉斯核滤波后，取T等于图像像素最高绝对值的90%。最高检测过程相当特殊，因为它基于单个像素位置的灰度突变。不满足这个条件时，其他方法更适合检测灰度变化。
Mat src = imread("./2.tif", 0);
Mat tmp, dst1;
Mat dst = Mat::zeros(src.size(), CV_8UC1);
// 利用拉普拉斯核滤波
cv::Laplacian(src, tmp, CV_32FC1);
// 取最高绝对值的90% 进行分割
double minValue, maxValue;
Point minIdx, maxIdx;
minMaxLoc(tmp, &minValue, &maxValue, &minIdx, &maxIdx);
cout << maxValue << endl;
for (int i = 0; i < src.rows; i++) {
for (int j = 0; j < src.cols; j++) {
if (tmp.at<float>(i, j) > maxValue*0.9)
dst.at<uchar>(i, j) = 255;
}
}
3 线检测
线检测的复杂度更高，二阶导数将导致更强的滤波器响应，并且产生比一阶导数更细的线。于是对于线检测我们也可以使用拉普拉斯核，但要记住必须正确处理二阶导数的双线效应。
示例：用拉普拉斯核进行线检测
下图为电子线路图模板的一部分。拉普拉斯图像包含负值，可以通过取绝对值来简单处理。然而这种方法会使线的宽度加倍。更合适的方法是只是用拉普拉斯图像的正值（在噪声情况下，我们使用超过阈值的那些值来消除由噪声导致的关于零的随机变化）。当线的宽度大于拉普拉斯核的大小时，这些线会被一个零值“山谷”分开。这是预料之中的，例如把3×3核居中放到一条宽为5像素的恒定灰度线上时，响应为零，于是产生了刚才提到的效应。当我们谈论线检测时，假设这些线要细于检测子的大小。我们最好把不满足这一假设的线当做区域，并用下一节中讨论的边缘检测方法来处理。
Mat src = imread("./2.tif", 0);
Mat tmp, dst1;
Mat dst2 = Mat::zeros(src.size(), CV_8UC1);
cv::Laplacian(src, tmp, CV_32FC1);
//取绝对值
convertScaleAbs(tmp, dst1);
//取正值
for (int i = 0; i < src.rows; i++) {
for (int j = 0; j < src.cols; j++) {
if (tmp.at<float>(i, j) > 250)
dst2.at<uchar>(i, j) = 255;
}
}
上面的拉普拉斯核是各向同性的，因此其响应和方向无关。通常，我们的兴趣在于检测规定方向的线，下图中的四个核，分别对水平、45°、垂直和-45°方向线的响应最好。每个核的首选方向用一个比其他方向更大的系数加权。每个核中的系数之和为零，这表明恒定灰度区域中的响应为零。
示例：寻找宽度为1像素，方向为水平/45°方向的线
Mat src = imread("./3.tif", 0);
Mat tmp, dst;
Mat k1 = (Mat_<double>(3, 3) << -1, -1, -1, 2, 2, 2, -1, -1, -1);
Mat k2 = (Mat_<double>(3, 3) << 2, -1, -1, -1, 2, -1, -1, -1, 2);
filter2D(src, tmp, CV_32FC1, k1);
normalize(tmp, tmp, 0, 255, NORM_MINMAX);
tmp.convertTo(tmp, CV_8UC1);
threshold(tmp, dst, 254, 255, THRESH_BINARY);
4 边缘检测
4.1 边缘模型
边缘模型可以根据它们的灰度剖面来分类。在实际工作中，数字图像都存在模糊且带有噪声的边缘，模糊程度主要取决于聚焦机制（如镜头）的限制，而噪声水平主要取决于成像系统的电子元件。在这些情况下，边缘被建模为一个更接近灰度斜坡的剖面线，如下图中间所示。斜坡的斜度与边缘的模糊程度成反比。在这一模型中，沿剖面线不再有单个“边缘点”，相反，边缘点是现在斜坡中包含的任何点，并且边缘线段是一组已连接起来的这样的点。
下面为实际工作中常见的斜坡边缘模型。沿灰度剖面线从左到右移动时，我们发现在斜坡上一阶导数为正，在恒定灰度区域的一阶导数为零。在斜坡的开始处，二阶导数为正；在斜坡上二阶导数为零。零灰度轴和二阶导数极值间的连线的交点，称为该二阶导数的零交叉点或过零点。由这些观察可以得出结论：
（1）一阶导数的幅度可以用于检测图像中的某个点处是否存在一个边缘；
（2）二阶导数的符号可用于确定一个边缘像素是位于边缘的暗侧还是亮侧；
（3）对图像中的每个边缘，二阶导数生成两个值；
（4）二阶导数的过零点可以用于确定粗边缘的中心位置
通常用于边缘检测的三个步骤如下：
1. 为了降低噪声，对图像进行平滑处理。微弱的噪声严重影响检测边缘所用的两个关键导数。
2. 检测边缘点。这是一个局部运算，它从图像中提取可能是边缘点的所有点（候选边缘点）。
3. 边缘定位。从候选边缘点中选择组成边缘的点。
4.2 基本边缘检测
4.2.1 图像梯度及其性质
求图像f中任意位置(x,y)处的边缘强度和方向的工具是梯度，梯度向量为
梯度向量在点(x,y)处的幅度M(x,y)由欧几里得向量范数给出：
梯度向量指出了f在(x,y)处的最大变化率的方向：
角度是相对于x周逆时针方向度量的。g(x,y)，M(x,y)，α(x,y)都是与f大小相同的图像。
示例：计算梯度
下面是包含边缘线段的一幅图像的部分放大。每个方格对应一个像素，我们需要计算框出点的边缘强度和方向。图中阴影像素假定为0，白色像素的值假定为1。以3×3邻域来计算x方向和y方向的导数，从底部一行的像素减去顶部一行邻域中的像素，得到x方向的偏导数；从右边邻域中的像素减去左列的像素得到y方向的偏导数。接下来，用这些差值作为计算的偏导数，在当前点处有，
由此得到梯度的幅值M=2根号2，梯度方向为-45度，边缘方向为α-90°=135°-90°=45°。梯度向量有时也称为边缘法线。当向量除以其幅值而归一化到单位长度时，结果向量通常称为边缘单位法线。
4.2.2 梯度算子
要得到一幅图像的梯度，就要在图像的每个像素位置计算偏导数。对于梯度，我们常使用前向或中心有限差分。使用前向差分得
利用下面的一维核进行滤波，可以实现这两个公式：
当我们感兴趣的是对角边缘方向时，需要一个二维核。罗伯特交叉梯度算子是最早使用具有对角性能的二维核的算子之一考虑下图的3×3区域，
罗伯特算子利用gx=z9-z5和gy=z8-z6，实现对角差分。2×2大小的核概念上很简单，但在计算边缘方向时，它们不如关于中心对称的核有用，中心对称核的最小尺寸为3×3。使用3×3的核时，对偏导数的最简单的数字近似为gx=(z7+2z8+z9)-(z1+2z2+z3)和gy=(z3+2z6+z9)-(z1+2z4+z7)，Prewitt算子可以实现该计算。Sobel算子在中心位置上使用了权值2，可以平滑图像，更好地抑制噪声。注意，所有核的系数之和为零，即对恒定灰度区域的响应为零。
由于平方和平方根的计算开销较大，因此我们常用绝对值来近似梯度幅度：
上式不仅计算简单，而且仍然能够保留灰度级的相对变化。但代价是，滤波器通常不再是各向同性（旋转不变）的。这对Prewitt和Sobel算子没有影响，因为这些核只为垂直边缘和水平边缘给出各向同性的结果。
下图为Kirsch罗盘核设计用于检测所有8个罗盘方向的边缘幅度和方向。Kirsch不在用上面的方法计算幅度和角度，此时的边缘幅值为该点处给出最强卷积值的核的响应，边缘角度为与核相关联的方向。比如，图像中某点处的最强响应是由北（N）核产生的，则该点的边缘幅值为这个核的响应，方向为0°
示例：二维梯度
Mat src = imread("./4.tif", 0);
Mat gx, gy, gxy;
Sobel(src, gx, CV_16S, 1, 0, 3);
Sobel(src, gy, CV_16S, 0, 1, 3);
convertScaleAbs(gx, gx);
convertScaleAbs(gy, gy);
gxy = gx + gy;
示例：梯度+平滑+阈值处理
上面墙砖对图像细节的贡献非常显著，这种细节在边缘中通常是我们要避免的，因为它往往表现为噪声，而导数计算会增强这种噪声，使得主要边缘的检测变得复杂。在计算梯度之前对图像进行平滑，可使得边缘检测更有选择性。实现同一目标的另一种方法是，对梯度图像进行阈值处理。实践中通常将两者结合，突出边缘的同时尽可能保留许多连通性。
Mat src = imread("./4.tif", 0);
blur(src, src, Size(5,5));
...(同上）
threshold(gx, gx, 120, 255, THRESH_BINARY);
threshold(gy, gy, 120, 255, THRESH_BINARY);
gxy = gx + gy;
4.3 更先进的边缘检测技术
4.3.1 Marr-Hildreth边缘检测子
Marr和Hildreth认为：（1）边缘检测中灰度变化和图像的尺寸无关，检测算子可以为不同尺度，大算子检测模糊边缘，小算子检测清晰的细节；（2）灰度变化梯度在一阶导数的极值点，或在二阶导数为零的交叉点。
标准差为σ的二维高斯函数，
代入拉普拉斯，得到
这个表达式称为高斯拉普拉斯（LoG）函数。
Marr和Hildreth算法如下：
（1）首先让LoG核与一幅输入图像卷积，
（2）然后寻找g(x,y)的过零点来确定f(x,y)中的边缘的位置。任意像素的过零点意味着，其左右，上下，两个对角的像素符号不同，且差值以及在这一点的值大于某个阈值。
因为拉普拉斯变换和卷积都是线性运算，可以将上式写为
即，先用高斯滤波器平滑图片（使用大于6σ的最小奇整数），然后计算结果的拉普拉斯变换。
示例：Marr-Hildreth边缘检测子
#include <iostream>
#include <opencv2/opencv.hpp>
using namespace std;
using namespace cv;
void MarrHildrethEdge(const Mat& srcImage, Mat& dstImage, int kSize, double sigma, int T)
{
//1. 计算卷积核
int kernel_dia = kSize * 2 + 1;
Mat kernel(kernel_dia, kernel_dia, CV_64FC1);
for (int i = -kSize; i <= kSize; i++)
{
for (int j = -kSize; j <= kSize; j++)
{
kernel.at<double>(i + kSize, j + kSize) = exp(-(pow(i, 2) + pow(j, 2)) / (pow(sigma, 2) * 2))
* ((pow(i, 2) + pow(j, 2) - 2 * pow(sigma, 2)) / (2 * pow(sigma, 4)));
}
}
//2. 卷积
Mat laplaceMat(srcImage.size(), CV_64FC1);
dstImage = Mat::zeros(srcImage.size(), CV_8UC1);
for (int i = kSize; i < srcImage.rows - kSize; i++)
{
for (int j = kSize; j < srcImage.rows - kSize; j++)
{
double sum = 0;
for (int x = -kSize; x <= kSize; x++)
{
for (int y = -kSize; y <= kSize; y++)
{
sum += srcImage.at<uchar>(i + x, j + y) * kernel.at<double>(x + kSize, y + kSize);
}
}
laplaceMat.at<float>(i - kSize, j - kSize) = sum;
}
}
//3. 查找过零点
for (int i = kSize; i < dstImage.rows - kSize; i++)
{
for (int j = kSize; j < dstImage.cols - kSize; j++)
{
float g1 = laplaceMat.at<float>(i - 1, j - 1);
float g2 = laplaceMat.at<float>(i - 1, j);
float g3 = laplaceMat.at<float>(i - 1, j + 1);
float g4 = laplaceMat.at<float>(i, j - 1);
float g5 = laplaceMat.at<float>(i, j);
float g6 = laplaceMat.at<float>(i, j + 1);
float g7 = laplaceMat.at<float>(i + 1, j - 1);
float g8 = laplaceMat.at<float>(i + 1, j);
float g9 = laplaceMat.at<float>(i + 1, j + 1);
if (g5 < -T)
if ((g1 * g9 < 0 && abs(abs(g1) - abs(g9)) > T) || (g2 * g8 < 0 && abs(abs(g2) - abs(g8)) > T) ||
(g3 * g7 < 0 && abs(abs(g3) - abs(g7)) > T) || (g4 * g6 < 0 && abs(abs(g4) - abs(g6)) > T))
dstImage.at<uchar>(i, j) = 255;
}
}
}
int main(int argc, char* argv[])
{
Mat srcImage = imread("./5.tif",0);
Mat dstImage;
MarrHildrethEdge(srcImage, dstImage, 13, 2, 30);
imshow("srcImage", srcImage);
imshow("dstImage", dstImage);
waitKey();
return 0;
}
4.3.2 Canny边缘检测子
虽然算法更复杂，但本节讨论的canny边缘检测子时迄今为止讨论的边缘检测算子中最优秀的。canny方法的步骤如下：
1. 使用圆形二维高斯函数平滑图像；
2. 计算梯度幅度和角度图像；
3. 利用非极大值抑制，细化梯度图像在局部极大值附近包含的一些宽脊；
4. 使用双阈值处理和连通性分析来检测与连接边缘，减少假边缘点。
具体算法与实现在下篇博客《Canny边缘检测实现（Opencv C++）》
4.4 连接边缘点
理想情况下，边缘检测应产生只位于边缘上的像素集合。在实际应用中，由于存在噪声、非均匀光照引起的边缘断裂及导致灰度值不连续的其他效应，因此这些像素很少能够完全表征边缘。因此，边缘检测后，通常要执行连接算法，将边缘像素组合为有意义的边缘和区域边界。本节讨论两种连接边缘的基本方法：第一种方法需要关于局部区域（如3×3邻域）中边缘点的知识；第二种方法是全局方法，它处理的是整个边缘图形。
4.4.1 局部处理
连接边缘点的一种简单方法是，在（已声明为边缘点的）每一个点(x,y)的一个小邻域内分析像素的特点。根据预定义的准则，将所有相似地点连接起来，形成具有相同性质的像素的一个边缘。假设以边缘点(x,y)为中心的邻域坐标集为Sxy，在这种局部分析中，用于建立边缘像素相似性的两个主要性质是
（1）梯度向量的幅度
如果满足上式，那么Sxy中坐标为(s,t)的边缘像素在幅度上与(x,y)处的像素是相似的。
（2）梯度向量的方向
如果满足上式，那么Sxy中坐标为(s,t)的边缘像素的角度与(x,y)处的像素是相似的。如果既满足幅度准则，又满足方向准则，那么Sxy中的坐标(s,t)的像素连接到(x,y)处的像素。对每个边缘像素重复这一过程。当邻域的中心逐个像素移动时，必须记录那些已连接的点。简单的记录过程是对每组已连接的像素分配不同的灰度值。
上述公式的计算开销很大，因为必须检查每个点的所有邻域。简化步骤如下：
1.计算梯度幅度和角度
2.形成二值图像
其中，TM是阈值，A是规定的角度方向，±TA定义了关于A的可接受方向的一个“条带”。
3.对二值图像g，逐行扫描并填充水平间隙。
4.对二值图像g，逐列扫描并填充垂直间隙。编程时可以将g转置后，再做一次逐行扫描水平填充来实现。
在必要时，还可以对二者图像 g 进一步在其它任何方向上扫描并填充间隙。
示例：局部处理（未填充）
Mat src = imread("./7.PNG", 0);
// 计算待测物体旋转角度
Mat thrImg;
threshold(src, thrImg, 128, 255, THRESH_BINARY_INV);
vector<vector<Point>> cons;
findContours(thrImg, cons, RETR_LIST, CHAIN_APPROX_SIMPLE);
int max_idx = 0; double max_area = 0.0;
for (int i = 0; i < cons.size(); i++) {
double area = contourArea(cons[i]);
if (area > max_area) {
max_idx = i;
max_area = area;
}
}
RotatedRect rect = minAreaRect(cons[max_idx]);
// 计算梯度幅度和角度
Mat grad_x, grad_y;
Sobel(src, grad_x, CV_32F, 1, 0, 3);
Sobel(src, grad_y, CV_32F, 0, 1, 3);
Mat directImg, gradXY;
divide(grad_y, grad_x, directImg);
gradXY = abs(grad_x) + abs(grad_y);
double minValue, maxValue;
Point minIdx, maxIdx;
minMaxLoc(gradXY, &minValue, &maxValue, &minIdx, &maxIdx);
// 二值化处理
double TM = 0.25 * maxValue;
float A = rect.angle; int TA = 5;
Mat edgeX = Mat::zeros(src.size(), CV_8UC1);
Mat edgeY = Mat::zeros(src.size(), CV_8UC1);
for (int i = 0; i < grad_x.rows; i++) {
for (int j = 0; j < grad_x.cols; j++) {
double angle = atan(directImg.at<float>(i, j)) / CV_PI * 180;
if ((gradXY.at<float>(i, j) > TM) && (angle > A - 90 - TA && angle < A - 90 + TA))
edgeX.at<uchar>(i, j) = 255;
if ((gradXY.at<float>(i, j) > TM) && (angle < A + TA && angle > A - TA))
edgeY.at<uchar>(i, j) = 255;
}
}
Mat dst;
bitwise_or(edgeX, edgeY, dst);
4.4.2 全局处理（霍夫变换）
上一节中讨论的方法适用于已知各个目标的像素的情况。通常，我们必须工作于非结构化环境中，这时我们所拥有的只是一幅边缘图，不知道感兴趣目标位于何处。在这种情况下，所有像素都是连接的候选点，因此必须根据预先定义的全局性质接受或删除某些点。本节介绍一种基于像素集是否在规定形状的曲线上的方法。一旦检测到这些曲线，它们就会形成边缘或感兴趣的区域边界。
已知一幅图像中的n个点，假设我们希望找到这些点中位于直线上的子集。一种可能的解决方法是，首先找到由每对点确定的所有直线，然后寻找靠近特定直线的那些点的所有子集。这种方法涉及寻找n(n-1)/2~n2条直线，然后将每个点与所有直线执行n(n(n-1))/2~n3次比较。在大多数应用中，这都是一项困难的计算任务。
Hough提出了一种替代方法，通常称为霍夫变换。具体概念与用法见链接。
参考：
1. 冈萨雷斯《数字图像处理（第四版）》Chapter 9 （所有图片可在链接中下载）
3. 【youcans 的 OpenCV 例程200篇】156. 边缘连接局部处理的简化算法
