Я пытаюсь реализовать алгоритм обнаружения углов Харриса с нуля, но не смог реализовать не максимальное подавление

Я пытаюсь реализовать алгоритм обнаружения углов Харриса с нуля. Предполагается, что на выходе этого алгоритма должен быть один пиксель, представляющий угол, но в моем коде я получаю несколько пикселей, представляющих угол. Это может быть из-за того, что не реализована заключительная часть алгоритма 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")