Harris角点检测、[公式推导+代码实现]

Harris角点检测

提示：有添加超链接，点击跳转

视频推荐：小破站，会飞的吴克：点击跳转

opencv：简述
强烈推荐：CSE486, Penn State

---

提示：写完文章后，目录可以自动生成，如何生成可参考右边的帮助文档

---

提示：以下是本篇文章正文内容，下面案例可供参考

一、基本原理

左侧图，窗口W在平坦区域的运动，当位移比较小时，没有差值
中间图，窗口W只有在edge垂直方向，会有较大差值
右侧图，窗口W向任何方向移动，都会有较大的差值

公式：
$$\begin{gathered} E(u,v)=\sum _{(x,y)ϵW}w(x,y)[I(x+u,y+v)-I(x,y){]}^2（公式1） \\ 对I(x+u,y+v)进行泰勒展开。 \\ I(x+u,y+v)=I(x,y)+\frac{\partial I}{\partial x}u+\frac{\partial I}{\partial y}u+o() \\ I(x+u,y+v)\approx I(x,y)+\frac{\partial I}{\partial x}u+\frac{\partial I}{\partial y}u\approx I(x,y)+\left[\begin{array}{cc}{I}_x& I_y\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right](公式3) \end{gathered}$$
将公式3代入公式1,得到：
$$\begin{gathered} E(u,v)\approx \sum _{(x,y)ϵW}w(x,y)[\left[\begin{array}{cc}{I}_x& I_y\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right]{]}^2 \\ =\sum _{(x,y)ϵW}w(x,y)\left[\begin{array}{cc}u& v\end{array}\right]\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_y{I}_x& I_y^2\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right](公式4) \end{gathered}$$
对公式4接着进行展开，实对称矩阵相似对角化
$$\begin{gathered} E(u,v)=\left[\begin{array}{cc}u& v\end{array}\right]\sum w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_y{I}_x& I_y^2\end{array}\right]\left[\begin{array}{c}u\\ v\end{array}\right] \\ =\left[\begin{array}{cc}u& v\end{array}\right]P\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]P^T{\left[\begin{array}{cc}u& v\end{array}\right]}^T \\ =\left[\begin{array}{cc}{u}^{\prime }& v^{\prime }\end{array}\right]\left[\begin{array}{cc}{\lambda }_1& 0\\ 0& {\lambda }_2\end{array}\right]{\left[\begin{array}{cc}{u}^{\prime }& v^{\prime }\end{array}\right]}^T \\ ={\lambda }_1{u}^{\prime 2}+{\lambda }_2{v}^{\prime 2}上式展开 \\ =\frac{u^{\prime 2}}{\frac{1}{{\lambda }_1}}+\frac{v^{\prime 2}}{\frac{1}{{\lambda }_2}}（公式5） \end{gathered}$$
由公式4可知，当中间的二维矩阵（M）特征值最大时，即得到角点。
由公式5，就是一个椭圆公式，当椭圆成圆的时候面积是最大的，即偏差E(u,v)最大。

$$M=\sum w(x,y)\left[\begin{array}{cc}{I}_x^2& I_x{I}_y\\ I_y{I}_x& I_y^2\end{array}\right]$$

如何更方便的表示E最大呢？

Corner Response Measure
$$\begin{gathered} R=det(M)-k[trace(M){]}^2 \\ det(M)={\lambda }_1\ast {\lambda }_2 \\ trace(M)={\lambda }_1+{\lambda }_2 \end{gathered}$$
k为一个经验值，k = 0.04~0.06
所以，只用求得矩阵M的行列式和迹，便可以表示R

最后，在对R进行非极大值抑制，便得到key point

二、代码实现

1.自己实现的高斯函数和卷积

def gaussian_kernel_2d(size, sigma_x, sigma_y=None):
    if sigma_y is None:
        sigma_y = sigma_x

    x, y = np.meshgrid(np.linspace(-(size//2), size//2, size),
                       np.linspace(-(size//2), size//2, size))

    kernel_2d = (1 / (2 * np.pi * sigma_x * sigma_y)) * np.exp(-(x ** 2 / (2 * sigma_x ** 2) + y ** 2 / (2 * sigma_y ** 2)))
    # 归一化
    kernel_2d /= kernel_2d.sum()

    return kernel_2d
    
def conv2d(image, kernel):
    image = np.array(image, dtype=np.float64)
    kernel = np.array(kernel, dtype=np.float64)

    im_height, im_width = image.shape
    ker_height, ker_width = kernel.shape

    padded_image = np.zeros((im_height + ker_height - 1, im_width + ker_width - 1))
    padded_image[:im_height, :im_width] = image

    result = np.zeros_like(image)

    for i in range(im_height):
        for j in range(im_width):
            # 计算当前像素周围与内核大小一致区域内的卷积结果
            window = padded_image[i:i+ker_height, j:j+ker_width]
            conv_result = np.sum(window * kernel)
            result[i, j] = conv_result

    return result

2.计算R, Response Measure

代码如下（示例）：

def harris_response(im, k=0.04, window_size=3, threshold=0.1):
    im.astype(np.float32)

    x_kernel = np.array([[-1, 0, 1],
                         [-2, 0, 2],
                         [-1, 0, 1]])
    y_kernel = np.array([[-1, -2, -1],
                         [0, 0, 0],
                         [1, 2, 1]])
    gaussian_k = gaussian_kernel_2d(3, 2)

    Ix = conv2d(im, x_kernel)
    Iy = conv2d(im, y_kernel)

    Ixx = Ix * Ix
    Iyy = Iy * Iy
    Ixy = Ix * Iy
    Ixx = conv2d(Ixx, gaussian_k)
    Iyy = conv2d(Iyy, gaussian_k)
    Ixy = conv2d(Ixy, gaussian_k)

    R_ = np.zeros(im.shape)

    for i in range(im.shape[0]):
        for j in range(im.shape[1]):
            M = np.array([[Ixx[i, j], Ixy[i, j]], [Ixy[i, j], Iyy[i, j]]])
            R_[i, j] = np.linalg.det(M) - k * np.trace(M) ** 2
    return R_

3. main

import cv2
import numpy as np

if __name__ == '__main__':
    in_filename = "./126007.jpg"
    im = cv2.imread(in_filename,cv2.IMREAD_ANYCOLOR)
    R = harris_response(cv2.cvtColor(im, cv2.COLOR_BGR2GRAY))
	print(R.mean())
    im[R > R.mean()] = [0,0,255]
    cv2.imwrite("harris.jpg", im)

原图

im[R > R.mean()] = [0,0,255] 效果图：

这里先简单的im[R > 10*R.mean()] = [0,0,255] 进行处理。简单粗暴，harris基本上掌握了～

也可以通过归一化、滤波器对相应函数R进行处理，再进行非极大值抑制，减少冗余检测点。

---

参考

提示：这里对文章进行总结：

https://blog.csdn.net/LeoEdwin/article/details/127539497
https://www.cse.psu.edu/~rtc12/CSE486/lecture06.pdf