Harris角点

最新推荐文章于 2023-01-06 21:30:08 发布

自由人_SZ

最新推荐文章于 2023-01-06 21:30:08 发布

阅读量538

点赞数

Harris角点

1. 不同类型的角点

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

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

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

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

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

2. Harris角点

2.1 基本原理

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

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

c (x, y; Δ x, Δ y) = \sum (u, v) \in W (x, y) w (u, v) (I (u, v) - I (u + Δ x, v + Δ y)) 2

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

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

I (u + Δ x, v + Δ y) = I (u, v) + I x (u, v) Δ x + I y (u, v) Δ y + O (Δ x 2, Δ y 2) \approx I (u, v) + I x (u, v) Δ x + I y (u, v) Δ y

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

c (x, y; Δ x, Δ y) \approx \sum w (I x (u, v) Δ x + I y (u, v) Δ y) 2 = [Δ x, Δ y] M (x, y) [Δ x Δ y]

其中

M (x, y) = \sum w [I x (x, y) 2 I x (x, y) I y (x, y) I x (x, y) I y (x, y) I y (x, y) 2] = [\sum w I x (x, y) 2 \sum w I x (x, y) I y (x, y) \sum w I x (x, y) I y (x, y) \sum w I y (x, y) 2] = [A C C B]

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

c (x, y; Δ x, Δ y) \approx A Δ x 2 + 2 C Δ x Δ y + B Δ y 2

其中

A = \sum w I 2 x, B = \sum w I 2 y, C = \sum w I x I y

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

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

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

图像中的直线。一个特征值大，另一个特征值小， λ1≫λ2 或 λ2≫λ1 。自相关函数值在某一方向上大，在其他方向上小。
图像中的平面。两个特征值都小，且近似相等；自相关函数数值在各个方向上都小。
图像中的角点。两个特征值都大，且近似相等，自相关函数在所有方向都增大。

根据二次项函数特征值的计算公式，我们可以求 M(x,y) 矩阵的特征值。但是Harris给出的角点差别方法并不需要计算具体的特征值，而是计算一个角点响应值 R 来判断角点。 R 的计算公式为：

R = d e t M - α (t r a c e M) 2

式中， detM 为矩阵 M=[ABBC] 的行列式； traceM 为矩阵 M 的直迹； α 为经常常数，取值范围为0.04~0.06。事实上，特征是隐含在 detM 和 traceM 中，因为，

d e t M = λ 1 λ 2 = A C - B 2

t r a c e M = λ 2 + λ 2 = A + C

2.2 Harris角点算法实现

根据上述讨论，可以将Harris图像角点检测算法归纳如下，共分以下五步：

1. 计算图像 I(x,y) 在 X 和 Y 两个方向的梯度 Ix、Iy 。

I x = \partial I \partial x = I \otimes (- 101) ， I y = \partial I \partial x = I \otimes (- 101) T

2. 计算图像两个方向梯度的乘积。

I 2 x = I x \cdot I y ， I 2 y = I y \cdot I y ， I x y = I x \cdot I y

3. 使用高斯函数对 I2x、I2y和Ixy 进行高斯加权（取 σ=1 ），生成矩阵 M 的元素 A、B 和 C 。

A = g (I 2 x) = I 2 x \otimes w ， C = g (I 2 y) = I 2 y \otimes w ， B = g (I x, y) = I x y \otimes w

4. 计算每个像素的Harris响应值 R ，并对小于某一阈值 t 的 R 置为零。

R = {R : d e t M - α (t r a c e M) 2 < t}

5. 在 3×3 或 5×5 的邻域内进行非最大值抑制，局部最大值点即为图像中的角点。

Harris角点检测的C++实现代码：https://github.com/RonnyYoung/ImageFeatures/blob/master/source/harris.cpp

2.3 Harris角点的性质

1. 参数 α 对角点检测的影响

假设已经得到了矩阵 M 的特征值 λ1≥λ2≥0 ，令 λ2=kλ1,0≤k≤1 。由特征值与矩阵 M 的直迹和行列式的关系可得：

d e t M = \prod i λ i t r a c e M = \sum i λ i

从而可以得到角点的响应

R = λ 2 λ 2 = α (λ 2 + λ 2) 2 = λ 2 (k - α (1 + k) 2)

假设 R≥0 ，则有：

0 \leq α \leq k ( 1 + k ) 2 \leq 0.25

对于较小的 k 值， R≈λ2(k−α),α<k 。

由此，可以得出这样的结论：增大 α 的值，将减小角点响应值 R ，降低角点检测的灵性，减少被检测角点的数量；减小 α 值，将增大角点响应值 R ，增加角点检测的灵敏性，增加被检测角点的数量。

2. Harris角点检测算子对亮度和对比度的变化不敏感

这是因为在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之，对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。

3. Harris角点检测算子具有旋转不变性

Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值 R 也不发生变化，由此说明Harris角点检测算子具有旋转不变性。

4. Harris角点检测算子不具有尺度不变性

如下图所示，当右图被缩小时，在检测窗口尺寸不变的前提下，在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。

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;
}

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

非极大值抑制原理是，在一个窗口内，如果有多个角点则用值最大的那个角点，其他的角点都删除，窗口大小这里我们用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角点检测算子具有部分图像灰度变化的不变性和旋转不变性，但它不具有尺度不变性。但是尺度不变性对图像特征来说至关重要。人们在使用肉眼识别物体时，不管物体远近，尺寸的变化都能认识物体，这是因为人的眼睛在辨识物体时具有较强的尺度不变性。在图像特征提取：尺度空间理论这篇文章里就已经讲到了高斯尺度空间的概念。下面将Harris角点检测算子与高斯尺度空间表示相结合，使用Harris角点检测算子具有尺度的不变性。

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

M = μ (x, σ I, σ D) = σ 2 D g (σ I) \otimes [L 2 x (x, σ D) L x L y (x, σ D) L x L y (x, σ D) L 2 y (x, σ D)]

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