Harris角点检测

一、原理

1. 请见这篇博客

2. Harris角点检测流程为：

① 计算图像所有像素的梯度Ix，Iy；
②计算Ix^2 ，Iy^2，IxIy
③使用窗函数对Ix^2 ，Iy^2，IxIy滤波，生成矩阵M的元素A,B,C

④计算每个像素的Harris响应值R，并对小于某阈值的R置零

⑤在局部小邻域内进行非极大值抑制，局部最大值即为角点

二、代码

（一）调用opencv库函数

import cv2

img = cv2.imread("D:/home.jpg")
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

# Harris 角点检测
dst = cv2.cornerHarris(gray, 2, 3, 0.04)

# 角点标记
dst = cv2.dilate(dst, None)
img[dst > 0.04 * dst.max()] = [0, 0, 255]

cv2.imshow("Harris Corners", img)
cv2.waitKey(0)
cv2.destroyAllWindows()

检测结果如下图所示：

dst是高为384，宽为512的数组

（二）不调用库函数

参考了这位博主的代码

检测结果如下：

输出的Ix、Iy如下图所示：

（三）参考源码编写

1.计算像素梯度Ix、Iy

本次实验使用的图片大小为 W * H = 512 * 384.
可以看到计算Ix、Iy调用的是cv2.sobe()函数，源码请见这篇博客。在使用sobel算子进行梯度求解时涉及到一个问题，就是我们用 3 * 3 的卷积核进行卷积时有一个像素的边缘没法被卷积操作，使用 5 * 5 的卷积核时有两个像素的边缘没法进行卷积操作（请见这篇博客）而cv2.sobe()函数输出的梯度矩阵Ix、Iy却是与图片一样大小的矩阵，因此我们需要对图片进行边界填充。
边界填充有四种方式，详细请见这篇博客。
sobel算子中使用的卷积函数为sepFilter2D()，sepFilter2D()采用的边界填充方式为 BORDER_DEFAULT，即采用镜像填充。我将sobel函数写为如下形式：

# 利用Sobel函数计算图像梯度
def sobel(image):
    Gx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    Gy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])

    # 对图像进行边界填充
    img_border = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    H, W = img_border.shape
    Grx = np.zeros((H-2, W-2), np.int)
    Gry = np.zeros((H-2, W-2), np.int)

    for i in range(H-2):
        for j in range(W-2):
            dx = np.sum(img_border[i:i + 3, j:j + 3] * Gx)
            Grx[i, j] = np.abs(dx)
            dy = np.sum(img_border[i:i + 3, j:j + 3] * Gy)
            Gry[i, j] = np.abs(dy)

    return Grx, Gry

Ix, Iy = sobel(img_gray)

输出的Ix、Iy如下图所示：

显然与调用库函数的结果不同。
又仔细阅读了cv2.sobe()函数的讲解，发现在这个函数中超出0~255范围的数会被截断（截断指的是，当数值大于255保留为255，当数值小于0保留为0，其余不变），详情请见这篇博客。
接下来我将函数改为如下形式：

# 利用Sobel函数计算图像梯度
def sobel(image):
    Gx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    Gy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])

    # 对图像进行边界填充
    img_border = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    H, W = img_border.shape
    Grx = np.zeros((H-2, W-2), np.int)
    Gry = np.zeros((H-2, W-2), np.int)

    for i in range(H-2):
        for j in range(W-2):
            dx = np.sum(img_border[i:i + 3, j:j + 3] * Gx)
            if dx < 0:
                dx = 0
            Grx[i, j] = dx
            dy = np.sum(img_border[i:i + 3, j:j + 3] * Gy)
            if dy < 0:
                dy = 0
            Gry[i, j] = dy

    return Grx, Gry

Ix, Iy = sobel(img_gray)

输出矩阵正确。

2.计算Ix^2 ，Iy^2，IxIy

实际就是两个矩阵对应位置相乘。

3.高斯滤波，生成矩阵M

这里用到的高斯滤波器实际上就是替代了窗函数 w(x, y) ，用来对一个窗口中的灰度变化值进行滤波，以便我们下一步对窗口中心的像素点是否是角点进行判断。
在（二）中的代码中，用到的是cv2.GaussianBlur()函数，关于这个函数的讲解请看这篇博客，讲得非常详细。
总的来说，首先就是通过getGaussianKernel()函数生成两个一维高斯核kernelX和kernelY，然后传入 cv2.sepFilter2D函数，是的又用到了这个函数，这个函数首先对镜像边缘填充之后图像的每一行以kernelX为卷积核做卷积，再对结果的每一列以kernelY为卷积核做卷积，最后归一化。这里我们不这么麻烦，直接指定一个3*3高斯核进行卷积。
对于这里的高斯窗口，我个人有一些理解，不知道对不对。窗口越大，所包含的像素点就越多，对于正常的窗口来说，每个像素点前面的系数都是1，算出来的灰度值差就是实际的灰度值差，然而我们在这里只是对窗口中心像素点是否为角点进行判断，故可以对中心像素点的权值设置的大一些，对周围像素点的权值设置的小一些，这样在窗口移动时，该点灰度值的变化就能贡献大一些。（link）

# 使用高斯平滑滤波进行加权计算
def GaussianBlur(array, K_size=3, sigma=1.3):
    # 对数组进行边界填充
    array_border = cv2.copyMakeBorder(array, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    H, W = array_border.shape
    array_Gaussian = np.zeros((H - 2, W - 2), np.int)

    # 准备卷积核
    pad = K_size // 2
    kenel = np.zeros((K_size, K_size), dtype=np.float)
    for x in range(K_size):
        for y in range(K_size):
            kenel[x, y] = np.exp(-(x ** 2 + y ** 2) / (2 * (sigma ** 2)))
    kenel /= (2 * np.pi * sigma * sigma)
    kenel /= kenel.sum()  # 归一化

    # 卷积
    for x in range(H-2*pad):
        for y in range(W-2*pad):
            array_Gaussian[x, y] = int(np.sum(kenel * array_border[x: x + K_size, y: y + K_size]))

    return array_Gaussian

Ix2 = GaussianBlur(Ix2)
Ixy = GaussianBlur(Ixy)
Iy2 = GaussianBlur(Iy2)

4.计算R值，并在局部小邻域内进行非极大值抑制

① 原本将代码写为如下形式：

# 计算R值
def R_calculate(array1, array2, array3, k = 0.04):
    H, W = array1.shape
    R = np.zeros((H, W))  # 定义空的R矩阵
    for i in range(H):
        for j in range(W):
            R[i, j] = (array1[i, j]*array2[i, j]) - (array3[i, j]*array3[i, j]) - k * ((array1[i, j] + array2[i, j])**2)
    return R

R = R_calculate(Ix2, Iy2, Ixy)
out = img[R > 0.04 * R.max()] = [0, 0, 255]

然后程序有警告：

D:\project\intelligent image analysis\Harris.py:88: RuntimeWarning: overflow encountered in long_scalars
  R[i, j] = (array1[i, j] * array2[i, j]) - (array3[i, j] * array3[i, j]) - k * ((array1[i, j] + array2[i, j])**2)
Traceback (most recent call last):
  File "D:\project\intelligent image analysis\Harris.py", line 94, in <module>
    cv2.imwrite('D:/project/intelligent image analysis/data/out.png', out)
cv2.error: OpenCV(4.6.0) :-1: error: (-5:Bad argument) in function 'imwrite'
> Overload resolution failed:
>  - img is not a numpy array, neither a scalar
>  - Expected Ptr<cv::UMat> for argument 'img'

说是 R[i, j] 这里有溢出。另外out那里表示有误。

② 我又将代码改为如下形式：

# 计算R值
def R_calculate(array1, array2, array3, k = 0.04):
    H, W = array1.shape
    R = np.zeros((H, W))  # 定义空的R矩阵
    A = np.zeros((H, W))
    B = np.zeros((H, W))
    C = np.zeros((H, W))
    for i in range(H):
        for j in range(W):
            A[i, j] = array1[i, j]*array2[i, j]  # 即M矩阵中的AB
            B[i, j] = array3[i, j]*array3[i, j]  # 即M矩阵中的C^2
            C[i, j] = (array1[i, j] + array2[i, j])**2  # 即M矩阵中的（A+B）^2
            R[i, j] = A[i, j] - B[i, j] - k * C[i, j]
    return R

R = R_calculate(Ix2, Iy2, Ixy)
img[R > 0.04 * R.max()] = [0, 0, 255]

还是有警告：

③ 继续改：

def R_calculate(array1, array2, array3, k = 0.04):
    H, W = array1.shape
    R = np.zeros((H, W))  # 定义空的R矩阵
    A = np.zeros((H, W), np.int)
    B = np.zeros((H, W), np.int)
    C = np.zeros((H, W), np.int)
    for i in range(H):
        for j in range(W):
            A[i, j] = array1[i, j]*array2[i, j]  # 即M矩阵中的AB
            B[i, j] = array3[i, j]**2  # 即M矩阵中的C^2
            C[i, j] = (array1[i, j] + array2[i, j])**2  # 即M矩阵中的（A+B）^2
            R[i, j] = A[i, j] - B[i, j] - k * C[i, j]
    return R

这里将 B 的表达式改了一下，B 的警告没有了：

④ 继续改：

def R_calculate(array1, array2, array3, k = 0.04):
    H, W = array1.shape
    R = np.zeros((H, W))  # 定义空的R矩阵
    A = np.zeros((H, W), np.float)
    B = np.zeros((H, W), np.float)
    C = np.zeros((H, W), np.float)
    for i in range(H):
        for j in range(W):
            A[i, j] = array1[i, j]*array2[i, j]  # 即M矩阵中的AB
            B[i, j] = array3[i, j]**2  # 即M矩阵中的C^2
            C[i, j] = (array1[i, j] + array2[i, j])**2  # 即M矩阵中的（A+B）^2
            R[i, j] = A[i, j] - B[i, j] - k * C[i, j]
    return R

这里将A、B、C定义为 float 型的数组，R 的警告消失了：

⑤ 继续：
对 A 的警告实在不知道怎么改，所以最后调用了个函数。。。。。。

def R_calculate(array1, array2, array3, k = 0.04):
    H, W = array1.shape
    R = np.zeros((H, W))  # 定义空的R矩阵
    A = np.multiply(array1, array2)
    B = np.zeros((H, W), np.float)
    C = np.zeros((H, W), np.float)
    for i in range(H):
        for j in range(W):
            B[i, j] = array3[i, j]**2  # 即M矩阵中的C^2
            C[i, j] = (array1[i, j] + array2[i, j])**2  # 即M矩阵中的（A+B）^2
            R[i, j] = A[i, j] - B[i, j] - k * C[i, j]
    return R

没有警告了。

5.检测结果

6.画圈

为了将角点用圆圈表示出来，我在代码最后做了如下更改：

R = r_calculate(Ix2, Iy2, Ixy)
# cv2.circle(img, img[R > 0.04 * R.max()], 5, (0, 0, 255), 1)
for t in range(h):
    for q in range(w):
        if R[t, q] > c:
            cv2.circle(img, (t, q), 5, (0, 0, 255), 1)

然而检测结果是这样的：

可以看到全部都乱了。

对下面这行代码进行分析：

img[R > 0.04 * R.max()] = [0, 0, 255]

（1）令：

a = R.max()
c = 0.04*a
b = R > 0.04 * R.max()

结果如下：


b是一个 384 * 512 大小的矩阵，在满足 R > 0.04 * R.max() 的条件的位置值为 truth，其他时候值为 false。
（2）将代码改为如下形式：

for t in range(h):
    for q in range(w):
        if R[t, q] > c:
            cv2.circle(img, (t, q), 5, (0, 0, 255), 1)
            out = img
            cv2.imwrite('D:/out.png', out)

进行分步运行后，发现在值为 truth 的位置根本没有画圈。

（3）

for t in range(h):
    for q in range(w):
        if R[t, q] > c:
            img[t, q] = [0, 0, 255]
            cv2.circle(img, (t, q), 5, (0, 0, 255), 1)
            out = img
            cv2.imwrite('D:out.png', out)

然而在这种情况下，值为 truth 的位置是被标成了红色的。说明是 cv2.circle() 函数出了问题。
（4）在查了资料之后，有些博主说要对图片进行如下处理：

img = cv2.cvtColor(img, cv2.COLOR_RGB2BGR)

照做之后还是不对。
(5)然后，我就看到了这篇帖子：
cv2.circle() 中圆心的（x，y）和原图中的（x，y）是反着的！
改了就对了：

R = r_calculate(Ix2, Iy2, Ixy)
for t in range(h):
    for q in range(w):
        if R[t, q] > 0.3 * R.max():
            cv2.circle(img, (q, t), 3, (0, 0, 255), 1, lineType=3)

有点乱，我将参数又改了一下，效果如下图所示：

示例效果如下：

效果比示例效果要差很多。而且。。。。。我传到实验平台上去之后发现代码运行时间超时了。。。。可见循环真的很消耗时间。。。于是我用ChitGPT写了一个代码，参考了它后面的画圈部分（可见它的画圈坐标也是 x 和 y 反着的）：

R = r_calculate(Ix2, Iy2, Ixy)

corners = np.argwhere(R > 0.3 * R.max())
for corner in corners:
    x, y = corner
    cv2.circle(img, (y, x), 3, (0, 0, 255), 1)

其中 corners 保存的是满足 R > 0.3 * R.max() 像素点的坐标。

7.改进

ChitGPT代码如下：

import numpy as np
import cv2


# 读取图像
img = cv2.imread('D:/home.jpg')


# 转化为灰度图像
gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)


# 设置窗口大小和卷积核大小
window_size = 3
kernel_size = 3


# 计算图像梯度
dx = cv2.Sobel(gray, cv2.CV_64F, 1, 0, ksize=kernel_size)
dy = cv2.Sobel(gray, cv2.CV_64F, 0, 1, ksize=kernel_size)


# 计算Harris矩阵的三个分量
Ixx = dx**2
Ixy = dx*dy
Iyy = dy**2


# 对Harris矩阵的三个分量进行高斯滤波
ksize = (window_size, window_size)
sigma = 1.5
Ixx = cv2.GaussianBlur(Ixx, ksize, sigma)
Ixy = cv2.GaussianBlur(Ixy, ksize, sigma)
Iyy = cv2.GaussianBlur(Iyy, ksize, sigma)


# 计算Harris响应函数
k = 0.04
det = Ixx*Iyy - Ixy**2
trace = Ixx + Iyy
R = det - k*trace**2


# 设置阈值并进行非极大值抑制
threshold = 0.05*np.max(R)
R[R < threshold] = 0
R_nms = cv2.dilate(R, np.ones((3, 3), np.uint8))
R_nms[R < R_nms] = 0


# 显示角点
corners = np.argwhere(R_nms > 0)
for corner in corners:
    x, y = corner
    cv2.circle(img, (y, x), 3, (0, 0, 255), 1)


# 显示图像
cv2.imshow('image', img)
cv2.waitKey(0)
cv2.destroyAllWindows()

它的检测结果如下图所示：

虽然还是没有示例结果效果好，但显然比我的好太多。我和它的主要区别在于 sobel 算子和高斯滤波函数，这两个函数我也是手写的，而ChitGPT是调用的 opencv 中的库函数。

（1）调用转灰度图函数

img_gray = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)

只调用一个函数的效果图如下图所示（左边为全部手写效果图，右边为调用灰度转换图的效果图）：

可以见到效果差别不大，比较明显的是窗户楞中间多检测出了一个圈，然而这个圈并不是角点。。。应该算是误检，所以手写的灰度图转换函数应该对结果影响不大。

（2）调用Sobel算子

只将sobel算子改为直接调用库函数：

kernel_size = 3
Ix = cv2.Sobel(img_gray, cv2.CV_64F, 1, 0, ksize=kernel_size)
Iy = cv2.Sobel(img_gray, cv2.CV_64F, 0, 1, ksize=kernel_size)

更改这个函数后检测效果如下图所示（左边为原效果，右边为调用后的效果）：

可以看到这个对比比较明显，显然调用库函数多检测出了很多角点。窗口那里虽然漏掉了原来那个点，但是又多检测出了三个角点。
对比 ChitGPT 调用的函数和前面的博主调用的函数，我们发现二者设置的图像深度不同。 ChitGPT 设置的图像深度为 CV_64F （ float64 ），而前面的博主所设置的 ddepth 值为 -1，即与原图像保持相同的图像精度。原图像像素精度都为整数型。故二者的输出梯度矩阵也有所不同。
ChitGPT 输出的梯度矩阵如下图所示：

均为float64类型；
我的输出矩阵与博主相同，如下图所示：

还是我之前特意更改的。那么我再对sobel部分的代码进行更改，即不对像素点进行截断，并使输出梯度矩阵为float类型：

# 利用Sobel函数计算图像梯度
def sobel(image):
    """

    :param image:
    :return:
    """
    Gx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    Gy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])

    # 对图像进行边界填充
    img_border = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    H, W = img_border.shape
    Grx = np.zeros((H - 2, W - 2), float)
    Gry = np.zeros((H - 2, W - 2), float)

    for i in range(H - 2):
        for j in range(W - 2):
            dx = np.sum(img_border[i:i + 3, j:j + 3] * Gx)
            Grx[i, j] = dx
            dy = np.sum(img_border[i:i + 3, j:j + 3] * Gy)
            Gry[i, j] = dy

    return Grx, Gry

输出矩阵如下图所示：

可见与 ChitGPT 输出矩阵基本相同。
再次运行，结果如下图所示（左边为原效果图，右边为更改梯度矩阵精度后的效果图）：

可见效果好了很多。再与 ChitGPT 的生成效果图进行比较（左边 ChitGPT 效果图，右边我的效果图）：

可见还是 ChitGPT 效果更好，基本没有重复检测。

（3）调用高斯滤波函数

在前面改了sobel算子的基础上，将高斯滤波函数改为调用库函数：

window_size = 3
ksize = (window_size, window_size)
sigma = 1.5
Ixx = cv2.GaussianBlur(Ix2, ksize, sigma)
Ixy = cv2.GaussianBlur(Ixy, ksize, sigma)
Iyy = cv2.GaussianBlur(Iy2, ksize, sigma)

输出结果如下图（左边为 ChitGPT 效果图，右边为我的效果图）：

可见差别非常大。
由于这个地方我们的输入数组是一样的，调用函数也是一样的，那么我们的输出数组也一定是一样的。造成结果差异大的原因只能是后面的Harris响应函数计算的不同。这个等一下再说。我们先来看一下高斯滤波函数对结果的影响（左边为我原来函数的效果图，右边为调用高斯滤波函数的效果图）：

这里我将窗口函数的大小都设置为了 3，sigma 都设置为了1.5 。显然差别还是很大。而且与上图中 ChitGPT 的效果图相比，检测结果的重复率还是很高。
来看一下我自己写的高斯函数和库函数的输出有什么不同，以Ix2为例，左边为我自己写的，右边为库函数的：

可见输出结果还是有些差别的。另外数据类型也不同。
又仔细阅读了前面链接中博主所介绍的高斯滤波函数，发现我的滤波函数是有一些问题。直接放上 k-size = 3, sigma = 1.5 时的高斯卷积核：
权重矩阵如下：

归一化的权重矩阵如下：

由高斯公式可知卷积核权重矩阵是轴对称也是对角线对称的。而我的卷积核如下图所示：

只是对角线对称。然后我发现我的卷积核是这样配置的：

pad = k_size // 2
    kernel = np.zeros((k_size, k_size), dtype=float)
    for x in range(k_size):
        for y in range(k_size):
            kernel[x, y] = np.exp(-(x ** 2 + y ** 2) / (2 * (sigma ** 2)))
    kernel /= (2 * np.pi * sigma * sigma)
    kernel /= kernel.sum()  # 归一化

也就是说我的卷积核坐标左上角的点是（0，0），右下角的点是（2，2）。而真正的卷积核坐标应该如下图所示：

所以计算出来的数值都是不对的。然后我将代码改为如下形式：

# 准备卷积核
    pad = k_size // 2
    kernel = np.zeros((k_size, k_size), dtype=float)
    for x in range(-pad, -pad + k_size):
        for y in range(-pad, -pad + k_size):
            kernel[y + pad, x + pad] = np.exp(-(x ** 2 + y ** 2) / (2 * (sigma ** 2)))
    kernel /= (2 * np.pi * sigma * sigma)
    kernel /= kernel.sum()  # 归一化

这样卷积核左上角的坐标就变为了（-1，-1）。这样输出的 Ix2 矩阵就与调用函数输出的矩阵相同了。

对比效果如下图所示，左边为我原来的代码，中间为 ChitGPT 的代码，右边为更改sobel和高斯滤波函数之后的代码：

可见角点基本都检测出来了，但是还是有很多点都有重复检测。接下来再来对比一下Harris响应函数。

（4）Harris响应函数

对于重复检测这个点，我又回去看了一下 Harris 角点检测的流程，如下图：

我应该是因为没有进行第五步，即在局部小邻域内进行非极大值抑制。在 ChitGPT 的代码中，它采用了 cv2.dilate() 函数进行邻域内的非极大值抑制，即膨胀操作。
我添加的代码如下：

R[R < 0.1 * R.max()] = 0  # 对小于阈值的R置0——步骤①
R_nms = cv2.dilate(R, np.ones((3, 3), np.uint8))  # 利用膨胀操作，在邻域内进行非极大值抑制——步骤②
R_nms[R < 0.1 * R.max()] = 0  # 对小于阈值的R_nms置0——步骤③

然后我的运行结果是下面这样的：

不太对劲。
又去复习了一下膨胀。具体请看下面链接：
link1，link2
总之膨胀就是取局部最大值。以此处的 R 为例，在执行完步骤 ① 的小于阈值置 0 之后，输出的数组如下图所示：

执行完步骤 ② 膨胀之后的 R_nms 数组如下图所示：

值得注意的是膨胀操作是在原数组上进行的，而不是卷积核平移一次保存一次，否则就会膨胀起来没完没了。。。。
根据我个人的理解，这一步执行膨胀操作实际上是把一个小邻域内检测的角点近似看作是同一个点，为了筛选出检测的最准确的那个点（即 R 值最大的那个点），我们把小邻域内的所有 R 值都置为最大值，然后再通过与最大值进行比较过滤掉重复检测的点。
所以接下来的步骤 ③ 应该是过滤重复检测结果的过程，我上面写的步骤 ③ 并没有什么实际意义。正确的应该是如下这种形式：

R[R < 0.1 * R.max()] = 0  # 对小于阈值的R置0
R_nms = cv2.dilate(R, np.ones((3, 3), np.uint8))  # 利用膨胀操作，在邻域内进行非极大值抑制
R_nms[R < R_nms] = 0  # 对小于阈值的R_nms置0

过滤之后的数组 R_nms 结果如下图所示：

最终输出结果为：

可见比一开始的效果好了不少，虽然还是不如示例效果。

我将这份代码传到平台上，运行时间为11.45s。

重新看了一下代码，感觉中间有些重复的地方，而且我写的数组相乘与 np.multiply() 函数和 a * a 输出的结果是一样的，所以我把我写的 multiply 函数屏蔽掉了，直接改成了数组相乘的形式（即 Ix * Ix），去掉了几个 for 循环，另外计算 R 值的函数中的小循环也去掉了。完整代码如下所示：

import cv2
import numpy as np

img = cv2.imread('D:/home.jpg')  # 读图
########## Begin ##########
# 输入图像变为单通道
h, w = img.shape[:2]  # 获取图片的high和wide
print('h, w=', h, w)

img_gray = np.zeros([h, w], img.dtype)  # 创建一张和当前图片大小一样的单通道图片
for i in range(h):
    for j in range(w):
        m = img[i, j]
        img_gray[i, j] = int(m[0] * 0.11 + m[1] * 0.59 + m[2] * 0.3)  # 将BGR坐标转换为gray坐标

# 利用Sobel函数计算图像梯度
def sobel(image):
    """

    :param image:
    :return:
    """
    Gx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
    Gy = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]])

    # 对图像进行边界填充
    img_border = cv2.copyMakeBorder(image, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    H, W = img_border.shape
    Grx = np.zeros((H - 2, W - 2), float)
    Gry = np.zeros((H - 2, W - 2), float)

    for i in range(H - 2):
        for j in range(W - 2):
            dx = np.sum(img_border[i:i + 3, j:j + 3] * Gx)
            # if dx < 0:
            #     dx = 0
            Grx[i, j] = dx
            dy = np.sum(img_border[i:i + 3, j:j + 3] * Gy)
            # if dy < 0:
            #     dy = 0
            Gry[i, j] = dy

    return Grx, Gry


Ix, Iy = sobel(img_gray)

Ix2 = Ix * Ix
Ixy = Ix * Iy
Iy2 = Iy * Iy


# 使用高斯平滑滤波进行加权计算
def gaussianblur(array, k_size=3, sigma=1.5):
    # 对数组进行边界填充
    array_border = cv2.copyMakeBorder(array, 1, 1, 1, 1, cv2.BORDER_DEFAULT)
    H, W = array_border.shape
    array_Gaussian = np.zeros((H - 2, W - 2), float)

    # 准备卷积核
    pad = k_size // 2
    kernel = np.zeros((k_size, k_size), dtype=float)
    for x in range(-pad, -pad + k_size):
        for y in range(-pad, -pad + k_size):
            kernel[x + pad, y + pad] = np.exp(-(x ** 2 + y ** 2) / (2 * (sigma ** 2)))
    kernel /= (2 * np.pi * sigma * sigma)
    kernel /= kernel.sum()  # 归一化

    # 卷积
    for x in range(H - 2 * pad):
        for y in range(W - 2 * pad):
            array_Gaussian[x, y] = np.sum(kernel * array_border[x: x + k_size, y: y + k_size])

    return array_Gaussian


Ix2 = gaussianblur(Ix2)
Ixy = gaussianblur(Ixy)
Iy2 = gaussianblur(Iy2)

# 计算R值
def r_calculate(array1, array2, array3, k=0.04):
    A = array1 * array2  # 即M矩阵中的AB
    B = array3 ** 2  # 即M矩阵中的C^2
    C = (array1 + array2) ** 2  # 即M矩阵中的（A+B）^2
    R_array = A - B - k * C

    return R_array


R = r_calculate(Ix2, Iy2, Ixy)
# b = np.ones((3, 3), np.uint8)
R[R < 0.08 * R.max()] = 0  # 对小于阈值的R置0
R_nms = cv2.dilate(R, np.ones((3, 3), np.uint8))  # 利用膨胀操作，在邻域内进行非极大值抑制
R_nms[R < R_nms] = 0  # 对小于阈值的R_nms置0
corners = np.argwhere(R_nms > 0)  # 保存所有角点坐标
for corner in corners:
    x, y = corner
    cv2.circle(img, (y, x), 2, (0, 0, 255), 1)

out = img
########## End ##########
cv2.imwrite('D:/out.png', out)

传到平台上之后发现，运行时间为11.82s，变慢了。。。。
本人还是太菜了，python也不熟悉，很多函数都不知道，只能用循环来做，做了好几天的一个小实验也没有 ChitGPT 几秒钟写的代码效果好。希望继续努力吧。