电脑基础 · 2023年3月31日

Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)

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、基本实现步骤

(1)相似判断

Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)(式1)
I(r0)为中心位置像素,I( r)为周边位置像素,t像素偏差阈值,c为布尔结果,表示相似或不相似。
关于阈值t:可见,t越大,周边像素越容易与中心位置像素“相似”,越不容易被判定为边缘或角点,即将获得较少的角点。相关文献说,t值影响的是角点检测的数量[1]。

(2)平滑曲线进行相似判断

实际中,使用平滑曲线来求c
Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)(式2)

(3)累计相似度

对于圆形模板区域,共计37个像素,周边36个像素与中心位置像素进行对比,可以得到该中心位置的36个c。对这36个c进行累加,得到该中心位置的累计相似度:
Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制) (式3)
显然,n越大,越不可能是边缘、角点。

(4)初始边缘响应/角点

这里,再引入一个阈值g,通过判断n与g的大小,来得到该中心位置属于边缘或角点的可能性。
Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)(式4)
可见,若n大于g,则R为0,表示该中心位置为非边缘(非角点);如果n小于g,R取值为(g-n),含有取负操作。故R值越大,该中心位置为边缘或角点的可能性就越大。

关于阈值g:g的取值越小,对角点检测越挑剔。文献[1]阐述,g影响角点检测的质量。g的取值,文献[1]阐述可以取max(n)的0.75倍;有的程序代码中,直接取值37/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()

Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)
Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)

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实现 (边缘检测、角点检测、重心计算、非极大值抑制)
与原始角点的局部对比
Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)

5、非极大值抑制

非极大值抑制,也可对角点进行剔除。
思想:对于已获取的角点矩阵,使用一3x3的矩形区域进行逐个判断。在3x3的区域中,若中心位置不是最大,则认为该中心位置不是局部极值,即认为该中心位置不是角点,将其抑制。
问题:若3x3区域内,有两个位置的值相等且均为最大,则如何处理?目前,本文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()

5.3 检测结果

Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)
与重心法处理后的局部对比
Susan角点检测python实现 (边缘检测、角点检测、重心计算、非极大值抑制)

6. 源码下载:

https://download.csdn.net/download/xiaohuolong1827/85030525

7. 其他

2022.03.22,黄宁然:“SUSAN,你竟然在看这个,我以前用过,你信不信”。呃,我当然信。