【CV】Harris角点检测算法-点特征提取
1 理论部分
1.1 特征
在数据挖掘、计算机视觉等领域经常会提起“特征”一词，那么究竟什么才是特征？数据挖掘中我们可以称被分析数据集的某一属性（一行为一个样本，一列为一个属性）为特征，因为它反映了该数据集的某一种特性，是该数据集的一种固有标识。在计算机视觉中，一张图像的特征其实也有类似的概念，这种特征也能够描述这张图像的全部（或局部）内容或含义。
图像特征可以是①有意义的图像区域 ，该区域具有独特的特征和良好的识别性；也可以是②对图像信息的一种数据抽取和表示 。
- 针对①可以是图像中的角点、边缘、斑点、直线、曲线及高密度区域等。大多数特征检测算法都会涉及图像的角点、边和斑点的识别，也有一些涉及脊向的概念，可以认为脊向是细长物体的对称轴，例如识别图像中的一条路。角点和边都好理解，那什么是斑点呢？斑点通常是指与周围有着颜色和灰度差别的区域。在实际地图中，往往存在着大量这样的斑点，如一颗树是一个斑点，一块草地是一个斑点，一栋房子也可以是一个斑点。由于斑点代表的是一个区域，相比单纯的角点，它的稳定性要好，抗噪声能力要强，所以它在图像配准上扮演了很重要的角色。
- 对于②，则是采用直方图向量、投影向量等方法提取出新形式的数据特征，使得所提取出来的特征具有一些较好的特性，比如旋转不变性、光照不变性等，使得模型能更好地检测或识别出图像的特征。
图像的特征应用很广泛，我们通过图像特征能过够将实现图像的分类、目标检测、图像拼接等。而特征就需要特征检测来提出，最终达到的目标是：通过特征检测我们实现的就是能够让计算机“读懂”图像，将人类看到的视觉信息转化成计算机能够识别和处理的定量形式。
1.2 角点
角点是常见的点特征，它具有旋转不变形，光照不变性，视角不变性等特性，可以在目标匹配、目标跟踪、三维重建等应用中使用。常见的点特征提取算法主要有Harris角点检测算法、FAST特征检测算法、SIFT特征检测算法、SURF特征检测算法等，本文主要介绍Harris算法。
用一个滑动的窗口在图像内移动，分别会出现以下三种情况（如图1所示）：
- 平坦区域 （左图）：当窗口向任意方向移动时，窗口内的像素值不会发生太大变化；
- 边缘 （中图）：当窗口沿着平行于边缘的方向移动时，窗口内的像素值不会发生太大变化；
- 角点 （右图）：不管窗口向着哪个方向移动，窗口内的像素值都会发生很大的改变。
图1. 角点的识别
1.3 图像梯度
图像的梯度，指的是图像灰度值的变化率。由于图像灰度变化并不连续，因此采用差分法来求梯度。如点(x,y)处在x方向和y方向的梯度按sobel算子计算，具体见Sobel算子简介。
1.4 Harris角点检测算法思想
该算法可以分为三个主要的步骤：
-
建立数学模型描述窗口灰度变化
当窗口在x和y方向分别移动了u和v该窗口的灰度值变化为I(x+u,y+v)-I(x,y)。
设w(x,y)为窗口函数，即窗口内各点之间的权重，表示该点对窗口灰度变化的贡献大小。当各点权重都为1时，此窗口函数也称为一个均值滤波核 ：
\left[ \begin{matrix} 1&1&1\\ 1&1&1\\ 1&1&1 \end{matrix} \right]当窗口函数为以窗口中心为原点并呈高斯分布时，此窗口函数也称为高斯滤波核 ，如：
\frac 1 {16} \left[ \begin{matrix} 1&2&1\\ 2&4&2\\ 1&2&1 \end{matrix} \right]高斯滤波实质是一种加权平均滤波，为了实现平均，核还带有一个系数，例如上面的十六分之一，这些系数等于矩阵中所有数值之和的倒数。因为高斯滤波是以空间距离来确定权重大小的，但并没有考虑用颜色距离来确定权重，这样就导致了高斯滤波在去除噪声的同时，也一定程度上模糊了边界。关于高斯滤波可参考图像平滑去噪之高斯滤波器。
图2. 移动窗口导致的灰度变化
将图2中公式进行简化，由泰勒公式得：
I(x+u,y+v)=I(x,y)+uI_x+vI_y I(x+u,y+v)-I(x,y)=uI_x+vI_y
窗口的灰度变化E(u,v)：
E(u,v)=\Sigma_{(x,y)} w(x,y)\times(u^2I_x^2+v^2I_y^2+2uvI_xI_y)
将上式写成矩阵形式：
E(u,v)=\sum_{(x,y)} w(x,y)\times\left[ \begin{matrix} u & v \end{matrix} \right] \left[ \begin{matrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{matrix} \right] \left[ \begin{matrix} u \\ v \end{matrix} \right]=\left[ \begin{matrix} u & v \end{matrix} \right] \boldsymbol{M} \left[ \begin{matrix} u \\ v \end{matrix} \right] \boldsymbol{M}矩阵为一个是实对称矩阵，将其对角化后，可得：
\boldsymbol{M}=\sum_{(x,y)} w(x,y)\times\left[ \begin{matrix} I_x^2 & I_xI_y \\ I_xI_y & I_y^2 \end{matrix} \right]=R\left[ \begin{matrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{matrix} \right]R^{-1}其中，R为正交矩阵，R^T=R^{-1}，因此E(u,v)=\left[ \begin{matrix} u & v \end{matrix} \right] R \left[ \begin{matrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{matrix} \right] R^{T} \left[ \begin{matrix} u & v \end{matrix} \right]^T=\left[ \begin{matrix} u' & v' \end{matrix} \right] \left[ \begin{matrix} \lambda_1 & 0 \\ 0 & \lambda_2 \end{matrix} \right] \left[ \begin{matrix} u' \\ v' \end{matrix} \right]=\lambda_1 ({u'})^2+\lambda_2 ({v'})^2=\frac {({u'})^2} {\frac 1 {\lambda_1}} + \frac {({v'})^2} {\frac 1 {\lambda_2}}我们得到 E(u,v)就是一个椭圆，椭圆的长轴为\lambda_2^{-0.5}，短轴为\lambda_1^{-0.5}（假设\lambda_1>\lambda_2）。因为u和v为每次窗口移动的位移，要使E(u,v)最大，则\lambda_1和\lambda_2必须同时最大。
我们还可以从另一个角度来理解：矩阵代表着运动，且当矩阵维度没有发生改变时，一个矩阵就包含了两个运动：旋转和拉伸。其中，对角化得到的正交矩阵 R可看做旋转因子 ，对角矩阵 包含的两个特征值λ_1和λ_2可看做缩放因子 ，也即是两个正交方向上的变化分量。椭圆的大小与旋转因子无关，只与缩放因子（即特征值）有关。
上述旋转和缩放过程可以直观地放在椭圆上观察。下图中，对角矩阵\left[ \begin{matrix} 3 & 0\\ 0 & 1 \end{matrix} \right]为包含两个特征值3和1，决定了红色椭圆短轴和长轴的长度。若在对角矩阵的基础上与正交矩阵相乘，即乘以旋转因子，则得M=\left[ \begin{matrix} 2 & -1 \\ -1 & 2 \end{matrix} \right]，红色椭圆也随之旋转为蓝色椭圆。
图3.椭圆的旋转 -
构造角点响应函数 R
窗口灰度的变化可由E(u,v)表示，而E(u,v)的大小主要由矩阵\boldsymbol{M}决定。上面提到，要使E(u,v)最大，则矩阵\boldsymbol{M}的特征值必须同时最大。
为了在实际中更好地应用于编程，我们构造角点响应函数R：
R=\text{det}(\boldsymbol{M})-k(\text{trace}(\boldsymbol{M}))^2 \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ (1)其中，\text{det}(\boldsymbol{M})=λ_1λ_2为矩阵的行列式，\text{trace}(\boldsymbol{M})=λ_1+λ_2为矩阵的迹，k为一个经验常数，通常取在0.04~0.06之间。
我们可以看出，当特征值同时取最大时，R才取最大，所以可以用R来代替E(u,v)。
另外还有其他形式的角点响应函数R：
R=\lambda_minR = \lambda_1-k\lambda_2,\ \ \ k=0.05R=\frac {det(M)} {trace(M)}=\frac {\lambda_1\lambda_2} {\lambda_1+\lambda_2}关于这些角点函数，本文不做过多的讨论，下文中R所指均为式（1）。 -
角点的判定
根据R值的大小，可以将窗口所在的区域划分为平面、边缘和角点（图3）：
a. 平面 ：由于灰度值变化较小，所以R值很小，并且沿任意方向的灰度变化都很小，所以I_x和I_y都较小，所以λ_1和λ_2也都较小；
b. 边缘 ：R为负数，仅在垂直于边缘的方向上有较大的变化量，I_x和I_y中只有一个变化较大，因此λ_1\ggλ_2或者λ_1\llλ_2；另外，对于斜边缘，I_x和I_y可能都比较大，但I_x和I_y成一定的比例，同样对应λ_1\ggλ_2或者λ_1\llλ_2；
c. 角点 ：R值很大，并且I_x和I_y都很大，所以λ_1和λ_2也都较大。
更详细的解释可见2.2.3 矩阵MM的理解。
图4. 角点的判定
1.5 Harris算法的特性
参数k对对角点检测结果的影响
设矩阵\boldsymbol{M}的特征值\lambda_1\ge\lambda_2\ge0，令\lambda_2=n\lambda_1,0\le n\le 1，则角点响应为：
R=\lambda_1\lambda_2-k(\lambda_2+\lambda_2)^2=\lambda_1^2(n-k(1+n)^2)假设R\ge0，则有：
0\le k\le\frac{n}{(1+n)^2}\le0.25对于较小的n值，R\approx\lambda_1^2(n-k),k
所以，可以得出这样的结论：
增大k的值，将减小角点响应值R，降低角点检测的灵性，减少被检测角点的数量；
减小k值，将增大角点响应值R，增加角点检测的灵敏性，增加被检测角点的数量。
对亮度和对比度的变化不敏感
Harris算法采用图像梯度或微分来进行运算，而梯度或微分对图像密度的拉升或收缩，以及对亮度的升高或降低不明感。也就是说，对亮度和对比度的仿射变换并不改变极值点出现的位置。但是阈值选择的不同会影响检测角点的数量。
图5. 部分不变性
旋转不变性
Harris角点检测算子使用的是角点附近的区域灰度二阶矩矩阵。而二阶矩矩阵可以表示成一个椭圆，椭圆的长短轴正是二阶矩矩阵特征值平方根的倒数。当特征椭圆转动时，特征值并不发生变化，所以判断角点响应值R也不发生变化，由此说明Harris角点检测算子具有旋转不变性。
图6. 用椭圆表示二阶矩矩阵
不具有尺度不变性
当图像被缩小时，所检测出的角点可能是完全不一样的。
图7. 不具有尺度不变性
2 实践部分
OpenCV中cornerHarris函数
OpenCV中有可直接调用的Harris角点检测的函数：
cv2.cornerHarris(src, blockSize, ksize, k[, dst[, borderType]])参数解释：
- src – 输入的灰度图像，float32类型；
- blockSize - 用于角点检测的邻域大小，就是上面提到的窗口的尺寸；
- ksize - 用于计算梯度的Sobel算子的尺寸；
- k - 用于计算角点响应函数的参数k，取值在0.04~0.06之间；
- dst - 存储着Harris角点响应的图像矩阵，大小与输入图像大小相同，是一个浮点型矩阵；
- borderType - 边界处理的类型。
下面基于python进行Harris角点检测：
import cv2 as cv
from matplotlib import pyplot as plt
import numpy as np
# detector parameters
block_size = 3
sobel_size = 3
k = 0.06
image = cv.imread('Building.jpg')
print(image.shape)
height = image.shape[0]
width = image.shape[1]
channels = image.shape[2]
print("width: %s height: %s channels: %s" % (width, height, channels))
# Convert BGR to Gray
gray_img = cv.cvtColor(image, cv.COLOR_BGR2GRAY)
# modify the data type setting to 32-bit floating point
gray_img = np.float32(gray_img)
# detect the corners with appropriate values as input parameters
corners_img = cv.cornerHarris(gray_img, block_size, sobel_size, k)
# result is dilated for marking the corners, not necessary
kernel = cv.getStructuringElement(cv.MORPH_RECT, (3, 3))
# 取核中像素值的最大值，代替核的中心位置（锚点）处的像素值
dst = cv.dilate(corners_img, kernel)
# Threshold for an optimal value, marking the corners in Red
for r in range(height):
for c in range(width):
pix = dst[r, c]
if pix > 0.05 * dst.max():
cv.circle(image, (c, r), 5, (0, 0, 255), 0)
image = cv.cvtColor(image, cv.COLOR_BGR2RGB)
plt.imshow(image)
plt.show()
图8. 角点检测结果
2.2 基于python实现Harris算法
def harris_det(img, block_size=3, ksize=3, k=0.04, threshold = 0.01, WITH_NMS = False):
'''
params:
img:单通道灰度图片
block_size:权重滑动窗口
ksize：Sobel算子窗口大小
k:响应函数参数k
threshold:设定阈值
WITH_NMS:非极大值抑制
return：
corner：角点位置图，与源图像一样大小，角点处像素值设置为255
'''
h, w = img.shape[:2]
# 1.高斯权重
gray = cv2.GaussianBlur(img, ksize=(ksize, ksize), sigmaX=2)
# 2.计算梯度
grad = np.zeros((h,w,2),dtype=np.float32)
grad[:,:,0] = cv2.Sobel(gray,cv2.CV_16S,1,0,ksize=3)
grad[:,:,1] = cv2.Sobel(gray,cv2.CV_16S,0,1,ksize=3)
# 3.计算协方差矩阵
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]
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__":
image = cv2.imread('./timg.jpg')
height, width, channel = image.shape
print('image shape --> h:%d w:%d c:%d' % (height, width, channel))
cv2.imshow('mount', image)
cv2.waitKey(2000)
cv2.destroyAllWindows()
gray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
# gray = np.float32(gray)
dst = harris_detect(gray)
image_dst = image[:, :, :]
image_dst[dst > 0.01 * dst.max()] = [0, 0, 255]
cv2.imwrite('./dsti8_1.jpg', image_dst)
cv2.imshow('dst', dst)
cv2.waitKey(0)
cv2.destroyAllWindows()