Harris角点检测

已于 2023-05-27 22:18:33 修改
阅读量373 收藏 3
点赞数 3
分类专栏： 智能图像分析 文章标签： opencv 计算机视觉 人工智能
于 2023-04-26 21:47:30 首次发布
版权声明：本文为博主原创文章，遵循 CC 4.0 BY-SA 版权协议，转载请附上原文出处链接和本声明。
本文链接：https://blog.csdn.net/SharbbyBai/article/details/130381728
版权

一、原理

  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 几秒钟写的代码效果好。希望继续努力吧。

一个没有本领的人
关注
  • 3
    点赞
  • 踩
  • 3
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
专栏目录
【老生谈算法】Harris角点检测算法原理及其MATLAB编程实现.doc
07-02
【老生谈算法】Harris角点检测算法原理及其MATLAB编程实现.doc
Harris角点检测算法
10-02
harris角点检测 算法代码和说明文档。采用matlab编写的。压缩文件里面包含m文件，dox说明为文档，ppt说明文档，以及检测用的图片
Harris角点检测算法示例
10-29
使用OPENCV+VS2013实现Harris角点检测，可以准确检测出图像中的Harris角点
基于Harris角点检测的改进算法研究
10-16
经典的Harris算法在提取图像的角点上具有计算简单、适应性强等优势，但该方法由于人为设定单一阈值，容易出现伪角点、漏检点及运行速度不理想等现象。针对这一情况，文章在传统Harris算法基础上提出一种新的检测方法，采用多阈值的圆形非极大值抑制法提取角点, 以此降低算法检测时间并增强图像旋转不变性，再借鉴SUSAN思想消去大部分伪角点。通过实验对比，该算法具有更好的角点检测性，为后期的图像配准奠定了良好的基础。
Harris角点检测（OpenCV）
冰冰bing的博客
10-31 1486
文章目录声明正文一、关于角点二、Harris角点检测原理三、opencv中的Harris角点检测四、亚像素级精确度的角点五、Shi-Tomasi角点检测感谢观看！ 声明 声明：本系列博客是我在学习OpenCV官方教程中文版(For Python)(段力辉 译)所做的笔记。所以，其中的绝大部分内容引自这本书，博客中的代码也是其配套所附带的代码或书中的代码，侵删。其中部分代码可能会因需要而改动。在本系列博客中，其中包含书中的引用，也包括我自己对知识的理解，思考和总结。本系列博客的目的主要有两个，一个是可以作为
OpenCV入门学习笔记之Harris角点检测与SIFT特征匹配算法
Miracle8070
04-05 9137
1. 写在前面 这篇文章整理两个图像处理中非常重要的算法，一个是Harris角点检测算法，另一个是SIFT特征匹配算法，这两个算法本质上还是去找图像里面的关键特征点，帮助我们后续更好的理解图像以及做各种各样的分析。 由于这两个算法涉及到的数学原理会比较多，而我刚入门，所以只是从使用的角度，简单的描述到底在做什么事情，至于详细的数学细节或者推导，这里不过多整理，以掉包能完成任务为首要目的啦。 首先，先介绍Harris角点检测算法，角点在图像中是很重要的特征，信息含量很高，那么如何找到一个图像里面的角点呢？
Harris角点检测数学计算过程与CornerHarris方法参数的一点说明【DataWhale学习记录】
m0_38052500的博客
06-24 3801
1 基础知识 1.1 图像梯度的解释 图像的梯度与数学中的梯度在形式上是由差异的，原因在于图像的特点，图像是一个离散的二维函数，接下来会尝试说明这一点，但从意义上来说是相同的，是为了表现图像灰度的变化率。 在数学微积分中，一维函数的一阶微分的基本定义是： 二维函数的一阶微分的基本定义是： 对于灰度图这样二维数组的图像，它其实就是一个离散的二位函数，说它离散，是因为每个灰度值取值为整数而不是小数，因此ϵ不能无限小，而ϵ的最小单位即是1像素。因此当ϵ取最小值1时，灰度图（离散的二维函数）的一阶微分基本定义
C语言实现harris角点检测
xiaoyang0307的博客
01-19 3736
纯C语言实现harris角点检测，代码忽略高斯滤波部分 int sobel_kernel_x[9]={-1,0,1,-2,0,2,-1,0,1}; int sobel_kernel_y[9]={1,2,1,0,0,0,-1,-2,-1}; extern uchar sobel_res_x[1000*1000]; extern uchar sobel_res_y[1000*1000]; uchar xx[1000*1000]; uchar yy[1000*1000]; uchar xy[1000
Harris角点检测原理详解
Rolandxxx的博客
08-21 4844
如果是平坦区域，那么像素点的梯度所构成的点集比较集中在原点附近，因为窗口区域内的像素点的梯度幅值非常小，此时矩阵M的对角化的两个特征值比较小；如果是边缘区域，在计算像素点的x、y方向上的梯度时，边缘上的像素点的某个方向的梯度幅值变化比较明显，另一个方向上的梯度幅值变化较弱，其余部分的点都还是集中原点附近，这样M对角化后的两个特征值理论应该是一个比较大，一个比较小，当然对于边缘这种情况，可能是呈45°的边缘，致使计算出的特征值并不是都特别的大，总之跟含有角点的窗口的分布情况还是不同的。
计算机视觉（角点检测）- 1 - Harris角点检测
热门推荐
qq_40374812的博客
05-18 1万+
计算机视觉（角点检测）-1.1 Harris角点检测学习前言 学习前言   今天更新一下计算机视觉的东西吧，哈哈哈o(￣︶￣)o   关于角点检测呢，就是检测图像中角点的算法，常见的有四种常见的算法：Moravec角点检测算法、 Harris角点检测算法、Shi-Tomasi 算法、FAST角点检测算法，不过常用的就两种吧——Harris角点检测算法、Shi-Tomasi 算法（反正我基本上就只用到这两种，先着重介绍前面两种，后面再花世界介绍吧）(～￣▽￣)～ ...
系统性综述：特征点检测与匹配
小白学视觉
02-10 1233
点击上方“小白学视觉”，选择加"星标"或“置顶” 重磅干货，第一时间送达 作者丨梦寐mayshine@知乎来源丨https://zhuanlan.zhihu.com/p/1...
Harris角点检测python实现及基于opencv实现
xiaohuolong的博客
07-01 5468
Harris角点检测python实现及基于opencv实现，（角点检测、非极大值抑制）
基于Python手动实现Harris角点检测
weixin_46585836的博客
12-01 3856
基于Python手动实现Harris角点检测算法
【角点检测】基于 Harris算法实现角点特征检测含Matlab源码
qq_59747472的博客
03-13 1183
1 简介 Harris 算子是一种计算简单、应用广泛的角点检测算子，只使用灰度的一阶差分和滤波，可以定量地提取特征角点并且提取的角点特征均匀。由于计算过程用到图像的一阶导数，故对存在灰度变化、图像旋转、视点变换和噪声干扰的图像也能稳定地提取角点。因此本文提出的算法采用 Harris 角点作为兴趣点，并提取每个兴趣点为中心的邻域内空域的五个均值特征形成特征向量，通过位置转移向量的频率统计连接标志匹配点。 2 完整代码 clear;​filename = 'Lena.jpg';X = imread(
OpenCV 双目相机标定
最新发布
dayuhaitang1的博客
06-08 242
OpenCV 双目相机标定
harris角点检测
11-28
Harris角点检测是一种常用的图像处理技术，用于检测图像中的角点。Harris角点检测算法通过计算图像中像素的局部灰度变化来识别角点，从而找到图像中的特征点。这种方法首先计算出每一个像素周围局部区域的灰度变化...

“相关推荐”对你有帮助么？

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论
添加红包

请填写红包祝福语或标题

个

红包个数最小为10个

元

红包金额最低5元

当前余额3.43元 前往充值 >
需支付：10.00元
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付元
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明：

1.余额是钱包充值的虚拟货币，按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载，可以购买VIP、付费专栏及课程。

余额充值