I am trying to implement the Harris corner detection algorithm from scratch, but could not able to implement the non-maximum suppression

Question

I am trying to implement Harris corner detection algorithm from the scratch. The output of this algorithm should supposed to get one single pixel representing the corner, but in my code I am getting multiple pixel representing the corner. This is may be because of not implemented the final part of the algorithm that is non-maximum suppression. This I could not able to implement because I did not understand it properly. How to implement this? Along with this I am also trying to find the coordinates of these corner, how to do this with out using cv2 library?

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")

Answer 1

First of all, there are a couple of bugs in your code:

• The R[i][j]==0 part in the final thresholding loop should be R[i][j]=0. Note thought that you don't have to go through a loop, you can just do something like R[R>thresh]=255 etc.

• If I'm not mistaken, the R values that corresponds to corners in Harris' algorithm are the large positive ones. When I run your code, I get R values that are negative for edges and corners, so I suspect that there is a bug somewhere there.

At this point, I don't think that the main issue in your code is non-maxima suppression, but in case it still is, here is a quick explanation of non maxima suppression and the paper that we discussed in the comments:

Basically, the idea of non-maximal suppression is very simple: if the (corner) response of a point x is not the highest in a neighborhood (that you are free to define depending on your needs), then you don't keep it. In your case, it will probably be simply sufficient to compare the response of each of your detected interest points with the response of its closest neighbors and keep them only if they are higher with respect to all of them. As for the paper that we discussed, its aim is to suppress keypoints (that are not local maxima) in a way that results in a more uniform spatial distribution. Let S be the keypoints lits, sorted in decreasing order of corner response. The idea is to assign each of them to a "suppression radius", that is, a radius in which those points wont be considered a local maximum. As S[0] has the highest corner response in the image, it will never be suppressed, so you can set its radius of suppression radius_list[0]=inf. Next, let's look at S[1]. As the list S is sorted, the only point with highest response than S[1] is S[0], and from that, it follows that the radius at which S[1] stops being a local maximum is Dist(S[1], S[2]). That is, once we include S[0] in the local neighborhood of S[1], since response[S[0]]>response[S[1]], S[0] will become the maximum in that neighborhood. Note that as you continue like this, the radii that you consider will become smaller and smaller. Once you have computed radius_list, assuming you need N feature points, you will just select the N points that have the highest radius_list values. In pseudo-code:

#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]
    

Hope this helps.