【CV】Harris Corner Detection

已实现的功能简述及运行简要说明

功能简述

对输入的一张彩色图像，自己写代码实现Harris Corner检测算法，显示中间的处理结果及最终的检测结果，包括最大特征值图，最小特征值图，R图（可以考虑彩色展示），原图上叠加检测结果等，并将这些中间结果都输出成图像文件。

运行简要说明

原始图像位置在img_path变量中存储，将在origin窗口中进行显示；最大特征值图在lambda_max窗口中显示，在相同目录下存储为lambda_max.jpg；最小特征值图在lambda_min窗口中进行显示，在相同目录下存储为lambda_min.jpg；R图在R窗口中显示，在相同目录下存储为R.jpg；设置threshold之后的R图在R_threshold窗口中显示，存储为R_threshold.jpg；局部极大值点图在local_max窗口中显示，存储为local_max.jpg；原图叠加显示结果在cover窗口中显示，并存储为cover.jpg。

开发与运行环境

操作系统：Windows 10，64位

开发环境：Python 3.7.3

库环境：numpy 1.16.4, cv2 4.1.1

算法基本思路、原理

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

当窗口发生[u,v]移动时，那么滑动前与滑动后对应的窗口中的像素点灰度变化描述如下：

其中$[u,v]$是窗口的偏移量，$(x,y)$是窗口内所对应的像素坐标位置，窗口有多大，就有多少个位置。$w(x,y)$是窗口函数，有以下两种表示方法，最简单情形就是窗口内的所有像素所对应的w权重系数均为1。但有时候，我们会将$w(x,y)$函数设定为以窗口中心为原点的二元正态分布或高斯分布。如图所示：

根据上述表达式，当窗口处在平坦区域上滑动，可以想象的到，灰度不会发生变化，那么E(u,v) = 0；如果窗口处在比纹理比较丰富的区域上滑动，那么灰度变化会很大。算法最终思想就是计算灰度发生较大变化时所对应的位置，当然这个较大是指针任意方向上的滑动，并非单指某个方向。

我们对$E(w,v)$表达式进行泰勒展开，结果如下：

可以将$E(u,v)$表达式更新为：

其中

求解M的两个特征值$\lambda_1$和$\lambda_2$，可以根据两个特征值直接的关系得出对应像素点的位置信息（边缘or角点or平坦区域）

corner：在水平、竖直两个方向上变化均较大的点，即Ix、Iy都较大；

edge ：仅在水平、或者仅在竖直方向有较大的点，即Ix和Iy只有其一较大；

flat ：在水平、竖直方向的变化量均较小的点，即Ix、Iy都较小；

可以定义角点响应函数R来表示：

其中k是介于0.04-0.06间的常数。

针对三种不同区域的点，R的取值情况如下：

corner：R为大数值整数

edge：R为大数值负数

flat：绝对值R是小数值

具体实现

梯度计算

计算各像素点位置的x、y方向梯度值及梯度方向，实现方式如下：

for i in range(width-1):
    for j in range(height-1):
        ix[i,j] = (int(img[i,j]) - int(img[i+1,j]) + int(img[i,j+1]) - int(img[i+1,j+1])) / 2
        iy[i,j] = (int(img[i,j]) - int(img[i,j+1]) + int(img[i+1,j]) - int(img[i+1,j+1])) / 2
        if ix[i,j] != 0:
            theta[i,j] = np.arctan(iy[i,j] / ix[i,j]) * 180 / np.pi
            if(theta[i,j] < 0):
                theta[i,j] = theta[i,j] + 180
        else:
            theta[i,j] = 90

M矩阵构建及计算

根据原理中的公式，对前一步计算得到的x、y方向梯度值进行计算，得到 $I_X^2$ 、$I_Y^2$ 和 $I_{XY}$，并进行高斯滤波，减少噪声带来的影响。

1
2
3

ixx = cv2.GaussianBlur(ix*ix, (5,5), 1.4) #ix^2
iyy = cv2.GaussianBlur(iy*iy, (5,5), 1.4) #iy^2
ixy = cv2.GaussianBlur(ix*iy, (5,5), 1.4) #ix*iy

计算M的行列式、迹、特征值

detM = ixx * iyy - ixy * ixy #行列式
traceM = ixx + iyy #迹
lambda_max = (traceM + np.sqrt(np.abs(traceM**2-4*detM)))/2 #较大特征值
lambda_min = (traceM - np.sqrt(np.abs(traceM**2-4*detM)))/2 #较小特征值

角点响应函数

计算角点响应函数，并进行处理，过滤掉R值小于0（即平坦区域的点）

R = detM - k * traceM**2
for i in range(width-1):
    for j in range(height-1):
        if R[i,j] < 0:
            R[i,j] = 0

设置阈值，去除边缘点，只保留角点

threshold = np.max(R) / 100000
for i in range(width-1):
    for j in range(height-1):
        if R[i][j] > threshold:
            Rthreshold[i][j] = R[i][j]
        else:
            Rthreshold[i][j] = 0

局部极大值

在之前选择的角点中进一步进行筛选，只保留局部极大值点

for i in range(15, width-15):
    for j in range(15, height-15):
        if R[i,j] == np.max(R[i-15:i+15, j-15:j+15]) and R[i][j] > threshold:
            local_max[i][j] = 255
        else:
            local_max[i][j] = 0

叠加显示

在原图上叠加显示角点检测结果

for i in range(15,width-15):
    for j in range(15,height-15):
        if local_max[i,j] == 255:
            img_cover[i-2:i+2,j-2:j+2,0] = 255
            img_cover[i-2:i+2,j-2:j+2,1] = 0
            img_cover[i-2:i+2,j-2:j+2,2] = 0