Harris角点检测和ShiTomasi角点检测

zoukankan html css js c++ java

Harris角点检测和ShiTomasi角点检测
https://www.cnblogs.com/zyly/p/9508131.html

目录
OpenCV可以检测图像的主要特征，然后提取这些特征、使其成为图像描述符，这类似于人的眼睛和大脑。这些图像特征可作为图像搜索的数据库。此外，人们可以利用这些关键点将图像拼接起来，组成一个更大的图像，比如将许多图像放在一块，然后形成一个360度全景图像。

这里我们将学习使用OpenCV来检测图像特征，并利用这些特征进行图像匹配和搜索。我们会选取一些图像，并通过单应性，检测这些图像是否在另一张图像中。

一特征检测算法

有许多用于特征检测和提取的算法，我们将会对其中大部分进行介绍。OpenCV最常使用的特征检测和提取算法有：
- Harris：该算法用于检测角点；
- SIFT：该算法用于检测斑点；
- SURF：该算法用于检测角点；
- FAST：该算法用于检测角点；
- BRIEF：该算法用于检测斑点；
- ORB：该算法代表带方向的FAST算法与具有旋转不变性的BRIEF算法；
通过以下方法进行特征匹配：
- 暴力(Brute-Force)匹配法；
- 基于FLANN匹配法；
可以采用单应性进行空间验证。

二特征定义

那么，究竟什么是特征呢？为什么一副图像的某个特定区域可以作为一个特征，而其他区域不能呢？粗略的讲，特征就是有意义的图像区域，该区域具有独特特征和易于识别性。因此角点及高密度区域都是很好的特征，而大量重复的模式或低密度区域(例如图像中的蓝色天空)则不是很好的特征。边缘可以将图像分为两个区域，因此也可以看做好的特征。斑点是与周围有很大差别的像素区域，也是有意义的特征。

大多数特征检测算法都会涉及图像的角点、边和斑点的识别，也有一些涉及脊向的概念，可以认为脊向是细长物体的对称轴，例如识别图像中的一条路。角点和边都好理解，那什么是斑点呢？斑点通常是指与周围有着颜色和灰度差别的区域。在实际地图中，往往存在着大量这样的斑点，如一颗树是一个斑点，一块草地是一个斑点，一栋房子也可以是一个斑点。由于斑点代表的是一个区域，相比单纯的角点，它的稳定性要好，抗噪声能力要强，所以它在图像配准上扮演了很重要的角色。

由于某些算法在识别和提取某类型特征的时候有较好的效果，所以知道输入图像是什么很重要，这样做有利于选择最合适的OpenCV工具包。

三 Harris检测角点特征

在这之前其实我们已经接触过角点检测了，在相机标定的时候，我们就利用到了角点检测。不过那时候没有深入的去研究。在这里，我们将会深入原理去学习角点检测。

下面我们从使用cornerHarris()函数讲起。
cornerHarris(src, blockSize, ksize,k[,dst[,borderType]]);
参数详解:
- image:输入的单通道8位或者浮点图像；
- blockSize：就是扫描时候窗口的大小。
- ksize：cornerHarris函数会使用Sobel算子，该参数定义了Sobel算子的中孔。简单来说，该函数定义了角点检测的敏感度，其值必须介于3~31之间的奇数。
- k：harris 计算响应公式中的
- borderType：像素插值方法；
函数 cornerHarris 对输入图像进行 Harris 边界检测。输出是一幅浮点值图像，大小与输入图像大小相同，浮点值越高，表明越可能是特征角点(我们可以对图像进行阈值化)。
# -*- coding: utf-8 -*- """ Created on Mon Aug 20 20:17:34 2018 @author: lenovo """ ''' Harris角点检测 ''' import cv2 import numpy as np img = cv2.imread('./image/cali.bmp') img = cv2.resize(img,dsize=(600,400)) gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) gray = np.float32(gray) #角点检测第三个参数为角点检测的敏感度，其值必须介于3~31之间的奇数 dst = cv2.cornerHarris(gray,3,23,0.04) print(dst.shape) #(400, 600) img[dst>0.01*dst.max()] = [0,0,255] cv2.imshow('',img) cv2.waitKey(0) cv2.destroyAllWindows()
运行结果如下：

如果我们把第三个参数改为3:，可以看到：

四 Harris检测原理

上面我们已经通过实例演示了Harris检测的效果，相信你对Harris角点检测已经有了初步的认识。这里我将带你深入了解Harris角点检测的原理。

我们先来看一幅图片，了解一下什么是角点？

上图中E,F中的角我们通常称作角点（corner points），他们具有以下特征：
- 轮廓之间的交点；
- 对于同一场景，即使视角发生变化，通常具备稳定性质的特征；
- 该点附近区域的像素点无论在梯度方向上还是其梯度幅值上有着较大变化；
harris特征角最早在paper A Combined Corner and Edge Detector中被Chris Harris & Mike Stephens提出。

Harris角点检测的基本思想：算法基本思想是使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，都有着较大灰度变化，那么我们可以认为该窗口中存在角点。

1、灰度变化描述

当窗口发生

$E (u, v) = \sum_{(x, y) € W} w (x, y) [I (x + u, y + v) - I (x, y)]^{2}$

参数解释：
我们的窗口函数通常有如下两种形式：

根据上述表达式，当窗口在平潭区域上移动，可以想象得到，灰度不会发生什么变换。

2、

首先需要了解泰勒公式，任何一个函数表达式，均可有泰勒公式进行展开，以逼近原函数，我们可以对下面函数进行一阶展开(如果对泰勒公式忘记了，可以翻翻本科所学的高等数学)。

$f (x + u, y + v) \approx f (x, y) + u f_{x} (x, y) + v f_{y} (x, y)$

那么

$\sum_{(x, y) € W} w (x, y) [I (x + u, y + v) - I (x, y)]^{2}$

$\approx \sum_{(x, y) € W} w (x, y) [I (x, y) + u I_{x} + v I_{y} - I (x, y)]^{2}$

$= \sum_{(x, y) € W} w (x, y) [u^{2} I_{x}^{2} + 2 u v I_{x} I_{y} + v^{2} I_{y}^{2}]$

$= \sum_{(x, y) € W} w (x, y) [\begin{matrix} u & v \end{matrix}] [\begin{matrix} I_{x}^{2} & I_{x} I y \\ I_{x} I_{y} & I_{y}^{2} \end{matrix}] [\begin{matrix} u \\ v \end{matrix}]$

$= [\begin{matrix} u & v \end{matrix}] (\sum_{(x, y) € W} w (x, y) [\begin{matrix} I_{x}^{2} & I_{x} I y \\ I_{x} I_{y} & I_{y}^{2} \end{matrix}]) [\begin{matrix} u \\ v \end{matrix}]$

所以

$E (u, v) = [\begin{matrix} u & v \end{matrix}] M [\begin{matrix} u \\ v \end{matrix}]$

其中：

3、矩阵

难道我们是直接求上述的

下图是对这三种情况窗口中的对应像素的梯度分布进行绘制：

不知道大家有没有注意到这三种区域的特点:
- 平坦区域上的每个像素点所对应的
- 边缘区域有一坐标轴分布较散，至于是哪一个坐标上的数据分布较散不能一概而论，这要视边缘在图像上的具体位置而定，如果边缘是水平或者垂直方向，那么
- 角点区域的
我们是不是可以根据这些特征来判断哪些区域存在角点呢？

上面我们已经计算出了

我们可以将

$E (u, v) = A u^{2} + 2 C u v + B v^{2}$

其中

$A = \sum_{(x, y) € W} w (x, y) * I_{x}^{2}$

$B = \sum_{(x, y) € W} w (x, y) * I_{y}^{2}$

$C = \sum_{(x, y) € W} w (x, y) * I_{x} I_{y}$

二次项函数本质上就是一个椭圆函数。椭圆的长和宽是由

$[\begin{matrix} u & v \end{matrix}] M [\begin{matrix} u \\ v \end{matrix}] = 1$

如果使用椭圆进行数据集表示，则绘制图示如下：

虽然我们利用

把

$\bar{I_{x}} = \sum_{i = 1}^{m} I_{x i}$

$\bar{I_{y}} = \sum_{i = 1}^{m} I_{y i}$

我们仍然使用

$X = [\begin{matrix} I_{x 1} & I_{x 2} & . . . & I_{x m} \\ I_{y 1} & I_{y 2} & . . . & I_{y m} \end{matrix}]$

则对应协方差矩阵可以写成：

$C = \frac{1}{m} X X^{T} = \frac{1}{m} \sum_{i = 1}^{m} [\begin{matrix} I_{x}^{2} & I_{x} I y \\ I_{x} I_{y} & I_{y}^{2} \end{matrix}]$

但矩阵

$M = \sum_{(x, y) € W} [\begin{matrix} I_{x}^{2} & I_{x} I y \\ I_{x} I_{y} & I_{y}^{2} \end{matrix}]$

如果我们对协方差矩阵
- 如果两个字段
- 如果是平坦区域，那么像素点的梯度所构成的点集比较集中在原点附近，因为窗口区域内的像素点的梯度幅值非常小，此时矩阵
- 如果是边缘区域，在计算像素点的
注：

因此可以得出下列结论：
- 特征值都比较大时，即窗口中含有角点；
- 特征值一个较大，一个较小，窗口中含有边缘；
- 特征值都比较小，窗口处在平坦区域；
4、如何度量角点响应

通常用下面表达式进行度量，对每一个窗口计算得到一个分数

$R = d e t (M) - k (t r a c e (M))^{2}$

$d e t (M) = λ_{1} λ_{2}$

$t r a c e (M) = λ_{1} + λ_{2}$

这里

但是为什么会使用这样的表达式呢？一下子是不是感觉很难理解？其实也不难理解，函数表达式一旦出来，我们就可以绘制它的图像，而这个函数图形正好满足上面几个区域的特征。通过绘制函数图像，直观上更能理解。绘制的

所以说难点不在于理解这个函数表达式，而在于如何创造出这个函数表达式。Harris也许对很多函数模型非常了解，对于创造出这样的一个函数表达式，易如反掌，当然在我们看来感觉是很了不起的，那是因为我们见过的函数模型太少。

最后设定R的阈值，进行角点判断。当然其中还有些后处理步骤就不再说了，比如说角点的极大值抑制等。

非极大值抑制的原理，在一个窗口内，如果有很多角点则用值最大的那个角点，其他的角点都删除。

注意：有些朋友在问是如何通过矩阵判断角点的？其实上面，我们已经推导出

五自己实现Harris角点检测

Harris角点检测可以分为5个步骤：

1、计算图像

$I_{x} = \frac{\partial I}{\partial x} = I \times [\begin{matrix} - 1 & 0 & 1 \end{matrix}]$

$I_{y} = \frac{\partial I}{\partial y} = I \times {[\begin{matrix} - 1 & 0 & 1 \end{matrix}]}^{T}$

2、计算图像两个方向梯度的乘积;

$I_{x}^{2} = I_{x} * I_{x}$

$I_{y}^{2} = I_{y} * I_{y}$

$I_{x} I_{y} = I_{x} * I_{y}$

3、使用高斯函数对

$A = \sum_{(x, y) € W} g (I_{x}^{2}) = \sum_{(x, y) € W} I_{x}^{2} * w (x, y)$

$B = \sum_{(x, y) € W} g (I_{y}^{2}) = \sum_{(x, y) € W} I_{y}^{2} * w (x, y)$

$C = \sum_{(x, y) € W} g (I_{x} I_{y}) = \sum_{(x, y) € W} I_{x} I_{y} * w (x, y)$

其中

4、计算每个像素点

$R = d e t (M) - k (t r a c e (M))^{2}$

5、过滤大于某一阈值

$R = {R : d e t (M) - k (t r a c e (M))^{2} > t}$

如果需要在3×3或者5×5的邻域进行非最大值抑制，则局部最大值点即为图像中的角点。

以上所有公式*表示逐像素点乘。
# -*- coding: utf-8 -*- """ Created on Mon Aug 20 20:17:34 2018 @author: lenovo """ ''' Harris角点检测 ''' import cv2 import numpy as np def test(): ''' 调用系统库函数进行测试 ''' img = cv2.imread('./image/cali.bmp') img = cv2.resize(img,dsize=(600,400)) gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) gray = np.float32(gray) #角点检测第三个参数为角点检测的敏感度，其值必须介于3~31之间的奇数 dst = cv2.cornerHarris(gray,3,3,0.04) print(dst.shape) #(400, 600) img[dst>0.01*dst.max()] = [0,0,255] cv2.imshow('',img) cv2.waitKey(0) cv2.destroyAllWindows() def harris_detect(img,ksize=3): ''' 自己实现角点检测 params: img:灰度图片 ksize：Sobel算子窗口大小 return： corner：与源图像一样大小，角点处像素值设置为255 ''' k = 0.04 #响应函数k threshold = 0.01 #设定阈值 WITH_NMS = False #是否非极大值抑制 #1、使用Sobel计算像素点x,y方向的梯度 h,w = img.shape[:2] #Sobel函数求完导数后会有负值，还有会大于255的值。而原图像是uint8，即8位无符号数，所以Sobel建立的图像位数不够，会有截断。因此要使用16位有符号的数据类型，即cv2.CV_16S。 grad = np.zeros((h,w,2),dtype=np.float32) grad[:,:,0] = cv2.Sobel(img,cv2.CV_16S,1,0,ksize=3) grad[:,:,1] = cv2.Sobel(img,cv2.CV_16S,0,1,ksize=3) #2、计算Ix^2,Iy^2,Ix*Iy m = np.zeros((h,w,3),dtype=np.float32) m[:,:,0] = grad[:,:,0]**2 m[:,:,1] = grad[:,:,1]**2 m[:,:,2] = grad[:,:,0]*grad[:,:,1] #3、利用高斯函数对Ix^2,Iy^2,Ix*Iy进行滤波 m[:,:,0] = cv2.GaussianBlur(m[:,:,0],ksize=(ksize,ksize),sigmaX=2) m[:,:,1] = cv2.GaussianBlur(m[:,:,1],ksize=(ksize,ksize),sigmaX=2) m[:,:,2] = cv2.GaussianBlur(m[:,:,2],ksize=(ksize,ksize),sigmaX=2) m = [np.array([[m[i,j,0],m[i,j,2]],[m[i,j,2],m[i,j,1]]]) for i in range(h) for j in range(w)] #4、计算局部特征结果矩阵M的特征值和响应函数R(i,j)=det(M)-k(trace(M))^2 0.04<=k<=0.06 D,T = list(map(np.linalg.det,m)),list(map(np.trace,m)) R = np.array([d-k*t**2 for d,t in zip(D,T)]) #5、将计算出响应函数的值R进行非极大值抑制，滤除一些不是角点的点，同时要满足大于设定的阈值 #获取最大的R值 R_max = np.max(R) #print(R_max) #print(np.min(R)) R = R.reshape(h,w) corner = np.zeros_like(R,dtype=np.uint8) for i in range(h): for j in range(w): if WITH_NMS: #除了进行进行阈值检测还对3x3邻域内非极大值进行抑制(导致角点很小，会看不清) if R[i,j] > R_max*threshold and R[i,j] == np.max(R[max(0,i-1):min(i+2,h-1),max(0,j-1):min(j+2,w-1)]): corner[i,j] = 255 else: #只进行阈值检测 if R[i,j] > R_max*threshold : corner[i,j] = 255 return corner if __name__=='__main__': img = cv2.imread('./image/cali.bmp') img = cv2.resize(img,dsize=(600,400)) #转换为灰度图像 gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) dst = harris_detect(gray) print(dst.shape) #(400, 600) img[dst>0.01*dst.max()] = [0,0,255] cv2.imshow('',img) cv2.waitKey(0) cv2.destroyAllWindows()
运行效果如下：

六 Harries角点的性质

1、参数

假设已经得到了矩阵

$d e t M = \prod_{i} λ_{i}$

$t r a c e (M) = \sum_{i} λ_{i}$

从而可以得到角点的响应：

$R = λ_{1} λ_{2} - k (λ_{1} + λ_{2})^{2} = λ_{1}^{2} (α - k (1 + α)^{2})$

假设

$0 \leq k \leq \frac{α}{(1 + α)^{2}} \leq 0.25$

对于较小的

由此，可以得出这样的结论，增大

2、Harris角点检测算子对亮度和对比度的变化不灵敏

这是因为在进行Harris角点检测时，使用了微分算子对图像进行微分运算，而微分运算对图像密度的拉升或收缩和对亮度的抬高或下降不敏感。换言之，对亮度和对比度的仿射变换并不改变Harris响应的极值点出现的位置，但是，由于阈值的选择，可能会影响角点检测的数量。

左图表示亮度变化，右图表示对比度变化。

3、Harris角点检测算子具有旋转不变性

Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值

4.、Harris角点检测算子不具有尺度不变性

如下图所示，当图像被缩小时，在检测窗口尺寸不变的前提下，在窗口内所包含图像的内容是完全不同的。左侧的图像可能被检测为边缘或曲线，而右侧的图像则可能被检测为一个角点。

参考文献：

[1]OpenCV教程(43) harris角的检测(1)

[2]Harris角点检测原理详解

[3]PCA的数学原理

[4]图像处理基础(4)：高斯滤波器详解

[5]Harris角点

[6]『Python』Harris角点检测代码实现

[7]斑点检测

[8]Harris原理及源码解析(C++)

https://zhuanlan.zhihu.com/p/83064609

Shi-Tomasi 角点检测器

知道了什么是 Harris 角点检测，后来有大佬在论文《Good_Features_to_Track》中提出了它的改进版——Shi-Tomasi 角点检测，Shi-Tomasi 方法在很多情况下可以得到比 Harris 算法更好的结果。

Harris 角点检测中每个窗口的分数公式是将矩阵 M 的行列式与 M 的迹相减：

$R = λ_1λ_2 − k (λ_1 + λ_2)^2 \$

由于 Harris 角点检测算法的稳定性和 k 值有关，而 k 是个经验值，不好设定最佳值。

Shi-Tomasi 发现，角点的稳定性其实和矩阵 M 的较小特征值有关，于是直接用较小的那个特征值作为分数。这样就不用调整k值了。

所以 Shi-Tomasi 将分数公式改为如下形式：

$R = min (λ_1, λ_2) \$

和 Harris 一样，如果该分数大于设定的阈值，我们就认为它是一个角点。

我们可以把它绘制到 λ1 ～ λ2 空间中，就会得到下图：

Shi-Tomasi 角点检测

OpenCV 提供了 Shi-Tomasi 的函数： cv2.goodFeaturesToTrack()，来获取图像中前 N 个最好的角点。

函数原型如下：
goodFeaturesToTrack(image, maxCorners, qualityLevel, minDistance[, corners[, mask[, blockSize[, useHarrisDetector[, k]]]]])
其中的参数如下：
- image：输入灰度图像，float32类型
- maxCorners：返回角点的最大数目，值为0表表示没有设置最大值限制，返回所有检测到的角点。
- qualityLevel：质量系数（小于1.0的正数，一般在0.01-0.1之间），表示可接受角点的最低质量水平。该系数乘以最好的角点分数（也就是上面较小的那个特征值），作为可接受的最小分数；例如，如果最好的角点分数值为1500且质量系数为0.01，那么所有质量分数小于15的角都将被忽略。
- minDistance：角之间最小欧式距离，忽略小于此距离的点。
- corners：输出角点坐标
- mask：可选的感兴趣区域，指定想要检测角点的区域。
- blockSize：默认为3，角点检测的邻域大小（窗口尺寸）
- useHarrisDetector：用于指定角点检测的方法，如果是true则使用Harris角点检测，false则使用Shi Tomasi算法。默认为False。
- k：默认为0.04，Harris角点检测时使用。
设定好这些参数，函数就能在图像上找到角点。所有低于质量水平的角点都会被忽略，然后再把合格角点按角点质量进行降序排列。

然后保留质量最高的一个角点，将它附近（最小距离之内）的角点都删掉（类似于非极大值抑制），按这样的方式最后得到 N 个最佳角点。
import numpy as np import cv2 maxCorners = 100 qualityLevel = 0.01 minDistance = 10 img = cv2.imread('bird.jpg') gray = cv2.cvtColor(img,cv2.COLOR_BGR2GRAY) corners = cv2.goodFeaturesToTrack(gray, maxCorners, qualityLevel, minDistance) corners = np.int0(corners) for i in corners: x,y = i.ravel() cv2.circle(img,(x,y),2,(0,0,255),-1) cv2.imwrite('new_bird.jpg', img)
最后得到的结果如下：

三张图分别为原图，harris图和Shi-Tomasi图

Harris 和 Shi-Tomasi 都是基于梯度计算的角点检测方法，Shi-Tomasi 的效果要好一些。基于梯度的检测方法有一些缺点: 计算复杂度高，图像中的噪声可以阻碍梯度计算。

想要提高检测速度的话，可以考虑基于模板的方法：FAST角点检测算法。该算法原理比较简单，但实时性很强。
查看全文

相关阅读:
C++文件(夹)选择对话框
 BCB中选择文件对话框TOpenDialog过滤后缀名使用方法
 pjlib深入剖析和使用详解
 PJNATH介绍 -- 开源的用于NAT穿透的ICE, STUN和TURN
STUN, TURN, ICE介绍
 一个boost底下的线程池
 在Windows下编译WebRTC
FEC（Forward Error Correction）前向纠错 UDPRTP 中使用用于改善无线等网络丢包等问题--转
 FEC之我见四
 FEC之异或运算应用