图像特征提取（五）—— Harris角点检测（一）

客官°小女子只卖身不卖艺 2022-05-20 00:25 335阅读 0赞

# 转载自：http://www.cnblogs.com/ronny/p/4009425.html    #

# 1. 不同类型的角点 #

在现实世界中，角点对应于物体的拐角，道路的十字路口、丁字路口等。从图像分析的角度来定义角点可以有以下两种定义：

1.  角点可以是两个边缘的角点；
2.  角点是邻域内具有两个主方向的特征点；

前者往往需要对图像边缘进行编码，这在很大程度上依赖于图像的分割与边缘提取，具有相当大的难度和计算量，且一旦待检测目标局部发生变化，很可能导致操作的失败。早期主要有Rosenfeld和Freeman等人的方法，后期有CSS等方法。

基于图像灰度的方法通过计算点的曲率及梯度来检测角点，避免了第一类方法存在的缺陷，此类方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等。

![image][]

这篇文章主要介绍的Harris角点检测的算法原理，比较著名的角点检测方法还有jianbo Shi和Carlo Tomasi提出的Shi-Tomasi算法，这个算法开始主要是为了解决跟踪问题，用来衡量两幅图像的相似度，我们也可以把它看为Harris算法的改进。OpenCV中已经对它进行了实现，接口函数名为[GoodFeaturesToTrack()][GoodFeaturesToTrack]。另外还有一个著名的角点检测算子即SUSAN算子，SUSAN是Smallest Univalue Segment Assimilating Nucleus（最小核值相似区）的缩写。SUSAN使用一个圆形模板和一个圆的中心点，通过圆中心点像素与模板圆内其他像素值的比较，统计出与圆中心像素近似的像元数量，当这样的像元数量小于某一个阈值时，就被认为是要检测的角点。我觉得可以把SUSAN算子看为Harris算法的一个简化。这个算法原理非常简单，算法效率也高，所以在OpenCV中，它的接口函数名称为：[FAST()][FAST] 。

# 2. Harris角点 #

## 2.1 基本原理 ##

人眼对角点的识别通常是在一个局部的小区域或小窗口完成的。如果在各个方向上移动这个特征的小窗口，窗口内区域的灰度发生了较大的变化，那么就认为在窗口内遇到了角点。如果这个特定的窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点；如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条直线的线段。

![image][image 1]

对于图像I(x,y)I(x,y)，当在点(x,y)(x,y)处平移(Δx,Δy)(Δx,Δy)后的自相似性，可以通过自相关函数给出：

c(x,y;Δx,Δy)=∑(u,v)∈W(x,y)w(u,v)(I(u,v)–I(u\+Δx,v\+Δy))2c(x,y;Δx,Δy)=∑(u,v)∈W(x,y)w(u,v)(I(u,v)–I(u+Δx,v+Δy))2

其中，W(x,y)W(x,y)是以点(x,y)(x,y)为中心的窗口，w(u,v)w(u,v)为加权函数，它既可是常数，也可以是高斯加权函数。

![image][image 2]

根据泰勒展开，对图像I(x,y)I(x,y)在平移(Δx,Δy)(Δx,Δy)后进行一阶近似：

I(u\+Δx,v\+Δy)=I(u,v)\+Ix(u,v)Δx\+Iy(u,v)Δy\+O(Δx2,Δy2)≈I(u,v)\+Ix(u,v)Δx\+Iy(u,v)ΔyI(u+Δx,v+Δy)=I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy+O(Δx2,Δy2)≈I(u,v)+Ix(u,v)Δx+Iy(u,v)Δy

其中，Ix,IyIx,Iy是图像I(x,y)I(x,y)的偏导数，这样的话，自相关函数则可以简化为：

c(x,y;Δx,Δy)≈∑w(Ix(u,v)Δx\+Iy(u,v)Δy)2=\[Δx,Δy\]M(x,y)\[ΔxΔy\]c(x,y;Δx,Δy)≈∑w(Ix(u,v)Δx+Iy(u,v)Δy)2=\[Δx,Δy\]M(x,y)\[ΔxΔy\]

其中

M(x,y)=∑w\[Ix(x,y)2Ix(x,y)Iy(x,y)Ix(x,y)Iy(x,y)Iy(x,y)2\]=\[∑wIx(x,y)2∑wIx(x,y)Iy(x,y)∑wIx(x,y)Iy(x,y)∑wIy(x,y)2\]=\[ACCB\]M(x,y)=∑w\[Ix(x,y)2Ix(x,y)Iy(x,y)Ix(x,y)Iy(x,y)Iy(x,y)2\]=\[∑wIx(x,y)2∑wIx(x,y)Iy(x,y)∑wIx(x,y)Iy(x,y)∑wIy(x,y)2\]=\[ACCB\]

也就是说图像I(x,y)I(x,y)在点(x,y)(x,y)处平移(Δx,Δy)(Δx,Δy)后的自相关函数可以近似为二项函数：

c(x,y;Δx,Δy)≈AΔx2\+2CΔxΔy\+BΔy2c(x,y;Δx,Δy)≈AΔx2+2CΔxΔy+BΔy2

其中

A=∑wI2x,B=∑wI2y,C=∑wIxIyA=∑wIx2,B=∑wIy2,C=∑wIxIy

二次项函数本质上就是一个椭圆函数。椭圆的扁率和尺寸是由M(x,y)M(x,y)的特征值λ1、λ2λ1、λ2决定的，椭贺的方向是由M(x,y)M(x,y)的特征矢量决定的，如下图所示，椭圆方程为：

\[Δx,Δy\]M(x,y)\[ΔxΔy\]=1\[Δx,Δy\]M(x,y)\[ΔxΔy\]=1

![image][image 3]

椭圆函数特征值与图像中的角点、直线（边缘）和平面之间的关系如下图所示。共可分为三种情况：

*  图像中的直线。一个特征值大，另一个特征值小，λ1≫λ2或λ2≫λ1。自相关函数值在某一方向上大，在其他方向上小。
 *  图像中的平面。两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。
 *  图像中的角点。两个特征值都大，且近似相等，自相关函数在所有方向都增大。
 *  ![70][]  
    
 *  根据二次项函数特征值的计算公式，我们可以求M(x,y)M(x,y)矩阵的特征值。但是Harris给出的角点差别方法并不需要计算具体的特征值，而是计算一个角点响应值RR来判断角点。RR的计算公式为：
    
      
    R=detM−α(traceM)2
    
    式中，detMdetM为矩阵M=\[ABBC\]M=\[ABBC\]的行列式；traceMtraceM为矩阵MM的直迹；αα为经常常数，取值范围为0.04~0.06。事实上，特征是隐含在detMdetM和traceMtraceM中，因为，
    
    detM=λ1λ2=AC−B2
    
    traceM=λ2\+λ2=A\+C
    
    ## 2.2 Harris角点算法实现 ##
    
    根据上述讨论，可以将Harris图像角点检测算法归纳如下，共分以下五步：
    
    1. 计算图像I(x,y)I(x,y)在XX和YY两个方向的梯度Ix、IyIx、Iy。
    
    Ix=∂I∂x=I⊗(−1 0 1)，Iy=∂I∂x=I⊗(−1 0 1)T
    
    2. 计算图像两个方向梯度的乘积。
    
      
    Ix2=Ix⋅Iy，Iy2=Iy⋅Iy，Ixy=Ix⋅Iy
    
    3. 使用高斯函数对I2x、I2y和IxyIx2、Iy2和Ixy进行高斯加权（取σ=1σ=1），生成矩阵MM的元素A、BA、B和C。
    
      
    A=g(Ix2)=Ix2⊗w，C=g(Iy2)=Iy2⊗w，B=g(Ix,y)=Ixy⊗w
    
    4. 计算每个像素的Harris响应值RR，并对小于某一阈值tt的RR置为零。
    
    R=\{R:detM−α(traceM)2<t\}
    
    5. 在3×33×3或5×55×5的邻域内进行非最大值抑制，局部最大值点即为图像中的角点。
    
    > Harris角点检测的C++实现代码：[https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harris.cpp][https_github.com_RonnyYoung_ImageFeatures_blob_master_source_harris.cpp]
    
    ## 2.3 Harris角点的性质 ##
    
    ### 1. 参数αα对角点检测的影响 ###
    
    假设已经得到了矩阵MM的特征值λ1≥λ2≥0λ1≥λ2≥0，令λ2=kλ1,0≤k≤1λ2=kλ1,0≤k≤1。由特征值与矩阵MM的直迹和行列式的关系可得：
    
      
    detM=∏iλi      traceM=∑iλi
    
    从而可以得到角点的响应
    
    R=λ2λ2=α(λ2+λ2)2=λ2(k−α(1+k)2)
    
    假设R≥0R≥0，则有：
    
      
    0≤α≤k(1+k)2≤0.25
    
    对于较小的kk值，R≈λ2(k−α),α<kR≈λ2(k−α),α<k。
    
    由此，可以得出这样的结论：增大αα的值，将减小角点响应值RR，降低角点检测的灵性，减少被检测角点的数量；减小αα值，将增大角点响应值RR，增加角点检测的灵敏性，增加被检测角点的数量。
    
    ### 2. Harris角点检测算子对亮度和对比度的变化不敏感 ###
    
    这是因为在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之，对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。
    
    ![70 1][]  
    
    
                        ![70 2][]
 *  ### 3. Harris角点检测算子具有旋转不变性 ###
    
    Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值RR也不发生变化，由此说明Harris角点检测算子具有旋转不变性。
    
    ### 4. Harris角点检测算子不具有尺度不变性 ###
    
    如下图所示，当右图被缩小时，在检测窗口尺寸不变的前提下，在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。
    
    ![70 3][]

## 2.4 Harris的OpenCV接口 ##

OpenCV的Hairrs角点检测的函数为cornerHairrs()，但是它的输出是一幅浮点值图像，浮点值越高，表明越可能是特征角点，我们需要对图像进行阈值化。

C++: void cornerHarris(InputArray src, OutputArray dst, int blockSize, int apertureSize, double k, int borderType = BORDER_DEFAULT);

*  **src** – 输入的单通道8-bit或浮点图像。
 *  **dst** – 存储着Harris角点响应的图像矩阵，大小与输入图像大小相同，是一个浮点型矩阵。
 *  **blockSize** – 邻域大小。
 *  **apertureSize** – 扩展的微分算子大。
 *  **k** – 响应公式中的，参数αα。
 *  **boderType** – 边界处理的类型。
 *

int main()
    {
        Mat image = imread("../buliding.png");
        Mat gray;
        cvtColor(image, gray, CV_BGR2GRAY);
    
        Mat cornerStrength;
        cornerHarris(gray, cornerStrength, 3, 3, 0.01);
        threshold(cornerStrength, cornerStrength, 0.001, 255, THRESH_BINARY);
        return 0;
    }

![70 4][]![70 5][]![70 6][]

从上面上间一幅图像我们可以看到，有很多角点都是粘连在一起的，我们下面通过加入非极大值抑制来进一步去除一些粘在一起的角点。

非极大值抑制原理是，在一个窗口内，如果有多个角点则用值最大的那个角点，其他的角点都删除，窗口大小这里我们用3\*3，程序中通过图像的膨胀运算来达到检测极大值的目的，因为默认参数的膨胀运算就是用窗口内的最大值替代当前的灰度值。

int main()
    {
        Mat image = imread("buliding.png");
        Mat gray;
        cvtColor(image, gray, CV_BGR2GRAY);
    
        Mat cornerStrength;
        cornerHarris(gray, cornerStrength, 3, 3, 0.01);
    
        double maxStrength;
        double minStrength;
        // 找到图像中的最大、最小值
        minMaxLoc(cornerStrength, &minStrength, &maxStrength);
    
        Mat dilated;
        Mat locaMax;
        // 膨胀图像，最找出图像中全部的局部最大值点
        dilate(cornerStrength, dilated, Mat());
        // compare是一个逻辑比较函数，返回两幅图像中对应点相同的二值图像
        compare(cornerStrength, dilated, locaMax, CMP_EQ);
    
        Mat cornerMap;
        double qualityLevel = 0.01;
        double th = qualityLevel*maxStrength; // 阈值计算
        threshold(cornerStrength, cornerMap, th, 255, THRESH_BINARY);
        cornerMap.convertTo(cornerMap, CV_8U);
        // 逐点的位运算
        bitwise_and(cornerMap, locaMax, cornerMap);
    
        drawCornerOnImage(image, cornerMap);
        namedWindow("result");
        imshow("result", image);
        waitKey();
    
        return 0;
    }
    void drawCornerOnImage(Mat& image, const Mat&binary)
    {
        Mat_<uchar>::const_iterator it = binary.begin<uchar>();
        Mat_<uchar>::const_iterator itd = binary.end<uchar>();
        for (int i = 0; it != itd; it++, i++)
        {
            if (*it)
                circle(image, Point(i%image.cols, i / image.cols), 3, Scalar(0, 255, 0), 1);
        }
    }

现在我们得到的效果就比默认的函数得到的结果有相当的改善，如上面最右边效果图。

# 3. 多尺度Harris角点 #

## 3.1 多尺度Harris角点的原理 ##

虽然Harris角点检测算子具有部分图像灰度变化的不变性和旋转不变性，但它不具有尺度不变性。但是尺度不变性对图像特征来说至关重要。人们在使用肉眼识别物体时，不管物体远近，尺寸的变化都能认识物体，这是因为人的眼睛在辨识物体时具有较强的尺度不变性。在[图像特征提取：尺度空间理论][Link 1]这篇文章里就已经讲到了高斯尺度空间的概念。下面将Harris角点检测算子与高斯尺度空间表示相结合，使用Harris角点检测算子具有尺度的不变性。

仿照Harris角点检测中二阶矩的表示方法，使用M=μ(x,σI,σD)M=μ(x,σI,σD)为尺度自适应的二阶矩：

M=μ(x,σI,σD)=σD2g(σI)⊗\[Lx2(x,σD)LxLy(x,σD)LxLy(x,σD)Ly2(x,σD)\]

其中，g(σI)g(σI)表示尺度为sigmaIsigmaI的高斯卷积核，xx表示图像的位置。与高斯测度空间类似，使用L(x)L(x)表示经过高斯平滑后的图像，符号⊗⊗表示卷积，Lx(x,σD)Lx(x,σD)和Ly(x,σD)Ly(x,σD)表示对图像使用高斯g(σD)g(σD)函数进行平滑后，在xx或yy方向取其微分的结果，即Lx=∂xLLx=∂xL和Ly=∂yLLy=∂yL。通常将σIσI称为积分尺度，它是决定Harris角点当前尺度的变量，σDσD为微分尺度或局部尺度，它是决定角点附近微分值变化的变量。显然，积分尺度σIσI应该大于微分尺度σDσD。

## 3.2 多尺度Harris角点实现 ##

首先，检测算法从预先定义的一组尺度中进行积分尺度搜索，这一组尺度定义为

σ1…σn=σ0…knσ0

一般情况下使用k=1.4。为了减少搜索的复杂性，对于微分尺度σDσD的选择，我们采用在积分尺度的基础上，乘以一个比例常数，即σD=sσIσD=sσI，一般取s=0.7。这样，通常使用积分和微分的尺度，便可以生成μ(x,σI,σD)μ(x,σI,σD)，再利用Harris角点判断准则，对角点进行搜索，具体可以分两步进行。

1. 与Harris角点搜索类似，对于给定的尺度空间值σDσD，进行如下角点响应值计算和判断：

cornerness=det(μ(x,σn)−αtrace2(μ(x,σn)))>thresholdH

2. 对于满足1中条件的点，在点的8邻域内进行角点响应最大值搜索（即非最大值抑制）出在8邻域内角点响应最大值的点。对于每个尺度σn(1,2,…,n)σn(1,2,…,n)都进行如上搜索。

由于位置空间的候选点并不一定在尺度空间上也能成为候选点，所以，我们还要在尺度空间上进行搜索，找到该点的所谓特征尺度值。搜索特征尺度值也分两步。

1. 对于位置空间搜索到的每个候选点，进行拉普拉斯响应计算，并满足其绝对值大于给定的阈值条件：

F(x,σn)=σn2|Lxx(x,σn)+Lyy(x,σn)|≥thresholdL

2. 与邻近的两个尺度空间的拉普拉斯响应值进行比较，使其满足：

F(x,σn)>F(x,σl),   l∈\{n−1.n+1\}

满足上述条件的尺度值就是该点的特征尺度值。这样，我们就找到了在位置空间和尺度空间都满足条件的Harris角点。

> 多尺度Harris角点检测C++实现：[https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harrisLaplace.cpp][https_github.com_RonnyYoung_ImageFeatures_blob_master_source_harrisLaplace.cpp]

# 4. 更多的讨论 #

在上面描述的Harris角点具有光照不变性、旋转不变性、尺度不变性，但是严格意义上来说并不具备仿射不变性。Harris-Affine是一种新颖的检测仿射不变特征点的方法，可以处理明显的仿射变换，包括大尺度变化和明显的视角变化。Harris-Affine主要是依据了以下三个思路：

1.  特征点周围的二阶矩的计算对区域进行的归一化，具有仿射不变性；
2.  通过在尺度空间上归一化微分的局部极大值求解来精化对应尺度；
3.  自适应仿射Harris检测器能够精确定位牲点；

这篇文章不对Harris-Affine作进一步的描述，有时间会对这一算法做专门的分析，有兴趣的可以参考原论文：[Scale & Affine Invariant Interest Point Detectors][Scale _ Affine Invariant Interest Point Detectors].

# 5. 参考资料 #

\[1\] 《图像局部不变特征与描述》王永明，王贵锦。

\[2\] [Harris角点及Shi-Tomasi角点检测][Harris_Shi-Tomasi]

\[3\] [图像特征提取PPT][PPT]

\[4\] [Harris角点检测算法 1][Harris_ 1]

\[5\] [OpenCV Harris角点检测][OpenCV Harris]

\[6\] [Opencv学习笔记（五）Harris角点检测][Opencv_Harris]

[image]: /images/20220520/1d03534a95a74ebda6e0d9266ee7a7ac.png
[GoodFeaturesToTrack]: http://www.opencv.org.cn/opencvdoc/2.3.2/html/modules/imgproc/doc/feature_detection.html#goodfeaturestotrack
[FAST]: http://docs.opencv.org/modules/features2d/doc/feature_detection_and_description.html#fast
[image 1]: /images/20220520/dfd58dab14fe4c11ae118c033de2a5cb.png
[image 2]: /images/20220520/906006f8a46545788bc5e641c3a8e35b.png
[image 3]: /images/20220520/64392023215440aeb9ba2aff39a25fec.png
[70]: /images/20220520/793910294f8a45f9b5862b412bb3815f.png
[https_github.com_RonnyYoung_ImageFeatures_blob_master_source_harris.cpp]: https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harris.cpp
[70 1]: /images/20220520/a65c969297e74348abb0f0a6dc1d8e9e.png
[70 2]: /images/20220520/f1f26c7a1d9340fc8fec6e08aa8054a0.png
[70 3]: /images/20220520/49165e477d4142aebd2ad9f275ce6235.png
[70 4]: https://img-blog.csdn.net/20180709154237753?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM3NzY0MTI5/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
[70 5]: https://img-blog.csdn.net/20180709154224520?watermark/2/text/aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L3FxXzM3NzY0MTI5/font/5a6L5L2T/fontsize/400/fill/I0JBQkFCMA==/dissolve/70
[70 6]: /images/20220520/abd24a77a7604db699de3e903836d862.png
[Link 1]: http://www.cnblogs.com/ronny/p/3886013.html
[https_github.com_RonnyYoung_ImageFeatures_blob_master_source_harrisLaplace.cpp]: https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harrisLaplace.cpp
[Scale _ Affine Invariant Interest Point Detectors]: https://github.com/RonnyYoung/ImageFeatures/blob/master/papers/corners/Scale%26Affine_Invariant_Interest_Point_Detectors..pdf
[Harris_Shi-Tomasi]: http://blog.csdn.net/xw20084898/article/details/21180729
[PPT]: http://wenku.baidu.com/link?url=OoAlTQjV3_WSi9ycEopQNL0FyWqHG2Sh_whAhnwqFB-R6Hi7Z-Sjsv_eM-ei4iugs9qu8LUF6uw1dirWTQvFTmer06mr9-2yvWPAWmM2-Fm
[Harris_ 1]: http://www.cnblogs.com/ztfei/archive/2012/05/07/2487123.html
[OpenCV Harris]: http://docs.opencv.org/modules/imgproc/doc/feature_detection.html?highlight=harris#cv.CornerHarris
[Opencv_Harris]: http://blog.csdn.net/crzy_sparrow/article/details/7391511