susan算子python实现

文章是自己写的，仅供参考如有错误欢迎指正（河理社恐男孩李某注）

Susan算子代码思路：

SUSAN算子原理，SUSAN使用一个圆形模板，通过对比检测模板中的像素与中心像素的偏离程度，来判断中心像素是否为边缘或角点。

代码实现步骤：

1. 首先创建圆形模板，进行相似性判断，判断公式如下：

公式的大概意思就是将圆形模板的像素点像素和中心像素进行对比，c为布尔结果，如果他们之差小于返回值为1反之为0，

关于阈值t：t越大周边像素越容易被判定与中心像素相似所以返回的角点数量越少

2. 累计像素度，对于圆形模板区域共37个像素，周围36个像素与中心像素做比较，对这36个像素返回的布尔值c进行累加，得到该中心位置的累计相似度（这里记作n），累计相似度越大，这个点越不可能是边缘，角点。

3. 初始边缘响应：这里需要再引入一个阈值g，将某个像素的累计相似度与g做比较，若n>g，表示该点为非边缘点（非角点），如果n<g，n-g差值越大则改点为边缘或角点的可能性就越大。

4. 重心法进行伪角点去除，为了测量的准确性，我们需要进行伪角点去除，这里我采用的是重心法进行伪角点去除：

式中，r为像素点位置，我们在计算时先计算x方向重心，再计算y方向上的重心，简而言之：圆形模板周围有36个像素，对应的有36个c求x方向上的重心时将36个位置的x坐标与相应的c相乘后累加再初一sum（c）即可，y方向上的重心求法同理，在求得重心后，计算重心与中心位置的距离，如果距离较小，则认为重心位置不是角点

下面附上代码：

import numpy as np
import cv2
import time
start =time.perf_counter()
#首先读取图像并将其改为灰度图像
img=cv2.imread(r"C:\Users\lihuiliang\Desktop\test_pattern(1).tif",0)
image=cv2.imread(r"C:\Users\lihuiliang\Desktop\test_pattern(1).tif",1)
#创建一个和图像大小相等的二维矩阵#h为y,w为x
h,w=img.shape
Hwin=np.zeros((h,w),dtype=np.uint8)
#创建圆形模板
X=[-3,-3,-3,-2,-2,-2,-2,-2,
-1,-1,-1,-1,-1,-1,-1,0, 0, 0, 0, 0, 0, 0,1 ,1 ,1 ,1 ,1 ,1 , 1,2,2,2,2,2,3,3,3]
Y=[-1,0,-1,-2,-1,0,1,2,-3,-2,-1,0,1,2,3,-3,-2,-1,0,1,2,3,-3, -2, -1, 0, 1, 2, 3,-2,-2,0,1,2,-1,0,1]
#创建一个大小和原图片相同的二维数组
img_huidu=np.zeros((w,h))#将原图像的像素灰度值存入灰度矩阵
for i in range(w):for j in range(h):img_huidu[i][j]=img[j][i]
#门限阈值，几何阈值
t=24
g=18
xxx=[]
yyy=[]
#创建一个m矩阵
m=np.zeros((w,h))
#遍历所有像素点，注意范围
for x in range(3,w-3):for y in range(3,h-3):#中心坐标点为img[x,y],圆形点坐标为img[x+X[i],y+Y[i]]same = 0ggx = []ggy = []zzx = []zzy = []c=[]for i in range(37):xx=(x+X[i])yy=(Y[i]+y)if abs(img_huidu[(xx),(yy)])-int(img_huidu[x,y])<=t  :m[xx][yy]=1same=same+1if abs(img_huidu[(xx),(yy)])-int(img_huidu[x,y])>t :m[xx][yy] = 0if i==36:for j in range(37):#去除伪角点#分别求x，y方向上面的重心x_x=(x+X[j])y_y=(y+Y[j])#36个x坐标与c相乘然后相加gx=m[x_x][y_y]*x_xggx.append(gx)#求所有gx的和g_x=sum(ggx)#36个y坐标与c相乘后相加gy=m[x_x][y_y]*y_yggy.append(gy)#求所有gy的和g_y=sum(ggy)#计算所有c的和c.append(m[x_x][y_y])cc=sum(c)if (cc-m[x,y])==0:continueif j==36:Gx=(g_x-m[x][y]*xx)/(cc-m[x,y])Gy = (g_y - m[x][y]*yy) / (cc-m[x,y])juli=((Gx-x)**2+(Gy-y)**2)**0.5if juli <=1:same=g+1if(same<g):xxx.append(x)yyy.append(y)
print(m)
men=len(xxx)
print(men)
for i in range(men):image = cv2.drawMarker(image, (xxx[i], yyy[i]), (0,255,0),markerType=0)Hwin[yyy[i],xxx[i]]=255
end = time.perf_counter()
print('共运行时间为: %s 秒'%(end-start))
cv2.imshow('win',Hwin)
cv2.waitKey(0)
cv2.imwrite('susan.png', image, [cv2.IMWRITE_PNG_COMPRESSION, 0])cv2.imshow('img',image)
cv2.waitKey(0)