OpenCV-8-图像分析

1 摘要

上一章节介绍了OpenCV中可用的图像变换函数，这些技术本质上都是通过一种映射关系将图像转换为另一个图像，这种转换会保留输入图像的主要特征。本章将要介绍的图像处理技术则可能得到和输入图像完全不一样的结果，甚至得到的矩阵中元素的含义也会发生根本变化。如下文将会将到的离散傅立叶变换(Discrete Fourier Transform，DFT)处理图像后的到的是像素的颜色频域表示，而原始图像是像素的颜色强度表示。某些变换得到的结果甚至已经不再是矩阵，而是类似成分列表，如霍夫线变换(Hough Line Transform)。

在本章末尾还会介绍图像分割技术，用于分割图像中有意义的相连接区块。

2 频域-空域变换

2.1 离散傅立叶变换

傅立叶变换用于将连续的时域信号转换为连续的频域信号，而离散傅立叶变换用于处理离散的数据，如二维图像中的离散像素点。对于一维N个复数信号x0，x1…xn-1，其离散傅立叶变换公式定义如下。完整的傅立叶变换推导公式请参考文章数字信号处理

其中

二维的离散傅立叶变换公式定义如下，对于更高维度的傅立叶变换这里暂时不讨论。

从上述公式中可以看出，对于N个样本的信号集，得到其N个频域表示分量离散傅立叶变换需要的计算负责度为O(N2)，实际上快速傅立叶变换(Fast Fourier Trasform, FFT)能够将算法复杂度降至O(N logN)。

OpenCV提供的离散傅立叶变换函数原型如下。该函数使用了快速傅立叶变换算法，能够处理1维和2维矩阵信号。在处理2维矩阵的时候可以选择将矩阵的每一行都看作是一个1维信号集，分别计算它们的1维傅立叶变换，这种处理方式的效率会比多次调用该函数更高。

// src：输入矩阵，待处理的信号
// dst：计算结果
// flags：计算策略，下文介绍
// nonzeroRows：有效的数据行，下文介绍
void cv::dft(cv::InputArray src, cv::OutputArray dst,
             int flags = 0, int nonzeroRows = 0);

输入矩阵必须是单通道或者双通道的浮点型数据。如果输入矩阵为单通道矩阵，则输入信号被认为是实数，而输出矩阵则被包装为一种特别的节约空间的形式，称为复数共轭对(Complex Conjugate Symmetrical, CSS)。此时输出矩阵的尺寸和输入矩阵的尺寸相同，因为对于离散傅立叶变换而言，根据尼奎斯特采样定理，得到的有效频率数为N/2，因此这种格式正好利用了这部分空间。

而对于双通道的输入矩阵而言，两个通道会分别被当作是复数的实部和虚部，数据不会被打包成CSS格式，并且其输出矩阵和输入矩阵尺寸相同，但是这部分空间将会被0填充。使用双通道矩阵处理数字图像时，第二个通道的复数虚部必须设置为0，一个简单的方式是使用函数cv::Mat::zeros()创建一个全0矩阵，然后调用函数cv::merge()和已有的像素颜色强度矩阵合并得到完整的复数矩阵，然后再调用函数cv::dft()。

CCS格式包装的1维信号矩阵傅立叶变换结果表示如下。

其中Re表示的是复数的实数部分，而Im表示的是复数的虚数部分。注意观察元素下标，特定元素缺失虚部数据，从而确保这部分数据计算结果一定为实数。另外只有当N为偶数时，最后一列才存在，否则将会被0填充。

2维信号矩阵傅立叶变换结果的CCS格式包装如下图。需要注意虚线将整个矩阵分为4个区域，左下角的行索引表示规则是处理真正的二维矩阵信号的示意图，而右下角的行索引表示规则是处理Ny个1维矩阵信号的示意图。同样的对于二维信号矩阵而言，当Nx不为偶数时，最后1列不会存在，当Ny不为偶数时，最后1行不会存在。

参数flags可以指定函数内部的算法策略，OpenCV提供了多个可选值，可以通过逻辑与符号组合这些选项。默认情况下执行的是离散傅立叶变换，而通过设置标记cv::DFT_INVERSE可以执行离散傅立叶逆变换。离散傅立叶逆变换和离散傅立叶变换的定义类似，它们的区别仅在于一个缩放系数和和指数的符号。需要注意如果想要执行离散傅立叶逆变换，你的输入矩阵应当符合前文描述的正变换输出结果的规则。

逆变换的缩放系数可以通过cv::DFT_SCALE启用，对于一维变换而言，这个系数为N^-1，对于二维变换而言，这个系数为(Nx Ny)^-1。如果是对一个正变换的结果应用逆变换，这个选项必须开启，否则将无法重建原始数据。由于标记cv::DFT_INVERSE和cv::DFT_SCALE通常会合并使用，OpenCV将它们组合为一个特殊标记cv::DFT_INVERSE_SCALE，可简写为cv::DFT_INV_SCALE。另外该参数的可选标记还有cv::DFT_ROWS，表示需要将二维的输入矩阵当作是Ny行一维矩阵处理。此外通过该标记可以处理三维及更高维度的矩阵，这里不再详细介绍。

前文说到过该函数在处理单通道矩阵时会使用CCS格式压缩输出矩阵尺寸，通过标记cv::DFT_COMPLEX_OUTPUT可以强制该函数输出复数结果。相反的，对一个复数矩阵应用逆变换得到的结果也是复数矩阵，如果输入矩阵具有复数对称性(这里不是说CCS格式，而是指对称性，如实数组的正向变换就有这种对称性)，则可以通过设置标记cv::DFT_REAL_OUTPUT强制输出实数矩阵。

离散傅立叶变换强烈倾向于处理特定长度的输入向量，在大多数DFT的算法实现中，希望的向量长度都是2的指数。在OpenCV中，偏好的向量长度或者矩阵维度维2、3或者5的指数，因此一个实用技巧是创建一个比实际数据更大的矩阵，多余的部分使用0填充。函数cv::getOptimalDFTSize()可以返回一个不小于输入尺寸的适用于离散傅立叶变换的最佳尺寸。此时通过参数nonzero_rows显示指明那些行是用于填充数据对其尺寸的，该函数在计算时会针对这些行进行优化提升运算效率。

2.2 离散傅立叶逆变换

除了使用函数cv::dft()及特定的参数执行离散傅立叶变换外，OpenCV还额外提供函数来执行离散傅立叶变换，它们是等效的，其函数原型如下。

// src：输入矩阵
// dst：输出矩阵
// flags：算法策略
// nonzeroRows：有意义的数据行数
void cv::idft(cv::InputArray src, cv::OutputArray dst,
              int flags = 0, int nonzeroRows = 0);

2.3 频谱乘法

在很多计算离散傅立叶变化的程序中，还需要逐元素的计算两个矩阵的乘积。由于傅立叶变换的结果是复数，并且可以是以CCS格式打包的，因此处理这部分数据可能比较繁琐。OpenCV提供专门的函数用于执行频频的乘法运算，其函数原型如下。

// src1：第一个频域信号矩阵，可以是单通道CCS格式或者双通道复数
// src2：第二个频域信号矩阵，可以是单通道CCS格式或者双通道复数
// dst：计算结果
// flags：计算策略，仅支持cv::DFT_ROWS表示输入矩阵的每一行为1维傅立叶变换的结果，
//        设置为0时表示2维傅立叶变换的结果
// conj：是否对第二个矩阵取共轭后再执行乘法运算，false用于卷积，而true用于计算相关性，
//       在后面的文章中会介绍
void cv::mulSpectrums(cv::InputArray src1, cv::InputArray src2,
                      cv::OutputArray dst, int flags, bool conj = false);

2.4 使用离散傅立叶变换加速卷积计算

卷积定理可以将空域上的卷积运算转换到频域上的乘法运算，从而极大的加快其算法计算效率，而函数从空域相频域的转换可以由傅立叶变换完成，对于离散的图像信号可以使用离散傅立叶变换。这里不再做消息的卷积定理公式推断，记住这个结论即可。示例DFTConvolution使用这种方式执行卷积运算，这部分代码直接拷贝自OpenCV的参考资料。

int main(int argc, const char * argv[]) {
    // 读取图像
    cv::Mat A = cv::imread(argv[1], 0);

    cv::Size patchSize(100, 100);
    cv::Point topleft(A.cols / 2, A.rows /2);
    cv::Rect roi(topleft.x, topleft.y, patchSize.width, patchSize.height);
    // 选中其中的一小块图像，需要注意实际的数据段仍然是共用的
    cv::Mat B = A(roi);
    
    // 获取离散傅立叶变换的最佳向量长度
    int dft_M = cv::getOptimalDFTSize(A.rows + B.rows - 1);
    int dft_N = cv::getOptimalDFTSize(A.cols + B.cols - 1);
    // 使用最佳尺寸创建待处理的样本矩阵
    cv::Mat dft_A = cv::Mat::zeros(dft_M, dft_N, CV_32F);
    cv::Mat dft_B = cv::Mat::zeros(dft_M, dft_N, CV_32F);
    // 映射内部实际有效数据矩阵，共用数据段内存
    cv::Mat dft_A_part = dft_A(cv::Rect(0, 0, A.cols, A.rows));
    cv::Mat dft_B_part = dft_B(cv::Rect(0, 0, B.cols, B.rows));
    // 将原图A和选择的区块B映射到共用的数据段内存，即此时dft_A和dft_B被填充对应数据
    // 此次平移各自矩阵中的均值单位
    A.convertTo(dft_A_part, dft_A_part.type(), 1, -mean(A)[0]);
    B.convertTo(dft_B_part, dft_B_part.type(), 1, -mean(B)[0]);
    
    // 执行傅立叶变换
    cv::dft(dft_A, dft_A, 0, A.rows);
    // 需要注意dft_B在执行傅立叶变换之前有效的数据数为100*100，在执行完成后整个矩阵都包含
    // 有效数据，推测内部对原始的频域值做了插值运算，以便接下来方便计算卷积
    cv::dft(dft_B, dft_B, 0, B.rows);
    
    // 将最后一个参数设置伟false则计算卷积，否则计算相关性
    cv::mulSpectrums(dft_A, dft_B, dft_A, 0, true);
    // 此处卷积并未使用前文常见的任何滤波器，如高斯滤波器，暂不讨论具体含义
//    cv::mulSpectrums(dft_A, dft_B, dft_A, 0, false);
    cv::idft(dft_A, dft_A, cv::DFT_SCALE, A.rows + B.rows - 1);

    // 获取相关性矩阵的有效数据区域
    cv::Mat corr = dft_A(cv::Rect(0, 0, A.cols + B.cols - 1, A.rows + B.rows - 1));
    // 标准话矩阵原始
    cv::normalize(corr, corr, 0, 1, cv::NORM_MINMAX, corr.type());
    // 取三次幂，提高对比度便于查看
    cv::pow(corr, 3.0, corr);

    // 对分割出的小块区域执行异或运算，等效于黑白反色
    B ^= cv::Scalar::all(255);
    
    // 显示处理结果
    cv::imshow("Image", A);
    // 显示相关性，相关性在后面章节中会详细介绍，新增可以简单的理解为函数的相似程度，
    // 即图像的相似程度
    cv::imshow("Correlation", corr);
    
    return 0;
}

使用傅立叶变换实现卷积运算的算法中，最慢的部分是执行变换的过程，对于一个N✖️N的图像而言，傅立叶变换的算法时间复杂度为O(N²logN)。假定使用尺寸为M✖️M的卷积核，并且M<N，则可以使用这个值近似估计整个算法的时间复杂度。相较于直接进行卷积运算的时间复杂度为O(N²✖️M²)而言，这种策略能明显提高算法的效率。需要注意的是如果使用的卷积核很小，则没有必要做这种转换。

2.5 离散余弦变换

离散余弦变换在处理数据时会将原始样本向x轴负方向对称扩充，对于实数域而言，计算扩充后的傅立叶离散变换根据三角函数的奇偶性质仅仅需要计算一半的样本即可，并且得到的频域矩阵都是有效的。离散余弦变换的具体公式推断这里不再详细阐述，如感兴趣可以点击此处。其定义如下。

尽管此处仅展示了一维离散余弦变换的公式，但是更高维度的离散余弦变换也存在，此处暂时不详细展开讨论。离散傅立叶变换的基本思想也适用于离散余弦变换，但是需要注意此时系数都是实数。OpenCV提供的离散余弦变换函数原型如下。

// src：待处理的空域样本数据
// dst：处理后的频域数据
// flag：算法策略，下文详细介绍
void cv::dct(cv::InputArray src, cv::OutputArray dst, int flags = 0);

除要求元素必须是实数，该函数的输入矩阵规则与函数cv::dft()类似，由于计算结果都为实数，输出矩阵也不再需要处理复数的打包。另外和函数cv::dft()不同，函数cv::dct()要求输入矩阵的元素个数必须为偶数，不足位可以用0补齐。当参数flags设置为cv::DCT_INVERSE表示执行离散余弦逆变换，指定为cv::DCT_ROWS可以将输入的矩阵每一行当成是一个一维随机变量样本进行处理，这两个值可以用逻辑与符号同时选择。另外由于正变换和逆变换都包含标准化相关操作，因此不提供类似DFT变换使用的flag参数cv::DFT_SCALE。

和函数cv::dft()一样，函数cv::dct()的效率也取决于矩阵的尺寸，其最佳尺寸可以通过如下公式计算，这里暂时不展开讨论原因，因为这涉及到OpenCV内部的算法实现。

size_t optimal_dct_size = 2 * cv::getOptimalDFTSize((N+1) / 2);

2.6 离散余弦逆变换

同样为了提高代码的可读性，OpenCV额外提供如下函数来执行离散余弦变换。该函数的调用结果和函数cv::dct()使用cv::DCT_INVERSE作为参数flag的值的调用结果相同，当然flags是否设置cv::DCT_ROWS也需要相同。

// src：待处理的空域样本数据
// dst：处理后的频域数据
// flag：算法策略，是否将二维矩阵看作多M行一维向量
void cv::idct(cv::InputArray src, cv::OutputArray dst, int flags = 0);

3 积分图

积分图(Integral Image)是一种用于快速计算指定区域像素强度和的数据结构，这种结构在很多场景中都非常有用，如Haar微波(Haar wavelets)计算，该技术用于人脸识别和类似算法。

OpenCV支持技术三种形式的积分图，分别是和(Sum)、平方和(Square Sum)以及倾斜和(Tilted Sum)积分图。每种积分图的尺寸都在原始图像尺寸的x和y轴上各扩展1行或者1列。

标准的求和积分图的计算公式如下。

平方和积分图计算公式如下。

倾斜求和积分图计算公式如下，这种方式等效于将原图旋转45度。

使用这三组公式可以计算任意直立或者倾斜的矩形区域的和、平均值和标准差。从而快速的执行模糊、近似梯度等相关运算，即使对可变窗口的情况而言，也可以快速的执行块相关运算。

例如对于由顶点(x1, y1)、(x2, y2)，其中x2 > x1, y2 > y1定义的矩形区域，其内部像素和可以由对应的积分图求出，计算公式如下。

积分图的通过逐行逐列的方式计算，即积分图内的每个元素都是通过已经计算好的元素计算得出，其计算公式如下。

例如对于一副7✖️5的灰度图像而言，其像素统计结果如下左图，其标准和积分图计算结果如下右图，你可以在下图中验证上述公式。如计算原图中由4个20强度像素值围成的矩形区域内像素强度和，可以通过积分图中对应位置的值计算，即398-9-10+1 = 380，这种求和方式对于计算任意大小的矩形区域时间复杂度都是O(1)。

3.1 标准求和积分图

标准求和积分图的函数原型如下。在处理基本数据类型为cv::F32的矩阵时，尽管也可以指定该格式作为计算结果，但是考虑到数据溢出的可能性，推荐使用cv::F64。

// image：待处理的图像，尺寸为M✖️N
// sum：计算出的积分图，尺寸为M+1✖️N+1
// sdepth：计算结果的基本数据类型，如cv::F32、cv::S32、cv::F64
void cv::integral(cv::InputArray image, cv::OutputArray sum, int sdepth = -1);

3.2 平方和积分图

平方和积分图的函数原型如下。

// image：待处理的图像，尺寸为M✖️N
// sum：计算出的标准和积分图，尺寸为M+1✖️N+1
// sqsum：计算出的平方和积分图，尺寸为M+1✖️N+1
// sdepth：计算结果的基本数据类型，如cv::F32、cv::S32、cv::F64
void cv::integral(cv::InputArray image,
                 cv::OutputArray sum, cv::OutputArray sqsum,
                 int sdepth = -1);

3.3 倾斜积分图

倾斜积分图的函数原型如下。

// image：待处理的图像，尺寸为M✖️N
// sum：计算出的标准和积分图，尺寸为M+1✖️N+1
// sqsum：计算出的平方和积分图，尺寸为M+1✖️N+1
// tilted：计算出的倾斜积分图，尺寸为M+1✖️N+1
// sdepth：计算结果的基本数据类型，如cv::F32、cv::S32、cv::F64
void cv::integral(cv::InputArray image,
                  cv::OutputArray sum, cv::OutputArray sqsum,
                  cv::OutputArray tilted, int sdepth = -1);

4 Canny边缘检测器

尽管可以直接使用例如拉普拉斯滤镜等简单的滤镜来检测图像的边缘，但是这种算法仍有提升的空间。J. Canny在1986年对这种算法进行了优化，发明了Canny边缘检测器(Canny Edge Detector)。Canny优化的方式是将计算得到的x轴和y轴上的一阶导数组合为4个方向上的导数。如果一个点在这四个方向上的导数是局部的最大值，则这些点就是边缘的候选点。该算法的最独特的创新点就是将单个候选边缘像素组装成为轮廓(Countours)。在后面的章节中会详细的对轮廓展开叙述，目前仅需要知道Canny边缘检测器不会直接返回检测目标的轮廓类型，如果需要确定轮廓类型，还需要使用Canny边缘检测器返回的结果调用函数cv::findContours()。

该算法对像素应用滞后阈值(Hysteresis Threshold)从而形成轮廓，也就是该算法使用了两个阈值，一个较大值和一个较小值。规则A，如果一个像素包含比较大值更大的梯度，则其会被认为是边缘像素。规则B，如果一个像素的所有梯度值都低于较小阈值则认为不是边缘像素。规则C，如果一个像素不包含比较大值更大的阈值，但是有比较小值更高的阈值，此时如果它和一个通过规则 A确定的边缘像素相邻，则它也会被认为是边缘像素，否则不是。Canny建议的最大和最小阈值的比值为2:1到3:1之间。

下图是对图像应用上下阈值粉笔为50和10，即比值为5:1的Canny边缘检测效果。

下图是对图像应用上下阈值粉笔为150和100，即比值为3:2的Canny边缘检测效果。

OpenCV定义的Canny边缘检测函数原型如下。

// image：待处理的图像矩阵
// edges：处理后的边缘图像
// threshold1：低梯度信号阈值
// threshold2：高梯度信号阈值
// aperturSize：函数内部依赖的Sobel梯度运算使用的孔径大小
// L2gradient：计算方向梯度使用的范数类型，下文详细介绍
void cv::Canny(cv::InputArray image, cv::OutputArray edges,
               double threshold1, double threshold2, int apertureSize = 3,
               bool L2gradient = false);

该函数的参数L2gradient可以指定计算方向梯度使用的范数类型，默认值为False，即使用L1范数，这种方式算法的效率更高，但是得到的方向梯度精确度更低，其计算公式如下。

当参数L2gradient设置为true时，即表示使用L2范数，其计算公式如下。

5 霍夫变换

霍夫变换(Hough Transform)可以寻找图像中的线、圆和其他简单形状。最初的霍夫变换是一个线性变换，能快速的寻找二值图像中的直线，该变换经过近一步推广后可以用于寻找其他的简单形状。

5.1 霍夫线变换

霍夫线变换的基本思想是对于一个二值图像中的任意非0点，原图都有可能存在各个方向的直线穿越该点。如果将这些线参数化表示，如使用斜率a和y轴截距b来表示这些直线，则过该点的这些直线在ab组成的坐标系中将形成一条新的轨迹。对于输入图像xy平面上的所有非零像素点处理完成后，在ab坐标系中具有局部最大值的点就是在原始图像中存在的曲线，因为原图中大量的点都有可能被这条直线穿越。

在实际实现的时候并不会选择使用斜率和截距表示一条直线，因为使用斜率表示的直线分布密度是非线性的，另外斜率的取-∞到+∞，也不利于表示。因此OpenCV使用极坐标系中的点来表示一条直线。如下a图中对于原始图像中的任意一点P(x0, y0)，在b图中可能存在1、2、3、4等直线过这一点，对于直线1而言，可以用极坐标系中的点Polor(p, θ)来表示唯一的这条直线，直线过该点，并且与该点至圆心的连线锤子。则在c图中这些直线就可以通过一系列的点形成的轨迹表示。

OpenCV的霍夫算法并不会直接将这些直线返回，相反它返回了使用极坐标表示的定义这些直线的点，如Line(p, θ)。OpenCV支持三种不同类型的霍夫变换，分别是标准霍夫变换(Standard Hough Transform, SHT)、多尺度霍夫变换(Multiscale Hough Transform, MHT)和渐进概率霍夫变换(Progressive Probabilities Hough Transform, PPHT)。

标准霍夫变换就是上文所讨论到的标准算法，多尺度霍夫变换是多标准霍夫变换的改进，得到的匹配直线精确度更高，概率渐进霍夫变换是标准霍夫变换的一个变体，它除了计算在直线参数累加平面内计算可能的直线外，还会计算直线的长度。这种算法被称为概率的原因是他不会在参数累加平面内统计所有的点，而是只会计算部分。它的基本思想是如果峰值足够高，则花更少的时间足以寻找到结果。这种方式能够大量的减少计算时间。

5.1.1 标准霍夫变换和多尺度霍夫变换

标准霍夫变换和多尺度霍夫变换共用一个函数，区别在于最后两个参数是否为0，如果是则为标准霍夫变换，否则为多尺度霍夫变换，其函数原型如下。

// image：待处理图像
// lines：N✖️1矩阵，表示所有寻找到的直线
// rho：参数累加平面内极坐标的幅度分辨率，单位为像素
// theta：参数累加平面内极坐标的弧度分辨率，单位为弧度
// threshold：非标准化的阈值，即在参数累加平面内某个点应当被认为定义了原图中的某条直线的阈值
// srn：多尺度霍夫变换对幅度分辨率的改进
// stn：多尺度霍夫变换对弧度分辨率的改进
void cv::HoughLines(cv::InputArray image, cv::OutputArray lines,
                    double rho, double theta, int threshold,
                    double srn = 0, double stn = 0);

该函数的输入图像的像素数据位深度必须为8，但是内部会将其处理为二值图像，即对于函数而言，所有分零值将会被看作是同一种情况。返回的参数lines矩阵元素数据类型为双通道浮点型，两个通道分别用于存储定义原图中直线极坐标的幅度和弧度分量。参数rho和theta定义了直线参数累加平面的幅度和弧度分辨率，该平面可以被看成是一个尺寸为rho✖️theta的二维直方图。参数threshold是一个未标准化的值，因此你需要根据实际处理图像的尺寸来调整该值的大小，该值可以理解为元素图像中至少应当有多少个像素点在同一条直线上，才认为这条直线是实际存在于原图中的。

参数srn和stn是专用于多尺度霍夫变换的，它们指定了更高的分辨率。多尺度霍夫变换先用参数rho和theta定义的平面计算原图可能存在的线的参考点极坐标，然后再使用参数src和stn缩放坐标，即最终的直线参数累加平面的分辨率为rho/src✖️theta/stn。这里暂时不需要理解其内部细节，只需要知道多尺度霍夫变换相较于标准霍夫变换具有更高的精度，但是成本更高即可。

5.1.2 渐进概率霍夫变换

渐进概率霍夫变换可以寻找线段，这里暂时不关注其内部的算法实现，OpenCV提供的函数原型如下。

// image：待处理的图像，单通道矩阵
// lines：N✖️1的4通道矩阵
// rho：参数累加平面内极坐标的幅度分辨率，单位为像素
// theta：参数累加平面内极坐标的弧度分辨率，单位为弧度
// threshold：非标准化的阈值，即在参数累加平面内某个点应当被认为定义了原图中的某条直线的阈值
// minLineLength：线段的最低长度，只有大于这个长度的线段才会被通缉
// maxLineGap：线段的最低间隔，小于该间隔的共线线段将会被连接成一个线段
void cv::HoughLinesP(cv::InputArray image, cv::OutputArray lines,
                     double rho, double theta, int threshold,
                     double minLineLength = 0, double maxLineGap = 0);

该函数的输出矩阵lines和标准霍夫变换以及多尺度霍夫变换函数不同，该函数的输出矩阵是4通道的，分别表示线段的两个端点(x0, y0)和(x1, y1)。

对一副格式图像和实际图像应用渐进概率霍夫变换的效果如下，可以明显看到渐进霍夫变换能够找到图像中的明显线段。其中渐进概率霍夫变换使用的参数minLineLength和maxLineGap分别是50和10，并且其使用的图像是先使用Canny边缘检查函数对原图像的处理结果，Canny边缘检测函数使用的两个阈值参数分别是50和150。

5.2 霍夫圆变换

将霍夫线变换推广到霍夫圆变换(Hough Circle Transform)，直接推广到方法是使用三维圆参数累加空间替代霍夫线变换里使用到的二维参数累加平面，因为在累加空间中表示一个原图中的唯一圆需要使用圆心的x和y坐标分量，以及园半径这三个参数。但是这种方式会带来昂贵的内存和算法计算时间成本，因此OpenCV使用霍夫梯度技术(Hough Gradient Method)来解决这个问题。

霍夫圆变换的效果如下图，在测试图像中找到了一个圆，在风景照发中未发现圆。需要注意这里首先对元素图像进行了Canny边缘检测。

霍夫梯度法的工作流程如下，首先对原图进行边缘检测，对于上述图像的示例使用的是Canny边缘检测。然后对于边缘检测结果中的每个非零像素找到其在原图中对应像素，使用Sobel梯度计算公式计算该点在x和y轴上的一阶梯度。然后沿着这个梯度表示的斜率表示的斜率在累加器遍历从一个最小值到最大值的点(这部分是直译，从非原书讲述的角度看，非零像素点的索贝尔梯度计算出的是器在该点的切线，切线的垂线应当经过圆心，因此将切线的垂线在累加平面内绘制，则对于圆心而言会存在极大值)，并将其值递增1，同时记录边缘检测结果中所有非0值像素的坐标。

当处理完边缘检测图像中的所有非零点后，在累积平面内那些比直接临近点更大的超过指定阈值点局部最大值点，将会被认为是圆心的候选点。这些候选点以值的大小降序排列，然后遍历这些候选点，如果将边缘检测图像中提取的非零点以其到候选圆心的距离进行排序。权衡最小距离和最大半径，选择一个能够囊括更多非零像素的圆的半径。如果圆囊括了更多的非零像素，并且与之前的圆心具有足够的距离，则保留该圆心。

这种实现方式加速了算法的效率，更重要的是克服了三维累加器的稀疏群体问题，这个问题会导致更多的噪声并使得结果不稳定。另外一方面，这个算法的一些缺点也应当受到重视。

首先使用Sobel导数计算局部梯度估计切线的方式并不是一个数值稳定的命题，会产生一定的噪声。其次每个候选中心都需要考虑边缘检测结果中的所有非零像素值，因此如果设置的阈值过低，算法会花很长的时间去处理这些候选圆心。再次，该算法无法检测同心圆，只会返回其中的一个，通常返回的是半径更大的那个。返回更大圆这种结果不是必然现象的原因是Sobel导数计算会出现噪声，对于分辨率无限大的图像而言，则必然会返回最大圆。

霍夫圆检测的函数原型如下。

// image：待检测的单通道图像
// circles：检测到的圆，N✖️1的三通道矩阵或者是Vec3f组成的向量
// method：只能选择cv::HOUGH_GRADIENT
// dp：累加图像的分辨率
// minDist：两个圆心的最小间距，低于此值会被认为是同一个圆
// param1：Canny边缘检测器使用的上阈值，下阈值将被设置为此值的1/2
// param2：未标准化的累加阈值，用于确定候选圆心
// minRadius：寻找到圆的最小搜索半径
// maxRadius：寻找到圆的最大搜索半径
void cv::HoughCircles(cv::InputArray image, cv::OutputArray circles,
                      int method, double dp, double minDist,
                      double param1 = 100, double param2 = 100,
                      int minRadius = 0, int maxRadius = 0);

霍夫圆检测函数和霍夫线检测函数的一个明显区别是，函数cv::HoughCircles()不再限制输入图像是二值图像，因为该函数会计算Sobel导数，因此你可以提供灰度图像。输出矩阵circles可能是数据类型为CV::F32C3的N✖️1的三通道矩阵或者是Vec3f组成的向量，每个元素的三个分量分别表示检测到的圆的圆心以及其半径。如果使用的是向量，则类型必须是std::vector<Vec3f>。参数dp是累加图像的分辨率，我们可以指定使用比原图更小的分辨率，该值必须大于等于1，最终分辨率为原始图像分辨率/db。

示例HoughCircles使用函数cv::HoughCircles()寻找图像中的圆，其核心代码如下。

int main(int argc, const char * argv[]) {
    // 读取原始图像
    cv::Mat src = cv::imread(argv[1], cv::IMREAD_COLOR);
    
    cv::Mat image;
    // 转换为灰度图像
    cv::cvtColor(src, image, cv::COLOR_BGR2GRAY);
    
    // 应用霍夫圆变换
    std::vector<cv::Vec3f> circles;
    cv::HoughCircles(image, circles, cv::HOUGH_GRADIENT, 2, image.cols/4);

    // 使用霍夫圆变换的结果在原图中绘制圆形
    for (size_t i = 0; i < circles.size(); ++i) {
        cv::circle(src,
                   cv::Point(cvRound(circles[i][0]), cvRound(circles[i][1])),
                   cvRound(circles[i][2]),
                   cv::Scalar(0, 0, 255, 1),
                   2, cv::LINE_AA);
    }
    
    // 显示霍夫圆变换的结果
    cv::imshow("Hough Circles", src);
    return 0;
}

该示例的运行结果如下图。

值的思考的是无论我们采取了什么策略，总不能绕过描述一个圆需要三个自由度(x, y, r)，描述一条线只需要两个自由度(p, θ)的事实，因此寻找图像中的圆总会消费更多的内存和时间。所以需要将半径参数限定在合理的范围来控制算法的开销，实际上该算法能够很好的找到圆心，但是并不能十分准确的找到半径。在一些特定的应用中，如果使用其他算法来寻找准确的半径或者只需要圆心数据，可以护绿该算法返回的这部分信息。Ballard在1981年将对象看作是渐变边缘的集合，从而将该算法推广到寻找任意形状。

6 距离变换

距离变换(Distance Transform)的结果是一副新的图像，图像中的每个点都等于该点在输入图像中对应像素距离零像素点最近距离。可以明显看出这种变换的典型输入图像应该是某一种类型的边缘图像。在大多数应用中输入图像是诸如Canny边缘检测器等边缘检测器的输出结果的反转值，即0值表示边缘，非零值表示非边缘。

计算距离变换的方法有两种，第一种方式使用了3✖️3或者5✖️5的矩阵蒙板。矩阵中的每个点的值被定义为是其到蒙板中心的距离。随着不断的移动这个蒙板，直至找到非零像素就能构建出一个当前像素距离最近非零像素的距离了，因此使用更大的蒙板会得到更精确的距离。当使用这种方式时，对于特定的距离计算方式，OpenCV需要选择合适的蒙板，这也是Borgefors在1986年最初发明的算法。第二种方式计算精确的距离，由Felzenszwalb提出，这两种算法的时间复杂度都和总像素数呈线性相关，但是精确算法会更慢一点。

6.1 无标记的距离变换

未标记的距离变换函数原型如下，

// src：待处理的图像
// dst：处理完成的图像，矩阵数据类型为cv::F32
// distanceType：计算距离使用的策略，下文介绍
// maskSize：蒙板尺寸，3、5或者是cv::DIST_MASK_PRECISE
void cv::distanceTransform(cv::InputArray src, cv::OutputArray dst,
                           int distanceType, int maskSize);

参数distanceType的可选值及其含义如下表。

distanceType的取值	含义
cv::DIST_C	棋盘距离，即x和轴的距离中的最大值
cv::DIST_L1	街区距离，即x和轴的距离和
cv::DIST_L2	欧式距离，即x和y轴的距离平方和的开方

当参数distanceType设置为cv::DIST_C和cv::DIST_L1时，参数maskSize设置为3就可以得到准确的结果，但是当参数distanceType设置为cv::DIST_L2时，Borgefors方法计算的总是近似距离，因此使用尺寸为5✖️5的蒙板会得到更接近实际值的欧式距离，但是计算成本会轻微增加。此外，当参数maskSize设置为cv::DIST_MASK_PRECISE时表示需要使用Felzenszwalb算法，但是只适用于计算cv::DIST_L2类型的距离。

6.2 有标记的距离变换

OpenCV提供的函数还可以在计算距离变换的同时，标记出每个像素的目标像素(即其最近零像素)。这些对象被称为时连接组件(Connected Components)，将在后面的章节中介绍，现在只需要将他们理解为连续的零像素组成的结构。其函数原型如下。

// src：待处理的图像
// dst：处理完成的图像
// labels：连接组件的标识
// distanceType：计算距离使用的策略
// maskSize：蒙板尺寸，3、5或者是cv::DIST_MASK_PRECISE
// labelType：标识策略
void cv::distanceTransform(cv::InputArray src, cv::OutputArray dst, 
                           cv::OutputArray labels,
                           int distanceType, int maskSize,
                           int labelType = cv::DIST_LABEL_CCOMP);

返回的标识矩阵labels尺寸和输入图像相同，其中的每个元素表示该点在输入图像中对应点的最近连接组件的标识，或者是最近零像素的坐标，这取决于下文将要将到的参数labelType的取值。得到的label矩阵也被称为离散Voronoi图，或者是泰森多边形，或是诺图。

参数labelType有两个可选值，其中cv::DIST_LABEL_CCOMP表示算法会自动将输入图像中的相连接的非零像素被划分为一个连接组件，并且为每个连接组件分配不同的标识。cv::DIST_LABEL_PIXEL表示会为输入图像中的每个非零像素分配不同的标识。

连通分量或者非零像素的标签是由该算法自动生成的，由于零像素的最近零像素就是自己，当参数labelType设置为cv::DIST_LABEL_PIXEL时，输出矩阵label中该点的值直接等于这个标识。当参数labelType设置为cv::DIST_LABEL_CCOMP时，该零像素属于连通分量的一部分，因此它和连通分量共享同一个标签。所以在任何情况下，如果想要知道原图中任意零像素的标签，只需要查询函数输出参数labels中对应像素的标签即可。

7 图像分割

图像分割是一个很大的话题，我们已经在多个地方接触到它，在后面的文章中会单独安排章节详细探讨。这里先聚焦几个OpenCV提供的方法，这些技术要么本身就是图像分割方法，要么它们处理的结果将会服务于后面章节将到的更复杂的图像学技术。需要注意即使目前也没有一个通用的绝对灵验的图像分割方法，这个话题仍然是计算机视觉领域活跃的研究方向。尽管如果，目前仍在一些特定的领域发展处理可靠的技术，实际应用中它们也能得到较好的结果。

7.1 浸水填充

浸水填充(Flood fill)在标记或者分离图像中的区块场景中十分有用，通常使用该算法得到的结果用于进一步的图像处理。浸水填充也可以从输入图像中计算得到一个特殊的蒙板，该蒙板可以用于加速后续函数的计算速度或者限制其只处理被该蒙板标识的区域。

相比于其他计算机绘图程序，OpenCV中的浸水填充是一种更一般化的填充方法。这些填充算法首先都会在图像中确定一个种子点(Seed Point)，然后使用同一种颜色为其以及其邻近的相似像素点着色。这些填充算法的区别是，OpenCV的浸水填充算法并不严格要求这些邻近像素的颜色完全相同。浸水填充的结果总会是一个连续的区域，即原图中的单个连通区块。如果某个像素和其相邻像素，或者是与之前选定的种子点，的颜色强度值在指定的偏差范围内，则该算法会为该像素着色。通过一个额外的蒙板矩阵也可以限制浸水填充的范围。OpenCV提供的浸水填充函数有两个形式，它们的原型如下。

// image：待处理的输入图像，单通道或者三通道，8位整型或者32位浮点型数据
// seed：浸水填充开始的种子点
// newVal：填充的颜色值
// rect：输出参数，描述重绘区域的最小矩形
// lowDiff：向下的强度偏差，即【当前像素强度 > 参考像素强度 - lowDiff】时
//         并且满足highDiff的限制，则当前像素会被填充
// highDiff：向上的强度偏差，即【当前像素强度 < 参考像素强度 + highDiff】时
//          并且满足lowDiff的限制，则当前像素会被填充
// flags：算法内部计算策略，包括像素比较策略等等，下文介绍
int cv::floodFill(cv::InputOutputArray image,
                  cv::Point seed, cv::Scalar newVal, cv::Rect* rect,
                  cv::Scalar lowDiff = cv::Scalar(),
                  cv::Scalar highDiff = cv::Scalar(),
                  int flags);

// image：待处理的输入图像，单通道或者三通道，8位整型或者32位浮点型数据，尺寸为W✖️H
// mask：蒙板矩阵，尺寸为W+2✖️H+2，单通道，位深度为8，用于控制填充范围等，下文介绍
int cv::floodFill(cv::InputOutputArray image, cv::InputOutputArray mask, 
                  cv::Point seed, cv::Scalar newVal, cv::Rect* rect,
                  cv::Scalar lowDiff = cv::Scalar(),
                  cv::Scalar highDiff = cv::Scalar(),
                  int flags);

通常情况下，该函数会修改输入图像矩阵image的数据。该函数首先将指定的种子点涂成参数newVal指定的颜色，然后向周围扩散处理其他的像素。如果在参数flags指定的比较策略下，像素的颜色强度满足参数lowDiff和highDiff指定的条件，则当前像素也会被涂成指定的颜色。参数newVal，loDiff和upDiff的数据类型都是cv::Scalar，也就是我们可以分别限制RGB三个颜色通道的像素强度偏差范围。

参数mask的宽高都应当比输入图像的尺寸大2，可以看作是向上下左右四个方向各扩展一个维度。另外该参数可以同时被看作是输入和输出参数，作为输入参数而言只有mask矩阵中非0值对应的像素才有可能被上色。作为输出参数而言，当算法运行后实际被上色的像素对应的元素将会被设置为指定的非0值。

需要注意的是蒙板矩阵的基本数据类型为8位整型，如果你直接在窗口中展示该图，则有可能看见的是一副黑色的图像。这可能是因为其中的元素值太小，因为对于格式为8位整型的图像而言，其像素强度的取值空间为0到255，你只需要适当对原数据进行缩放即可。

参数flags包含了多个含义，它是一个32位的整型数据，其中低8位(0-7)用于指定算法在考虑相邻像素的连通策略。可选值有4和8，分别表示4连通和8连通，前者表示浸水算法应当考虑上下左右四个相邻像素，后者表示在前者基础上还应当考虑对角线上的4个共8个相邻像素。

中间8位(8-15)用于指定满足前文描述的相关条件时，蒙板元素应当被填充的值，如果设置为0或者未设置，则默认使用1。

高8位(16-23)则用于指定像素比较的策略以及是否需要修改原图，当flags包含cv::FLOODFILL_FIXED_RANGE时，比较的参考像是是选定的种子像素，否则比较相邻像素，当flags包含cv::FLOODFILL_MASK_ONLY时，该函数只会更新蒙板矩阵mask的值，而不会修改原图。

参数flags在3个位段段选项可以用逻辑与符号连接得到最终的值，例如如果需要使用8连通的相邻像素推断策略，并且仅和种子像素比较(原书中同时设置该策略和连通性，此时连通性的设置可能会被忽略，推测此处仅为演示flags参数的设置方法)，并且仅更新蒙板矩阵，并且使用47填充蒙板矩阵，则参数flags的设置方式如下。

flags = 8
        | cv::FLOODFILL_MASK_ONLY
        | cv::FLOODFILL_FIXED_RANGE
        | (47 << 8);

下图展示了一个浸水算法应用于图像处理的结果，黑色小圆圈表示的是选定的种子像素，它们的像素强度偏差范围参数都设置为7，其中上图的填充色为灰色，下图为白色。

下图是对通一副图像应用不同参数的浸水算法对处理的结果，与上一个示例不同的是此次选择了cv::FLOODFILL_FIXED_RANGE的比较策略，通过与种子点比较从而获得更大范围的填充效果(推断比较相邻点需要所有相邻像素满足限制条件，而与种子点比较只考虑当前像素，因此能够沿着细小裂隙渗透到更大的范围)。

7.2 分水岭算法

很多时候想要分割一副图像，但是却没有可用的图片分割背景蒙板，此时可以使用分水岭算法(Watershed Algorithm)来解决问题。该算法通过计算图像的像素强度梯度来确定边界，并将这些边界转化为山脊，被这些山脊分割的区域被看作是谷地或者盆地。用户需要预先指定一些点作为不同区块的标识点，算法随后会从这些标识点向这些盆地灌水直至填充满整个图像，并使用标识点的值来标记其被包含的区块，从而达到分割图像的目的。

此外你也可以在原始图像中绘制一系列的线条，每个线条都表示你想要分割出的一块区域，并在标记矩阵中使用不同的值表示每个线条。在已经通过梯度计算处理好的山脊-盆地图中，分水岭算法将会把这些线条跨越的以及临近的盆地合并，并赋予其与对应线条相同的标记，从而达到图像分割的目的。

下图是一个通过线条标识应用分水岭算法处理图像的示例，左图中每个线条都标识了一个独立的分区，在右图中分水岭算法成功将图像分割为多个区块。

分水岭算法的函数原型如下。

// image：待处理的图像，通道数为3，数据类型为8位整型
// markers：标记矩阵，通道数为1，数据类型为32位整型，尺寸和image相同
void cv::watershed(cv::InputArray image, cv::InputOutputArray markers);

该函数传入的参数markers矩阵中，除了认为属于同一个区域的像素点使用同一个正值表示，如上图中每条直线对应的点都使用同一个正值表示，其余元素都应设置为0。当算法执行完毕后，对于markers矩阵中的所有非零像素，如果为边界像素将会被设置为-1，否则将会被设置为其所属区域的标记值。理想情况下所有区域都应该被值为-1表示的边界分割，但是实际上并不能得到这样的结果。很简单的一个例子就是如果传入的标记矩阵中两个相邻像素的值不为0，并且不相同，则输出结果中它们会保持不变，并且也不会被-1标识的边界所分隔。

7.3 Crabcuts算法

Rother、Kolmogorov和Blake对Graphcuts算法进行改进得到Grabcuts算法，该算法能够用于处理用户主导的前景和背景图像分割问题。该算法能够得到很好的分割结果，通常不会超过待分割的前景图的限制矩形(调用者指定的)区域。

传统的Graphcuts算法使用用户标记的背景区域和前景区域建立分布直方图。然后假设在理想情况下未标记的背景和前景区域应当和已经标记的对应区域有相似的分布，当然理想情况下这些区域应当是平滑并且相连的，实际上却可能包含很多斑点。这些假设最终会组合成能量函数(Energy Functional)，该能量函数会给出一个遵守这些假设的低成本解决方案，以及一个不遵守这些假设的高成本解决方案。算法通过最小化能量函数计算最终的结果，通常使用一种称为Mincut的技术，这也是Graphcuts算法和Grabcuts命名的由来。

Grabcuts通过如下几个重要的方式改进了Graphcuts算法。首先它使用了高斯混合模型(Gaussian Mixture Model)替代了直方图模型，因此能够处理彩色图像。另外它使用迭代的方式来解决能量函数最小化问题，因此提供更好的全局结果，并且增加了用户提供标记的灵活性。此外这种方式使算法能够支持单一标记，即用户只需要前景图或者背景图，而使用Graphcuts算法需要同时标记这两个区域。

OpenCV提供的函数原型如下。

// img：待分割的图像
// mask：前景背景标识蒙板，通道数为1，数据类型为cv::U8
// rect：待分割的区域矩形，矩形外部都被认为是背景区域
// bgdModel：背景缓存，下文介绍
// fgdModel：前景缓存，下文介绍
// iterCount：内部Graphcuts算法的迭代次数
// mode：算法策略，下文介绍
void cv::grabCut(cv::InputArray img, cv::InputOutputArray mask, cv::Rect rect,
                 cv::InputOutputArray bgdModel, cv::InputOutputArray fgdModel,
                 int iterCount, int mode = cv::GC_EVAL);

参数mask可以同时看作是输入矩阵和输出矩阵，当参数mode包含cv::GC_INIT_WITH_MASK时，该值作为输入矩阵参与算法的计算，矩阵中的每个元素取值必须符合该函数的规定，它们可取值及其含义如下表。

蒙板矩阵元素取值常量	字面值	描述
cv::GC_BGD	0	确定的背景区域
cv::GC_FGD	1	确定的前景区域
cv::PR_GC_BGD	2	疑似的背景区域
cv::PR_GC_FGD	3	疑似的前景区域

此时算法会通过分析标记的确定区域识别疑似区域为确定的区域。作为输出矩阵而言，该矩阵用于保存函数识别的结果。当不使用mask矩阵初始化原图中的前景和背景区域时，需要提供参数rect，并且此时参数mode应当包含cv::GC_INIT_WITH_RECT。此时矩形的外部区域会被认为是确定的背景区域，而矩形内部区域会被认为是疑似的前景区域。

参数bgdModel和fgdModel分别是前景和背景缓存，当第一次调用该函数执行iterCount次迭代时可以直接传入空矩阵。当因为某个原因需要再次调用该函数再执行几次迭代时，可能是用户提供了额外的确定前景或背景信息，除了将上一次计算的到的蒙板矩阵作为本次函数调用的输入外，你还需要将上一次该函数调用后的返回的缓存矩阵作为参数传入到当前次的调用函数中。

Grabcuts算法内部会多次调用Graphcuts算法，当然包含前文提到的一些改进，并且每次调用都会重新计算混合模型。参数iterCount指定了迭代的次数，通常设置为10或者12，不过它的大小也还需取决于被处理图像的大小和性质。

7.4 Mean-Shift分割算法

Mean-Shift图像分割算法用于寻找颜色再空间分布上的峰值，它和Mean-Shift算法相关，后者将会在讨论运动和追踪时详细聊到。它们之间的区别是前者关注的是颜色在空间上的分别，而后者通过连续的帧在时间上追踪这些分布。

给定维度为(x, y, blue, green, red)的数据点集，mean-shift算法可以通过窗口扫描指定在空间找到密度最高的块。空间变量(x, y)的范围可以和颜色强度范围(blue, green, red)相差较大，因此该算法需要两个不同的窗口半径分别计算空间和颜色变量。随着mean-shift滑动窗口移动，窗口遍历到的所有收敛于峰值的点都应当连接或拥有该峰值。这种所有权从最密集的峰值处向外辐射，从而达到了分割图像的目的。

图像分割实际上是在图像金字塔中完成的，即使用道理前面章节中介绍的函数cv::pyrUp()和cv::pyrDown()。在图像金字塔的高层图像中的颜色聚类(Color Clusters)会在低层的图像中有更明确的边界。

Mean-Shift图像分割算法的输出是一副色调分离(Posterized)的新图像，即其中的纹理细节被移除，颜色梯度也变得更加平坦。你可以使用合适的算法进一步处理该算法的输出结果，如使用函数cv::Canny()和cv::findCountours()寻找图像中的轮廓。OpenCV提供的函数原型如下。

// src：待处理的图像，通道数为3，数据类型为8位整型
// dst：处理完成的图像，通道数为3，数据类型为8位整型，尺寸和src相同
// sp：滑动窗口的空间半径
// sr：滑动窗口的颜色半径
// maxLevel：图像金字塔的最大层数
// termcrit：迭代终止条件，不设置时使用默认值
void cv::pyrMeanShiftFiltering(cv::InputArray src, cv::OutputArray dst,
                               cv::double sp, cv::double sr, int maxLevel = 1,
                               cv::TermCriteria termcrit = TermCriteria(
                                cv::TermCriteria::MAX_ITER | cv::TermCriteria::EPS,
                                5, 1));

对于一个640✖️480点彩色图像而言，参数sp和sr的建议值分别为2和40，参数maxLevel建议值为2或者3。该算法内部会包含迭代逻辑，而参数termcrit可以指定迭代逻辑的终止条件，实际使用是可以将该参数留空使用默认值。

下图是一个Mean-Shift图像分割示例的运行结果，其中参数sp和sr分别设置为2和40，参数maxLevel设置为2。

8 小结

前一章节介绍了一般的图像变换技术，而在本章中详细讨论了图像分析技术，这些技术使我们对图像有了更深入的了解，睡着本系列文章的不断深入，我们将会学到更多更复杂的图像处理技术，而本章及上一章学到的这些知识将是这些高级图像处理技术的基础。