Harris-Laplace角点检测初探+python代码实现

数学学院外行入门，课程团队作业，刚学的Markdown，版排得很随意（丑），多多担待

预备知识

角点

数字图像特征之一,包含了大量的位置信息以及指向型信息
一般而言，就是物体边缘的一些拐点，如下图所示：
上图中红色点为角点
用途：图像匹配、后期处理抵消图像形变、三维重建

Harris算法检测角点（重头戏来了）

Harris核心思路（通过周边区域来判断角点）

一个窗口在一检测点处上下、左右或者沿对角线移动：

平坦点flat：无论怎么移动，方框内的灰度值不变，则平平无奇和飞机场一样，没有什么点；
边缘点edge：方框内的灰度值只沿着一个方向改变，其余方向不变，则是边缘点；
角点corner：无论怎么移动，上下左右or斜向，灰度值变化都比较大，雨露均沾，则是角点。

Moravec算子

其实Harris是在Moravec基础上改进而来的，所以在讨论Harris之前，先来看一下Moravec

缺陷：Moravec只考虑了纵横斜着四种方向，所以没有旋转不变性，当原图旋转一下可能原来检测是角点的点现在检测就不是了，算法对图像的旋转鲁棒性差。

改进：那就考虑多一点方向呗，方向细化 =》全微分 =》泰勒展开
（这个想法太妙了，摆脱了对平移量[u,v]的依赖，仅依靠自身信息来判断角点,而且还增强了旋转不变性）

Harris算子

二次型，这里是不是有点椭圆曲线的意思，详见下图：

其实从线性代数的角度很好理解，就是要矩阵M的两个特征值都很大，为了表示“两个特征值都很大”这一层意思，给出下列角点响应公式，角点响应R大的话，则表明两个特征值都很大，则该点判断为角点。从数学的角度来看，这个角点响应取得真好，先是乘积打头，让两个值都尽可能地大，同时减去矩阵迹（trace）的平方，排除出现一个远大于另一个的情况（这种情况是边缘点），太妙了~

注：这里的M矩阵省略了权重w（x,y）
前面关于Harris角点检测说了这么多，所以。。。Laplace在哪，前面不都是一阶梯度吗？？？我是不是找错资料了？？？？

Harris-Laplace角点检测（原来主角这才来啊，哭了）

前言

传统的Harris角点检测（由Harris等在1988年提出）虽然对旋转改变有很好的检测不变性，但是却不具有尺度不变性及仿射不变性。为了获得尺度不变性，就需要在传统的Harris角点检测中引入多尺度空间。Mikolajczyk和Schmid于2004年提出具有尺度不变性的Harris-Laplace检测方法[1]。

算法流程

代码实现

以下代码虽然能够正常运行，但是很多细节我也还没搞明白，看看大致框架就行，想深究细节的话看最后参考文献[1].

import cv2
import math
import numpy as np
from scipy import ndimage as ndim
import matplotlib.pyplot as plt
%matplotlib inline# 取最近的正奇数
def odd(val):if int(val)%2 == 0:return int(val)+1else:return int(val)# 生成归一化后的高斯内核
def gauss_kernel(size,sigma):'''size表示高斯内核的大小，为正的奇数sigma表述随机变量X和Y的标准差，此处假设X和Y具有相同的标准差'''xy = range(-(size-1)//2,(size+1)//2)X,Y = np.meshgrid(xy,xy)G = 1/(2*np.pi*sigma**2)*np.exp(-0.5*(X**2+Y**2)/sigma**2)return G/np.sum(G) # 生成归一化后的LOG内核
def log_kernel(size,sigma):'''size表示高斯内核的大小，为正的奇数sigma表述随机变量X和Y的标准差，此处假设X和Y具有相同的标准差'''xy = range(-(size-1)//2,(size+1)//2)X,Y = np.meshgrid(xy,xy)LOG = (X**2+Y**2-2*sigma**2)/(sigma**4)*np.exp(-0.5*(X**2+Y**2)/sigma**2)return LOG# 非极大值抑制
def find_local_maximum(val):'''功能：查找领域极大值；输入：val-矩阵输出：row-领域极大值的行坐标       col-领域极大值的列坐标     max_local-领域极大值'''footprint = np.ones((3,3))footprint[1,1] = 0maxima_around = ndim.maximum_filter(val,footprint = footprint)  # 8领域内最大值滤波mask = val > maxima_around   # 找出极大值点所在的位置points = np.argwhere(mask==1)return points # plt.imshow(img,cmap = 'gray')
# plt.axis('off')
# plt.showdef Har_Lap_corners_detect(img): '''功能：Harris Laplace角点检测参数：img-输入的图像输出：points = [l,c,radius] l,c-角点的横纵坐标，radius-特征尺度值所对应的圆形区域的半径'''# 图像参数height,width = img.shape # 尺度参数sigma_begin = 1.5sigma_step = 1.2sigma_nb = 13sigma_array = (sigma_step**np.arange(sigma_nb))*sigma_begin# 第一部分：提取Harris角点harris_pts = np.zeros((0,3),dtype = float)k = 0.06 for i in range(sigma_nb):# 尺度s_I = sigma_array[i]ksize = odd(6*s_I+1)s_D = 0.7*s_I#     # 微分掩模#     x = np.arange(-round(3*s_D),round(3*s_D)+1)#     Ix = np.array([x*np.exp(-x*x/(2*s_D*s_D))/(s_D*s_D*s_D*(2*np.pi)**0.5)])#     Iy = Ix.T#     # 分别计算x，y方向上的梯度#     dx = cv2.filter2D(img,cv2.CV_32F,Ix)#     dy = cv2.filter2D(img,cv2.CV_32F,Iy)#     dxx = np.array([-1,0,1])g = np.array([np.exp(-x*x/(2*s_D**2))/(s_D*(2*np.pi)**0.5)])g = g/np.sum(g)   # 归一化fx = np.array([[-1,0,1]])dx = cv2.filter2D(img,cv2.CV_32F,g*fx)dy = cv2.filter2D(img,cv2.CV_32F,g.T*fx.T)# 自相关矩阵g = gauss_kernel(max(1,ksize),s_I) # 生成高斯核dx2 = cv2.filter2D(dx**2,cv2.CV_32F,g)dy2 = cv2.filter2D(dy**2,cv2.CV_32F,g)dxy = cv2.filter2D(dx*dy,cv2.CV_32F,g)cim = dx2*dy2-dxy**2-k*(dx2+dy2)**2# 阈值tmp = 0.01*np.max(cim)cim[cim<tmp] = 0# 3*3领域非极大值抑制points = find_local_maximum(cim)n = points.shape[0]points = np.c_[points,i*np.ones((n,1))]harris_pts =  np.r_[harris_pts,points]# 第二部分：Laplace变换# 计算尺度空间log = np.zeros((height,width,sigma_nb))for i in range(sigma_nb):s_L = sigma_array[i]#img_g = cv2.filter2D(img,-1,gauss_kernel(max(1,int(6*s_I+1)),s_L))#log[:,:,i] = s_L*s_L*cv2.filter2D(img,-1,log_kernel(odd(6*s_L+1),s_L))log[:,:,i] = s_L*s_L*cv2.Laplacian(cv2.GaussianBlur(img,(odd(6*s_L+1),odd(6*s_L+1)),s_L,s_L),cv2.CV_32F,ksize = 3)# 检测每个特征点在某一尺度LoG相应是否最大n = harris_pts.shape[0]cpt = -1points_t = np.zeros((n,3))for i in range(n):l = int(harris_pts[i,0])c = int(harris_pts[i,1])s = int(harris_pts[i,2])val = log[l,c,s]if s>0 and s<sigma_nb-1:if val>log[l,c,s-1] and val>log[l,c,s+1]:cpt += 1;points_t[cpt,:] = harris_pts[i,:]points_t[cpt,2] = 3*sigma_array[s]elif s == 0:if val > log[l,c,1]:cpt += 1;points_t[cpt,:] = harris_pts[i,:]points_t[cpt,2] = 3*sigma_array[s]else:if val > log[l,c,s-1]:cpt += 1;points_t[cpt,:] = harris_pts[i,:]points_t[cpt,2] = 3*sigma_array[s]points = np.zeros((cpt,3))for i in range(cpt):points[i,:] = points_t[i,:]return pointsdef draw(img,points):'''功能：在原图上标注出角点并显示输入：img-图像   points = [l,c,radius] l,c-角点的横纵坐标，radius-特征尺度值所对应的圆形区域的半径'''for i in range(points.shape[0]):cv2.circle(img,(int(points[i,1]),int(points[i,0])),int(points[i,2]),(255,0,0),2)  #  原来这里要横纵坐标互换，害我改了半天才发现错误所在，哭了。plt.figure(figsize=(7,7))plt.imshow(img)plt.axis('off')plt.show

Harris-Laplace角点检测

使用上述代码对检测图像角点，结果如下：

img = cv2.imread('wujiaoxing.jpg')
img1 = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
img2 = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)draw(img2,Har_Lap_corners_detect(img1))

Perfect！！！这次大作业有底了，逃出生天，哈哈哈哈哈~
下面换一个难一点的图试试看

img = cv2.imread('door.jpg')
img1 = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY)
img2 = cv2.cvtColor(img,cv2.COLOR_BGR2RGB)draw(img2,Har_Lap_corners_detect(img1))

。。。。这画的什么鬼

此处是编程过程中的小发现，与主题无关可不看
* 破案了，其实在MATLAB中卷积计算conv2（A，G，‘same’）与图像滤波fitler2(G,A)的计算差不多，就是相差了一个G矩阵旋转180度。但是对于高斯核G而言，其旋转180度不变，所以在这种情况下，此两者等价。
* 矩阵的卷积运算和我想的差不多，矩阵对应元素的相乘之后求和。
* matlab中的filter2计算和CV2中的cv2.filter2D 差不多，除了CV2中会把小于0的点直接设置为0

Harris-Laplace和Harris角点检测算法（胜利就在眼前了，加油！！！）

设计实验，对图像进行缩放和旋转，比较Harris-Laplace和Harris角点检测算法的稳定性

实验数据

实验数据来自图像数据集网站，Affine Covariant Features，选择6张存在旋转和放缩的图片，如下图红框所示：

实验结果

分别采用Harris以及Harris-Laplace算法对这六张图片进行角点检测，统计出角点的数量，如下表所示：

可见，Harris-Laplace算法检测出的角点数量要比Harris算法少得多。

可以用重复率来衡量角点检测算法的稳定性，重复率定义如下：r=NbpreNbinir = \frac{{{N_{bpre}}}}{{{N_{bini}}}}r=NbiniNbpre
式中，NbpreN_{bpre}Nbpre表示图像变化后检测出的角点，NbiniN_{bini}Nbini是的原始图像检测出的角点，本实验中原始图像选择为img1.重复率表示图像变化后检测出来的角点与未发生变化时检测出来的角点之比，因此重复率越高，表示角点检测算子越稳定。实验结果如下图所示：

由上图可见，在同样的旋转和缩放条件下，Harris-Laplace检测算子的稳定性要高于Harris角点检测算子。

总结

角点一般可以认为是物体边缘的一些拐点，属于图像特征信息，在很多场景下有用
Harris-Laplace算子是在传统Harris算子基础上引入尺度空间得到的，增强了传统算子的尺度不变性
通过实验发现，Harris-Laplace算子所检测到的角点数目比传统Harris算子少得多，而且在同样的旋转和缩放条件下更稳定。
代码使人头大，这也太顶了吧！！！不管了，我肝不动了，就这样吧。
经验教训：这一次把大量的时间花在了代码的实现上，导致之后的时间很少，实验做得很水，很多场景以及参数的比较都没有做，下一次代码可能就没有这么细致了，能调库就调库，给前期的算法赏析（“赏析”这个词用得真好）以及后期的实验环节多一点时间。

参考文献

K. Mikolajczyk and C. Schmid. Scale & affine invariant interest point detectors. International Journal of Computer Vision, 2004
matlab实现：《现代数字图像处理技术提高及应用案例详解（MATLAB版）》赵小川
人工智能：harris角点检测

Harris-Laplace角点检测初探+python实现相关推荐

Harris的角点检测和特征匹配
一.特征检测(提取) 基于特征的图像配准方法是图像配准中最常见的方法之一.它不是直接利用图像像素值,二十通过像素值导出的符号特征(如特征点.特征线.特征区域)来实现图像配准,因此可以克服利用灰度信息进 ...
python 角点检测_opencv python Shi-Tomasi角点检测
cv2.goodFeaturesToTrack(image, maxCorners, qualityLevel, minDistance[, corners[, mask[, blockSize[, ...
SURF角点检测（python）
简述算是SIFT算法的改进算法. 原图算法图代码 import cv2imgpath = '5.jpg' img = cv2.imread(imgpath) gray = cv2.cvtColo ...
python+OpenCV笔记（二十四）：Shi-Tomasi角点检测
Shi-Tomasi角点检测原理 python+OpenCV笔记(二十二):角点检测原理(Harris角点检测原理.Shi-Tomasi角点检测原理)https://blog.csdn.net/qq ...
harris角点检测与ncc匹配
转自:http://zixuanjinan.blog.163.com/blog/static/11543032620097510122831/ file1:---------------------- ...
测角误差估计算法matlab,Harris角点检测及 Matlab实验
目录 1 基础知识 1.1 图像变化的类型 1.2 提取点特征的作用 1.3 什么是好的角点检测算法? 2 Harris 角点检测 2.1 Harris角点检测基本思想 2.2 Harris角点检测: ...
harris角点检测_角点检测(2) - harris算子 - 理论与Python代码
数字图像,图像=矩阵,[m*n]从[0,255]的灰度值角点检测:物体边缘的拐点 ->应用:图像匹配与检索.图像物体形变恢复(摄像机标定).三维重建 Harris角点检测(早期,原理简单,视频 ...
cv2.cornerHarris()详解　python+OpenCV 中的 Harris 角点检测
原文作者:aircraft 原文地址:https://www.cnblogs.com/DOMLX/p/8763369.html 参考文献----------OpenCV-Python-Toturial ...
Python+OpenCV：图像Harris角点检测(Harris Corner Detection)
Python+OpenCV:图像Harris角点检测(Harris Corner Detection) 理论 corners are regions in the image with large v ...

Harris-Laplace角点检测初探+python实现