1. 概述

因为角点在现实生活场景中非常常见，因此，角点检测算法也是一种非常受欢迎的检测算法，尤其本文要讲的Harris角点检测，可以说传统检测算法中的经典之作
目前角点检测算法主要可归纳为3类
- 基于灰度图像的角点检测
- 基于二值图像的角点检测
- 基于轮廓的角点检测

2. 角点

角点就是轮廓之间的交点
像素点附近区域像素无论是在梯度方向、还是在梯度幅值上都发生较大的变化
角点可以是两个边缘的角点
- 需要对图像边缘进行编码，这在很大程度上依赖于图像的分割与边缘提取，具有相当大的难度和计算量，且一旦待检测目标局部发生变化，很可能导致操作的失败。早期主要有Rosenfeld和Freeman等人的方法，后期有CSS等方法
角点是邻域内具有两个主方向的特征点
- 基于图像灰度的方法通过计算点的曲率及梯度来检测角点，避免了第一类方法存在的缺陷，此类方法主要有Moravec算子、Forstner算子、Harris算子、SUSAN算子等

3. Harris角点

3.1 算法思想

选取一个固定的窗口在图像上以任意方向的滑动，如果灰度都有较大的变化，那么认为这个窗口内部存在角点
如果这个特定的窗口在图像各个方向上移动时，窗口内图像的灰度没有发生变化，那么窗口内就不存在角点
如果窗口在某一个方向移动时，窗口内图像的灰度发生了较大的变化，而在另一些方向上没有发生变化，那么，窗口内的图像可能就是一条直线的线段

3.2 数学推导

对于图像$I(x, y)$，当在点$(x, y)$处平移$(\Delta x, \Delta y)$后的自相似性由自相关函数计算：

$c(x, y; \Delta x, \Delta y) = \sum_{(u, v)\in W(x, y)}w(u, v)(I(u, v) - I(u + \Delta x, v + \Delta y))^2$,

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

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

$I(u + \Delta x, v + \Delta y) = I(u, v) + I_x(u, v)\Delta x + I_y(u, v)\Delta y + O(\Delta x^2, \Delta y^2) \approx I(u, v) + I_x(u, v)\Delta x + I_y(u, v)\Delta y$,

其中，$I_x$和$I_y$是图像$I(x, y)$的偏导数。

进一步地，自相关函数简化为

$c(x, y; \Delta x, \Delta y) \approx \sum_{w}(I_x(u, v)\Delta x + I_y(u, v)\Delta y)^2 = [\Delta x, \Delta y]M(x, y) \left[\begin{matrix} \Delta x \\ \Delta y \end{matrix} \right]$,

其中，

$M(x, y) = \sum_{w}\left[\begin{matrix} 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 \end{matrix} \right] = \left[\begin{matrix} A & C \\ C & B \end{matrix} \right]$.

根据主成分分析(PCA)的原理可知，如果对矩阵$M(x, y)$对角化，那么，特征值就是主分量上的方差，矩阵是二维的方阵，有两个主分量，如果在窗口区域内是角点，那么梯度变化会较大，像素点的梯度分布比较离散，这体现在特征值上就是特征值比较大
- 如果矩阵对应的两个特征值都较大，那么窗口内含有角点
- 如果特征值一个大一个小，那么窗口内含有线性边缘
- 如果两个特征值都很小，那么窗口内为平坦区域
角点的检测转化为数学模型，就是求解窗口内矩阵的特征值并且判断特征值的大小
Harris给出的角点差别方法并不需要计算具体的特征值，而是计算一个角点响应值$R$来判断角点

$R = det M - \alpha(trace M)^2$,

其中，$det M$为矩阵$M$的行列式，$trance M$为矩阵$M$的迹，$\alpha$为常数，取值范围为$0.04\sim0.06$。回头来看，特征值包含在$R$的计算公式中，即

$det M = \lambda_1\lambda_2 = AC - B^2$,

$trace M = \lambda_1 + \lambda_2 = A + C$.

3.3 算法实现

因为Harris角点检测算法非常经典，因此，一些成熟的图像处理或视觉库都会直接提供Harris角点检测的算法
以OpenCV为例，代码如下

import cv2
import numpy as np

"""
gray为灰度图像，
blockSize为邻域窗口的大小，
ksize是用于Soble算子的参数，
k是一个常量，取值为0.04~0.06
"""

filename = 'test.jpg'

img = cv2.imread(filename)
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
gray = np.float32(gray)
dst = cv2.cornerHarris(gray, blockSize, ksize, k)

Harris角点检测算法实现流程
- 计算图像$I(x, y)$在X和Y两个方向上的梯度$I_x$和$I_y$
- 计算图像两个方向梯度的乘积
- 使用高斯函数对$I_x^2$、$I_y^2$和$I_{xy}$进行高斯加权，生成矩阵$M$
- 计算Harris响应值R
- 给定阈值，当一个位置的强度大于阈值则认为是角点
- 在3×3或5×5的邻域内进行非最大值抑制，局部最大值点即为图像中的角点
  - 在一个窗口内，如果有多个角点则用值最大的那个角点，其他的角点都删除
实现代码

from scipy.ndimage import filters
import cv2
from matplotlib.pylab import *

class Harris(object):
    def __init__(self, img_path):
        self.img = cv2.imread(img_path, 0)

    def calculate_corner_strength(self):
        # 计算图像在x、y方向的梯度
        # 用滤波器采用差分求梯度的方式
        scale = self.scale
        k = self.k
        img = self.img
        gradient_imx, gradient_imy = zeros(img.shape), zeros(img.shape)
        filters.gaussian_filter(img, (scale, scale), (0, 1), gradient_imx)
        filters.gaussian_filter(img, (scale, scale), (1, 0), gradient_imy)

        # 计算矩阵M的每个分量
        I_xx = filters.gaussian_filter(gradient_imx*gradient_imx, scale)
        I_xy = filters.gaussian_filter(gradient_imx*gradient_imy, scale)
        I_yy = filters.gaussian_filter(gradient_imy*gradient_imy, scale)

        # 计算矩阵的迹、特征值和响应强度
        det_M = I_xx * I_yy - I_xy ** 2
        trace_M = I_xx + I_yy
        return det_M + k * trace_M ** 2

    def corner_detect(self, img):
        # 首先对图像进行阈值处理
        _threshold = img.max() * self.threshold
        threshold_img = (img > _threshold) * 1
        coords = array(threshold_img.nonzero()).T
        candidate_values = [img[c[0], c[1]] for c in coords]
        index = argsort(candidate_values)

        # 选取领域空间，如果邻域空间距离小于min的则只选取一个角点
        # 防止角点过于密集
        neighbor = zeros(img.shape)
        neighbor[self.min:-self.min, self.min:-self.min] = 1
        filtered_coords = []
        for i in index:
            if neighbor[coords[i, 0], coords[i, 1]] == 1:
                filtered_coords.append(coords[i])
                neighbor[(coords[i, 0] - self.min):(coords[i, 0] + self.min),
                (coords[i, 1] - self.min):(coords[i, 1] + self.min)] = 0
        return filtered_coords

    def corner_plot(self, img, corner_img):
        figure()
        gray()
        imshow(img)
        plot([p[1] for p in corner_img], [p[0] for p in corner_img], 'ro')
        axis('off')
        show()
    
    def __call__(self, k=0.04, scale=3, min=15, threshold=0.03):
        self.k = k
        self.scale = scale
        self.min = min
        self.threshold = threshold
        strength_img = self.calculate_corner_strength()
        corner_img = self.corner_detect(strength_img)
        self.corner_plot(self.img, corner_img)

if __name__ == '__main__':
    harris = Harris("2007_002619.jpg")
    harris()

4. 特性分析

Harris角点检测算子对亮度和对比度的变化不敏感
- 在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感，由于阈值的选择，可能会影响角点检测的数量
Harris角点检测算子具有旋转不变性
- Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值𝑅R也不发生变化，由此说明Harris角点检测算子具有旋转不变性
Harris角点检测算子不具有尺度不变性
- 后续介绍Harris角点解决尺度变换场景

5. 多尺度Harris角点

Harris角点检测算子与高斯尺度空间表示相结合，使用Harris角点检测算子具有尺度的不变性

vSLAMNet（三）-特征点-Harris角点