我正在尝试从头开始实施 Harris 角点检测算法，但无法实施非最大抑制
I am trying to implement the Harris corner detection algorithm from scratch, but could not able to implement the non-maximum suppression
我正在尝试从头开始实施 Harris 角点检测算法。该算法的输出应该得到一个代表角的像素，但在我的代码中我得到了代表角的多个像素。这可能是因为未实现算法的最后部分 non-maximum suppression。这是我无法实施的，因为我没有正确理解它。如何实施？除此之外，我还试图找到这些角的坐标，如何在不使用 cv2 库的情况下做到这一点？
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as im
# 1. Before doing any operations convert the image into gray scale image
img = im.imread('OD6.jpg')
plt.imshow(img)
plt.show()
# split
R=img[:,:,0]
G=img[:,:,1]
B=img[:,:,2]
M,N=R.shape
gray_img=np.zeros((M,N), dtype=int);
for i in range(M):
for j in range(N):
gray_img[i, j]=(R[i, j]*0.2989)+(G[i, j]*0.5870)+(B[i, j]*0.114);
plt.imshow(gray_img, cmap='gray')
plt.show()
# 2. Applying sobel filter to get the gradients of the images in x and y directions
sobelx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype = np.float)
sobely = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype = np.float)
sobelxImage = np.zeros((M,N))
sobelyImage = np.zeros((M,N))
sobelGrad = np.zeros((M,N))
image = np.pad(gray_img, (1,1), 'edge')
for i in range(1, M-1):
for j in range(1, N-1):
gx = (sobelx[0][0] * image[i-1][j-1]) + (sobelx[0][1] * image[i-1][j]) + \
(sobelx[0][2] * image[i-1][j+1]) + (sobelx[1][0] * image[i][j-1]) + \
(sobelx[1][1] * image[i][j]) + (sobelx[1][2] * image[i][j+1]) + \
(sobelx[2][0] * image[i+1][j-1]) + (sobelx[2][1] * image[i+1][j]) + \
(sobelx[2][2] * image[i+1][j+1])
gy = (sobely[0][0] * image[i-1][j-1]) + (sobely[0][1] * image[i-1][j]) + \
(sobely[0][2] * image[i-1][j+1]) + (sobely[1][0] * image[i][j-1]) + \
(sobely[1][1] * image[i][j]) + (sobely[1][2] * image[i][j+1]) + \
(sobely[2][0] * image[i+1][j-1]) + (sobely[2][1] * image[i+1][j]) + \
(sobely[2][2] * image[i+1][j+1])
sobelxImage[i-1][j-1] = gx
sobelyImage[i-1][j-1] = gy
g = np.sqrt(gx * gx + gy * gy)
sobelGrad[i-1][j-1] = g
sobelxyImage = np.multiply(sobelxImage, sobelyImage)
# 3 Apply gaussian filter along x y and xy
size=3;# define the filter size
sigma=1; # define the standard deviation
size = int(size) // 2
xx, yy = np.mgrid[-size:size+1, -size:size+1]
normal = 1 / (2.0 * np.pi * sigma**2)
gg = np.exp(-((xx**2 + yy**2) / (2.0*sigma**2))) * normal
gaussx =gg
gaussy =gg
gaussxImage = np.zeros((M,N))
gaussyImage = np.zeros((M,N))
gaussxyImage = np.zeros((M,N))
gaussresult = np.zeros((M,N))
gaussimagex = np.pad(sobelxImage, (1,1), 'edge')
gaussimagey = np.pad(sobelyImage, (1,1), 'edge')
gaussimagexy = np.pad(sobelxyImage, (1,1), 'edge')
for i in range(1, M-1):
for j in range(1, N-1):
ggx = (gaussx[0][0] * gaussimagex[i-1][j-1]) + (gaussx[0][1] *gaussimagex[i-1][j]) + \
(gaussx[0][2] * gaussimagex[i-1][j+1]) + (gaussx[1][0] * gaussimagex[i][j-1]) + \
(gaussx[1][1] * gaussimagex[i][j]) + (gaussx[1][2] * gaussimagex[i][j+1]) + \
(gaussx[2][0] * gaussimagex[i+1][j-1]) + (gaussx[2][1] * gaussimagex[i+1][j]) + \
(gaussx[2][2] * gaussimagex[i+1][j+1])
ggy = (gaussy[0][0] * gaussimagey[i-1][j-1]) + (gaussy[0][1] * gaussimagey[i-1][j]) + \
(gaussy[0][2] * gaussimagey[i-1][j+1]) + (gaussy[1][0] * gaussimagey[i][j-1]) + \
(gaussy[1][1] * gaussimagey[i][j]) + (gaussy[1][2] * gaussimagey[i][j+1]) + \
(gaussy[2][0] * gaussimagey[i+1][j-1]) + (gaussy[2][1] * gaussimagey[i+1][j]) + \
(gaussy[2][2] * gaussimagey[i+1][j+1])
crossgg = (gg[0][0] * gaussimagexy[i-1][j-1]) + (gg[0][1] * gaussimagexy[i-1][j]) + \
(gg[0][2] * gaussimagexy[i-1][j+1]) + (gg[1][0] * gaussimagexy[i][j-1]) + \
(gg[1][1] * gaussimagexy[i][j]) + (gg[1][2] * gaussimagexy[i][j+1]) + \
(gg[2][0] * gaussimagexy[i+1][j-1]) + (gg[2][1] * gaussimagexy[i+1][j]) + \
(gg[2][2] * gaussimagexy[i+1][j+1])
gaussxImage[i-1][j-1] = ggx
gaussyImage[i-1][j-1] = ggy
gaussxyImage[i-1][j-1] = crossgg
blur = np.sqrt(ggx * ggx + ggy * ggy)
gaussresult[i-1][j-1] = blur
det = gaussxImage *gaussyImage - (gaussxyImage ** 2)
alpha = 0.04
trace = alpha * (gaussxImage +gaussyImage) ** 2
#finding the harris response
R = det - trace
# applying threshold
for i in range(1, M-1):
for j in range(1, N-1):
if R[i][j] > 140:
R[i][j]==0
else:
R[i][j]==255
f, ax1 = plt.subplots(1, 1, figsize=(5,5))
ax1.set_title('corners')
ax1.imshow(R, cmap="gray")
首先，您的代码中有几个错误：
最终阈值循环中的 R[i][j]==0 部分应该是 R[i][j]=0。请注意，您不必经过循环，只需执行 R[R>thresh]=255 等操作即可
如果我没记错的话，哈里斯算法中对应角点的R值是较大的正值。当我 运行 你的代码时，我得到 R 边角的负值，所以我怀疑那里有错误。
在这一点上，我不认为你代码中的主要问题是非极大值抑制，但如果它仍然是，这里是非极大值抑制的快速解释和我们在评论：
基本上，非最大抑制的想法很简单：如果一个点x的（角）响应不是邻域中最高的（你可以根据需要自由定义)，那么你就不要保留它。在您的情况下，将检测到的每个兴趣点的响应与其最近邻居的响应进行比较并仅在它们相对于所有兴趣点都更高时才保留它们可能就足够了。至于我们讨论的论文，其目的是以一种导致更均匀空间分布的方式抑制关键点（不是局部最大值）。让 S 成为关键点，按角响应的降序排列。这个想法是将它们中的每一个分配给一个“抑制半径”，即这些点不会被视为局部最大值的半径。由于S[0]是图像中角点响应最高的，永远不会被抑制，所以可以设置它的抑制半径radius_list[0]=inf。接下来我们来看S[1]。随着列表 S 的排序，唯一响应高于 S[1] 的点是 S[0]，由此得出 S[1] 停止为局部的半径最大值为 Dist(S[1], S[2])。也就是说，一旦我们将 S[0] 包含在 S[1] 的局部邻域中，由于 response[S[0]]>response[S[1]]，S[0] 将成为该邻域中的最大值。请注意，当您继续这样做时，您考虑的半径将变得越来越小。计算 radius_list 后，假设您需要 N 个特征点，您将只 select 具有最高 radius_list 值的 N 个点。在伪代码中：
#let S be the keypoints, sorted in decreasing corner response order
#Assume you want only to keep N keypoints at the end
radius=zeros(len(S))
radius[0]=inf
for i in range(len(S[1:])):
candidate_radii=[]
for j in range(0,i):
if response[i]<response[j]*some_const:#you can set some_const to something in [0.9,1]
candidate_radii.append(image_space_dist(S[i],S[j]))
radius[i]=min(candidate_radii)
sorted_indexes = argsort(radius)
kept_points = S[sorted_indexes][:N]
希望这对您有所帮助。
我正在尝试从头开始实施 Harris 角点检测算法。该算法的输出应该得到一个代表角的像素，但在我的代码中我得到了代表角的多个像素。这可能是因为未实现算法的最后部分 non-maximum suppression。这是我无法实施的，因为我没有正确理解它。如何实施？除此之外，我还试图找到这些角的坐标，如何在不使用 cv2 库的情况下做到这一点？
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.image as im
# 1. Before doing any operations convert the image into gray scale image
img = im.imread('OD6.jpg')
plt.imshow(img)
plt.show()
# split
R=img[:,:,0]
G=img[:,:,1]
B=img[:,:,2]
M,N=R.shape
gray_img=np.zeros((M,N), dtype=int);
for i in range(M):
for j in range(N):
gray_img[i, j]=(R[i, j]*0.2989)+(G[i, j]*0.5870)+(B[i, j]*0.114);
plt.imshow(gray_img, cmap='gray')
plt.show()
# 2. Applying sobel filter to get the gradients of the images in x and y directions
sobelx = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype = np.float)
sobely = np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype = np.float)
sobelxImage = np.zeros((M,N))
sobelyImage = np.zeros((M,N))
sobelGrad = np.zeros((M,N))
image = np.pad(gray_img, (1,1), 'edge')
for i in range(1, M-1):
for j in range(1, N-1):
gx = (sobelx[0][0] * image[i-1][j-1]) + (sobelx[0][1] * image[i-1][j]) + \
(sobelx[0][2] * image[i-1][j+1]) + (sobelx[1][0] * image[i][j-1]) + \
(sobelx[1][1] * image[i][j]) + (sobelx[1][2] * image[i][j+1]) + \
(sobelx[2][0] * image[i+1][j-1]) + (sobelx[2][1] * image[i+1][j]) + \
(sobelx[2][2] * image[i+1][j+1])
gy = (sobely[0][0] * image[i-1][j-1]) + (sobely[0][1] * image[i-1][j]) + \
(sobely[0][2] * image[i-1][j+1]) + (sobely[1][0] * image[i][j-1]) + \
(sobely[1][1] * image[i][j]) + (sobely[1][2] * image[i][j+1]) + \
(sobely[2][0] * image[i+1][j-1]) + (sobely[2][1] * image[i+1][j]) + \
(sobely[2][2] * image[i+1][j+1])
sobelxImage[i-1][j-1] = gx
sobelyImage[i-1][j-1] = gy
g = np.sqrt(gx * gx + gy * gy)
sobelGrad[i-1][j-1] = g
sobelxyImage = np.multiply(sobelxImage, sobelyImage)
# 3 Apply gaussian filter along x y and xy
size=3;# define the filter size
sigma=1; # define the standard deviation
size = int(size) // 2
xx, yy = np.mgrid[-size:size+1, -size:size+1]
normal = 1 / (2.0 * np.pi * sigma**2)
gg = np.exp(-((xx**2 + yy**2) / (2.0*sigma**2))) * normal
gaussx =gg
gaussy =gg
gaussxImage = np.zeros((M,N))
gaussyImage = np.zeros((M,N))
gaussxyImage = np.zeros((M,N))
gaussresult = np.zeros((M,N))
gaussimagex = np.pad(sobelxImage, (1,1), 'edge')
gaussimagey = np.pad(sobelyImage, (1,1), 'edge')
gaussimagexy = np.pad(sobelxyImage, (1,1), 'edge')
for i in range(1, M-1):
for j in range(1, N-1):
ggx = (gaussx[0][0] * gaussimagex[i-1][j-1]) + (gaussx[0][1] *gaussimagex[i-1][j]) + \
(gaussx[0][2] * gaussimagex[i-1][j+1]) + (gaussx[1][0] * gaussimagex[i][j-1]) + \
(gaussx[1][1] * gaussimagex[i][j]) + (gaussx[1][2] * gaussimagex[i][j+1]) + \
(gaussx[2][0] * gaussimagex[i+1][j-1]) + (gaussx[2][1] * gaussimagex[i+1][j]) + \
(gaussx[2][2] * gaussimagex[i+1][j+1])
ggy = (gaussy[0][0] * gaussimagey[i-1][j-1]) + (gaussy[0][1] * gaussimagey[i-1][j]) + \
(gaussy[0][2] * gaussimagey[i-1][j+1]) + (gaussy[1][0] * gaussimagey[i][j-1]) + \
(gaussy[1][1] * gaussimagey[i][j]) + (gaussy[1][2] * gaussimagey[i][j+1]) + \
(gaussy[2][0] * gaussimagey[i+1][j-1]) + (gaussy[2][1] * gaussimagey[i+1][j]) + \
(gaussy[2][2] * gaussimagey[i+1][j+1])
crossgg = (gg[0][0] * gaussimagexy[i-1][j-1]) + (gg[0][1] * gaussimagexy[i-1][j]) + \
(gg[0][2] * gaussimagexy[i-1][j+1]) + (gg[1][0] * gaussimagexy[i][j-1]) + \
(gg[1][1] * gaussimagexy[i][j]) + (gg[1][2] * gaussimagexy[i][j+1]) + \
(gg[2][0] * gaussimagexy[i+1][j-1]) + (gg[2][1] * gaussimagexy[i+1][j]) + \
(gg[2][2] * gaussimagexy[i+1][j+1])
gaussxImage[i-1][j-1] = ggx
gaussyImage[i-1][j-1] = ggy
gaussxyImage[i-1][j-1] = crossgg
blur = np.sqrt(ggx * ggx + ggy * ggy)
gaussresult[i-1][j-1] = blur
det = gaussxImage *gaussyImage - (gaussxyImage ** 2)
alpha = 0.04
trace = alpha * (gaussxImage +gaussyImage) ** 2
#finding the harris response
R = det - trace
# applying threshold
for i in range(1, M-1):
for j in range(1, N-1):
if R[i][j] > 140:
R[i][j]==0
else:
R[i][j]==255
f, ax1 = plt.subplots(1, 1, figsize=(5,5))
ax1.set_title('corners')
ax1.imshow(R, cmap="gray")
首先，您的代码中有几个错误：
最终阈值循环中的
R[i][j]==0部分应该是R[i][j]=0。请注意，您不必经过循环，只需执行R[R>thresh]=255等操作即可如果我没记错的话，哈里斯算法中对应角点的
R值是较大的正值。当我 运行 你的代码时，我得到R边角的负值，所以我怀疑那里有错误。
在这一点上，我不认为你代码中的主要问题是非极大值抑制，但如果它仍然是，这里是非极大值抑制的快速解释和我们在评论：
基本上，非最大抑制的想法很简单：如果一个点x的（角）响应不是邻域中最高的（你可以根据需要自由定义)，那么你就不要保留它。在您的情况下，将检测到的每个兴趣点的响应与其最近邻居的响应进行比较并仅在它们相对于所有兴趣点都更高时才保留它们可能就足够了。至于我们讨论的论文，其目的是以一种导致更均匀空间分布的方式抑制关键点（不是局部最大值）。让 S 成为关键点，按角响应的降序排列。这个想法是将它们中的每一个分配给一个“抑制半径”，即这些点不会被视为局部最大值的半径。由于S[0]是图像中角点响应最高的，永远不会被抑制，所以可以设置它的抑制半径radius_list[0]=inf。接下来我们来看S[1]。随着列表 S 的排序，唯一响应高于 S[1] 的点是 S[0]，由此得出 S[1] 停止为局部的半径最大值为 Dist(S[1], S[2])。也就是说，一旦我们将 S[0] 包含在 S[1] 的局部邻域中，由于 response[S[0]]>response[S[1]]，S[0] 将成为该邻域中的最大值。请注意，当您继续这样做时，您考虑的半径将变得越来越小。计算 radius_list 后，假设您需要 N 个特征点，您将只 select 具有最高 radius_list 值的 N 个点。在伪代码中：
#let S be the keypoints, sorted in decreasing corner response order
#Assume you want only to keep N keypoints at the end
radius=zeros(len(S))
radius[0]=inf
for i in range(len(S[1:])):
candidate_radii=[]
for j in range(0,i):
if response[i]<response[j]*some_const:#you can set some_const to something in [0.9,1]
candidate_radii.append(image_space_dist(S[i],S[j]))
radius[i]=min(candidate_radii)
sorted_indexes = argsort(radius)
kept_points = S[sorted_indexes][:N]
希望这对您有所帮助。