Susan角点检测（边缘检测、角点检测、重心计算、非极大值抑制）

写在前面

参考文献镇楼

[1]https://blog.csdn.net/tostq/article/details/49305615
[2]https://blog.csdn.net/qq_45613931/article/details/117819107
[3]陈丽莉《基于SUSAN算法的角点检测》
[4]王栋、朱明《SUSAN角点探测算法分析改进》

问题来源

an***** xue100: https://bbs.csdn.net/topics/*********?spm=1001.2014.3001.**77
1）如何计算USAN区域的重心？
2）非极大值抑制，指的是若模板中心对应的灰度小于模板内的任意像素灰度，就舍弃该像素点吗？
3）几何门限g如何选取？
因“当前发帖距今超过3年，不再开放新的回复”，故新建帖子。迟到的回复。

1、原理简介

SUSAN（Smallest UnivalueSegment Assimilating Nucleus）使用一个圆形模板，通过检测模板中的像素与中心位置像素的偏离程度，来判断中心位置像素是否为边缘或角点。
模板大小：半径为3.5像素，模板内共计有37个像素。
Susan角点检测python实现（边缘检测、角点检测、重心计算、非极大值抑制）
中心思想：如果周边像素与中心位置像素偏差较小，则认为周边像素与中心位置像素相似，中心位置为非边缘；如果周边像素与中心位置偏差较大，则认为中心位置为边缘或角点。

2、基本实现步骤

3、python代码实现

3.1 圆形掩模

def get_susan_mask():
   mask=np.ones((7,7))
   mask[0,0]=0;mask[0,1]=0;mask[0,5]=0;mask[0,6]=0;
   mask[1,0]=0;mask[1,6]=0;
   mask[5,0]=0;mask[5,6]=0;mask[6,0]=0;mask[6,1]=0;
   mask[6,5]=0;mask[6,6]=0;
   return mask

该程序返回的mask如下：
Susan角点检测python实现（边缘检测、角点检测、重心计算、非极大值抑制）
为1的位置，即是圆形模板区域。

3.2 susan角点检测

def susan_corner_detect(img_src,t=10):
   susan_mask = get_susan_mask()
   img = img_src.copy()
   row_s,col_s = 3,3
   row_e,col_e = img_src.shape[0]-3,img.shape[1]-3
   n_max = 0
   n_arr=37*np.ones(img.shape) # 初始认为没有角点
   for r in range(row_s,row_e):#遍历所有行
      for c in range(col_s,col_e):#遍历所有列
         susan_zone = img[r-3:r+3+1,c-3:c+3+1]#获取矩形区域
         susan_zone = susan_zone[susan_mask!=0]#使用mask截取圆形区域
         r0 = img[r,c]
         similarity = np.exp(-((1.0*susan_zone-r0)/t)**6 )
         n=np.sum(similarity)
         if n>n_max:
            n_max = n
         n_arr[r,c] = n
   g = n_max /2
   R = np.zeros(img.shape)
   index = n_arr<g       #小于g，认为是可能的角点，越小，可能性越大
   R[index] = g-n_arr[index] # 取反，所以R越大，是角点的可能性越大

   plt.figure()
   plt.title("edge")
   plt.imshow((6.37 * n_arr).astype(np.uint8), cmap=cm.gray)
   return R

3.3 检测结果

取t=10，g= n_max /2。
主程序中，调用算法：

if __name__ == '__main__':
    img_src = cv2.imread('susan_input1.png',-1)
    if len(img_src.shape)==3:
       img_src = cv2.cvtColor(img_src,cv2.COLOR_BGR2GRAY)
    corner = susan_corner_detect(img_src)

    img_show = cv2.cvtColor(img_src, cv2.COLOR_GRAY2BGR)
    img_show[corner != 0] = (255, 0, 0)
    plt.figure()
    plt.title("original corners")
    plt.imshow(img_show, cmap=cm.gray)
    plt.show()

4、重心法去除伪角点

从上图的角点检测结果来看，原始角点中存在较多冗余角点，或伪角点。
采用重心法去除。
根据文献1，求取重心
Susan角点检测python实现（边缘检测、角点检测、重心计算、非极大值抑制）（式5）
上式中，r为位置，实际为（x，y）两轴构成。计算时，先计算x方向重心，再计算y方向重心。对公式的简言之：圆形模板中心位置周边有36个像素，相应的有36个c，在求x方向重心Gx时，将36个位置的x坐标与相应的c相乘后累加，再除以sum( c)，即可得到Gx；同理得Gy。另外，对于r0，可假设其位置为(0,0)。
在求得重心后，计算重心与中心位置的距离，根据文献[3]描述，如距离较小，则认为中心位置不是角点。根据文献[4]，该距离阈值设置为1.8。本文采用1.5。另，重心法去除伪角点与初始角点的检测，在代码实现时，有些许内容的重复，所以在初始角点检测时，即可同步应用重心法。

4.1重心法代码

def gravity_filter(img_src,corner_src,t=10,F=1.5):
   x_label = np.zeros((7,7))
   y_label = np.zeros((7,7))
   x_label[:,0]=-3;x_label[:,1]=-2;x_label[:,2]=-1;
   x_label[:, -1] = 3;    x_label[:, -2] = 2;    x_label[:, -3] = 1;
   y_label[0,:]=-3;y_label[1,:]=-2;y_label[2,:]=-1;
   y_label[4, :] = 1;y_label[5, :] = 2;   y_label[6, :] = 3;
   print(x_label,"\r\n",y_label) #查看矩形区域内x、y轴信息

   img = img_src.copy()
   row_s, col_s = 3, 3
   row_e, col_e = img_src.shape[0] - 3, img.shape[1] - 3
   corner = corner_src.copy()
   susan_mask = get_susan_mask()
   for r in range(row_s,row_e):
      for c in range(col_s,col_e):
         if corner[r,c] ==0: #对于不是角点的位置，就没必要进行后续计算了
            continue
         susan_zone = img[r-3:r+3+1,c-3:c+3+1]#获取矩形区域
         r0 = img[r,c]
         similarity = np.exp(-((1.0*susan_zone-r0)/t)**6 )
         g_x = np.sum(similarity[susan_mask==1]*x_label[susan_mask==1] )/np.sum(similarity[susan_mask==1])#使用mask截取圆形区域
         g_y = np.sum(similarity[susan_mask == 1] * y_label[susan_mask == 1]) / np.sum(similarity[susan_mask == 1])#使用mask截取圆形区域
         distance = np.sqrt(g_x**2+g_y**2)
         if distance<F:
            corner[r,c] = 0
   return corner

4.2 主程序调用

if __name__ == '__main__':
    img_src = cv2.imread('susan_input1.png',-1)
    if len(img_src.shape)==3:
       img_src = cv2.cvtColor(img_src,cv2.COLOR_BGR2GRAY)
    corner = susan_corner_detect(img_src)

    img_show = cv2.cvtColor(img_src, cv2.COLOR_GRAY2BGR)
    img_show[corner != 0] = (255, 0, 0)
    plt.figure()
    plt.title("original corners")
    plt.imshow(img_show, cmap=cm.gray)


    cor_g = gravity_filter(img_src, corner)
    img_show2 = cv2.cvtColor(img_src, cv2.COLOR_GRAY2BGR)
    img_show2[cor_g != 0] = (255, 0, 0)
    plt.figure()
    plt.title("corners-gravity ")
    plt.imshow(img_show2, cmap=cm.gray)
	plt.show()

4.3 检测结果

Susan角点检测python实现（边缘检测、角点检测、重心计算、非极大值抑制）
与原始角点的局部对比

5、非极大值抑制

非极大值抑制，也可对角点进行剔除。
思想：对于已获取的角点矩阵，使用一3×3的矩形区域进行逐个判断。在3×3的区域中，若中心位置不是最大，则认为该中心位置不是局部极值，即认为该中心位置不是角点，将其抑制。
问题：若3×3区域内，有两个位置的值相等且均为最大，则如何处理？目前，本文2个均保留。
220719更新：该nms代码存在些许不足之处，详见后续的博文，https://blog.csdn.net/xiaohuolong1827/article/details/125859795

5.1 非极大值抑制代码

def corner_nms(corner,kernal=3):
   out = corner.copy()
   row_s = int(kernal/2)
   row_e = out.shape[0] - int(kernal/2)
   col_s,col_e = int(kernal/2),out.shape[1] - int(kernal/2)
   for r in range(row_s,row_e):
      for c in range(col_s,col_e):
         if corner[r,c]==0: #不是可能的角点
            continue
         zone = corner[r-int(kernal/2):r+int(kernal/2)+1,c-int(kernal/2):c+int(kernal/2)+1]
         index = corner[r,c]<zone
         (x,y) = np.where(index==True)
         if len(x)>0 : #说明corner[r,c]不是最大，直接归零将其抑制
            out[r,c] = 0
         else:
            out[r,c] = 255
   return out

5.2 主程序调用

if __name__ == '__main__':
    img_src = cv2.imread('susan_input1.png',-1)
    if len(img_src.shape)==3:
       img_src = cv2.cvtColor(img_src,cv2.COLOR_BGR2GRAY)
    corner = susan_corner_detect(img_src)

    img_show = cv2.cvtColor(img_src, cv2.COLOR_GRAY2BGR)
    img_show[corner != 0] = (255, 0, 0)
    plt.figure()
    plt.title("original corners")
    plt.imshow(img_show, cmap=cm.gray)


    cor_g = gravity_filter(img_src, corner)
    img_show2 = cv2.cvtColor(img_src, cv2.COLOR_GRAY2BGR)
    img_show2[cor_g != 0] = (255, 0, 0)
    plt.figure()
    plt.title("corners-gravity ")
    plt.imshow(img_show2, cmap=cm.gray)


    cor_g_nms = corner_nms(cor_g)
    img_show3 = cv2.cvtColor(img_src, cv2.COLOR_GRAY2BGR)
    img_show3[cor_g_nms != 0] = (255, 0, 0)
    plt.figure()
    plt.title("corners-gravity-nms ")
    plt.imshow(img_show3, cmap=cm.gray)
    plt.show()