【图像处理】（3）harris corner detection(角点检测)

本小段引用内部的内容与本篇文章无关，仅仅是我个人的思考。

我在想如果是上述窗体内的梯度总体变化与相邻窗体的梯度的总体变化的差异程度也可以作为角点。

• 平坦区域就不会有什么变化，
• 边缘，某个方向上亮度变化剧烈，而垂直的方向上变化很小
• 角点，任意方向移动，平均亮度变化都很剧烈

1.角点特征的优越性：（WHY）

按照之前对优秀的局部特征的定义，角点特征具有可重复性和显著性特征，所以是优秀的局部特征。

2.角点检测的基础思想：（WHAT）

（1）窗口：

对局部特征检测来说，基础的方法是滑动窗口在图向上滑动，检测窗口滑动的内容（如所有像素值梯度）或者内容变化情况（如像素密度变化的整体情况），观察窗口捕捉到的信息是否符合局部特征的特点。这种方法很符合对“局部性”的阐释。

（2）角点区域的特征：

对角点来说，它是边的交点，边的特征是梯度在某一方向上产生突变，那么这个交点应该表现为某一区域内有两个或多个方向上的“边缘性变化”的梯度信息。需要强调的是，寻找角点的目标或许是找到一个“点”，但是分析的对象并不能只是一个“点”（在图像处理中，点指的是像素点），而应该是一个区域所有像素点的整体所呈现出来的特征（先记住这个“点”与“点所在区域”的区别，之后会再提到）。所以我们需要在一定范围内观察这个局部的整体，之前提到的窗口就起到了作用。

如果用滑动窗口来观察角点区域的话，窗口向多个方向滑动，都能感知到像素密度的强变化（梯度）。反观平坦区域和边缘区域，前者在任意滑动方向上不会产生强变化，后者只会在一个滑动方向上产生强变化。如下图：

3.角点检测的量化与算法构建：（HOW）

3.1.量化“窗口滑动下的变化特征”：

由角点在滑动窗口方法下的特征可知，我们需要度量的是：

窗口向某方向滑动前后，窗口内整体像素值的变化情况。

这里有两个重要的量化重点：

• 窗口的滑动方向
• 像窗口内整体素值变化

这两个量化重点之间的关系应该是自变量和因变量的关系，也就是说，我们需要建立一个函数，表达在某个滑动方向上窗口内整体像素值的变化情况。这里给出一个定义：

这是一个二元函数，其要素解释如下：

 

函数的物理意义在于：

因此，这个函数与我们的量化目标相符。

3.2.用量化特征识别角点区域

我们需要用量化函数来评估一块区域是否是角点区域，那么就必须了解不同类型的区域对这个函数的不同响应，从而知道角点区域对这个函数有什么特别的响应。

首先需要解构这个函数，将函数进一步改造，使之可以更直观地体现角点区域与其他类型区域的不同。

（1）等效变换：

我们在不改变函数的物理意义的情况下可以将函数变得更直观。

首先看二元函数的泰勒展开公式：

（2）在评价函数与局部特征间建立桥梁：

我们需要把抽象的评价函数与角点区域或其他类型区域的特征建立关系，这样我们才能将评价函数的功能发挥出来。

我们会发现这个函数的函数值好像只能告诉我们变化程度的大小，并不能告诉我们变化的类型。

那么我们的讨论对象当然就不能是函数的值的大小。

那我们应该讨论什么呢？

我们应该讨论函数的形状！也就是讨论M这一个系数对函数形状的影响。

（其实在这里，我就有一个疑惑：这个函数被定义的逻辑是什么？定义这个函数的时候，我怎么样才能想到要来讨论他的形状呢？如果有朋友看到了这篇文章，能否帮我解答一下这个问题。）

那为什么要讨论M呢？因为M的定义式中包含了x和y两个方向上的梯度信息，这两个梯度信息不正是之前分析角点区域特征时所提到的“像素密度的变化情况”吗？

之前，在分析角点区域特征时，我们提到：

对角点来说，它是边的交点，边的特征是梯度在某一方向上产生突变，那么这个交点应该表现为某一区域内有两个或多个方向上的“边缘性变化”的梯度信息。

我们注意到角点区域的像素会有“两个或多个方向”上的“边缘性变化”的梯度信息，刚好和这里的“x和y两个方向上的梯度信息”产生了联系，我们可以先假设角点区域只在x和y两个相互垂直的方向上有较为明显的变化，那么梯度信息就会有如下表现：

 

由上图可知，在“假设角点区域只在x和y两个相互垂直的方向上有较为明显的变化”之后，如果一个区域是角点区域，那么其梯度方向就会集中在x和y两个垂直的方向上，对于窗口内绝大多数像素点来说，要么y方向的梯度趋近0，要么x方向的梯度趋近0，两个方向的梯度的乘积也就趋近0了。对所有像素点进行累加，因为绝大多数像素点的“两个方向的梯度的乘积趋近0”，M矩阵经累加操作后，反对角线上的两个值也就趋近于0。

这样，M矩阵的信息就集中在了λ1和λ2两个值上，其中λ1反应了角点区域在x方向上的梯度信息，λ2反应了角点区域在y方向上的梯度信息。

这时有朋友可能会问两个问题，我分别用我的理解加以解释：

问：如果角的两边不垂直咋办？答：两边不垂直，但两边的梯度信息也应该会相对集中在两个垂直的方向上，只是如果两边越垂直，IxIy就越趋近于0。
问：如果角的两边不在x和y方向上咋办？答：首先，我们可做旋转变换，旋转变换是一种仿射变换，我们做角点检测也不在乎角点的方向。其次，对矩阵M进行“特征值分解”的操作，是不是就是一种旋转变换的体现呢，我线性代数没学扎实，等我学扎实了再来谈谈这个......

所以当角点区域的梯度集中在任意两个垂直的方向上时，我们对M进行特征值分解，总能得到以下结果：

物理意义上：

R表示一种旋转操作，决定了两个垂直方向的整体朝向。
λ1和λ2分别表示这两个垂直方向上的区域整体梯度信息

这个时候，我们回到那个评价函数E的形状上去，看看M怎样决定了这个形状，我们又一个怎样将函数形状的特点与角点区域特征对应起来。

如上图，若让函数E取一个函数值，就对得到一个方程，这个方程表示了E函数图像的水平切片，根据上图中那个不太严谨的变换，可以看出这是一个椭圆方程。

我们重点来分析这个椭圆切片的形状特征与角点区域特征有什么关系，如下图：

那么这个椭圆有啥用呢？

椭圆水平切片就是E函数的形状表现，具体如下：

如上图，这里对λ1和λ2的大小关系进行分类讨论，强调：λ1和λ2的大小反应两个垂直方向的梯度信息集中强度

1.如果两个值都大，且相差不大，则“E函数水平切片椭圆”趋近于圆，区域梯度在两个垂直的方向都有较强集中度，所以是角点区域。
2.如果其中一个值远大于另一个值，则，则“E函数水平切片椭圆”趋近于线，区域梯度只在一个方向上有较强集中度，所以是边缘区域
3.如果两个值都很小，则“E函数水平切片椭圆”趋近于点，区域梯度信息在两个方向都较弱，所以是像素值平坦区域

这样，我们就把评价函数与局部特征间的桥梁单间起来了。

之后，我们可以进一步用一个函数的值来综合表达λ1和λ2的三种情况，来省略分类讨论，如下图：

det(M)为行列式，trace(M)表示迹。下面的R仅仅是一种方式，主要是为了让公式契合下图的情况，仅此而已，你也可以用其他方式，只要能契合下述情况就好了。

我们将M称为“Second moment matrix”，将这个函数R称为“角点响应函数”。

总之：

1.某个区域的角点响应函数值大于0且绝对值较大的时候，我们判断为角点区域。
2.某个区域的角点响应函数值小于0且绝对值较大的时候，我们判断为边缘区域。
2.某个区域的角点响应函数绝对值较小的时候，我们判断为平坦区域。

同时，需要注意的是，角点响应函数不仅略去了分类讨论λ的过程，还将E函数的形状特征直接构建为一个函数，我们求解函数R的时候就可以忘记函数E了（R是多么可爱(*╹▽╹*)），这样，函数E物理意义中关于“窗口移动”的操作，也被消隐了，因为函数R与（u, v）无关，至于图像本身的梯度信息有关。

终于把角点区域的量化判定方法给理清了，之后就可以讨论具体算法了，这里，我们使用的算法是Harris角点检测。

4.Harris角点检测：

4.1.算法流程：

1. 计算整张图像每一个像素的高斯一阶导数（高斯平滑+求梯度的综合）
2. 用每一个像素的高斯一阶导数信息计算每一个像素的“Second moment matrix”M
3. 用M计算每一个像素点的角点响应函数值R
4. 设定R的阈值（阈值为一个较大的正数）
5. 判定角点区域，用非极大值抑制得到“角点”。

4.2.算法解析：

（1）高斯一阶导数：

这个操作也是老朋友了，在边缘检测中就有解释过了，图像上的噪声会影响梯度信息，所以要先用高斯滤波器进行滤波后再求导，等价于对高斯滤波器求一阶导得到高斯一阶导滤波器，直接用高斯一阶导滤波器对图像进行滤波（相关correlation或卷积convolution操作）得到梯度信息。

（2）“Second moment matrix”M和角点响应函数值R

之前我们提到，M是针对一个窗口区域进行构建，Harris算法中，对每个像素都计算了一个M，相当于是对所有“以每个像素为中心的区域”计算了M。

同样的，对每个像素点计算了一个R，就相当于对以每个像素点为中心的区域都评估了是角点区域的可能性。

 还记得一开始我们强调的那句话吗？

需要强调的是，寻找角点的目标或许是找到一个“点”（在图像处理中，点指的是像素点），但是分析的对象并不能只是一个“点”，而应该是一个区域所有像素点的整体所呈现出来的特征（先记住这个“点”与“点所在区域”的区别，之后会再提到）

（3）R的阈值与非极大值抑制：

我们设定了R的阈值Rt，R值大于该阈值的像素点所在区域被判定为角点区域，该像素点被判定为角点。

这里还进行了非极大值抑制操作，将局部区域内角点响应函数的局部最大值像素点作为这个区域的角点，同一个区域中，响应函数值比这个局部最大值小的像素点都被忽略掉，这体现了局部特征的“表达紧凑而高效”要求。

5.角点变换的不变性与协变性：

最后，我们讨论角点特征的不变性与协变性：

（1）仿射不变性：

整个图像的像素值都进行仿射变换后，角点信息是仿射不变的，如下图：

总体上看，角点响应函数R的总体大小发生了变化，但R在图像上的分布没有发生改变。

只是在阈值一定的情况下，会出现部分低于阈值的极值点在仿射变换后超过了阈值，被判定为新的角点。

（2）平移不变性：

如下图：

（3）旋转协变性

如下图：

原因是，M矩阵的特征值不变，梯度信息不变，只是旋转矩阵发生了变化，仍能通过特征值计算响应函数，识别角点。

以上不变性与协变性体现了局部特征的可重复性。

（4）尺度不协变

如下图：

图片尺度变化而检测窗口尺度不变时，会导致检测窗口对特征的检测结果不同，我们不想让这种结果出现，于是，我们需要一种检测窗口与图像的尺度协变策略。这是下一节要重点描述的内容。 

代码

 

#coding:utf-8

import numpy as np
import matplotlib.pyplot as plt
import os
import math

def convolve(filter,mat,padding,strides):
    '''
    :param filter:卷积核，必须为二维(2 x 1也算二维) 否则返回None
    :param mat:图片
    :param padding:对齐
    :param strides:移动步长
    :return:返回卷积后的图片。(灰度图，彩图都适用)
    @author:bilibili-会飞的吴克
    '''
    result = None
    filter_size = filter.shape
    mat_size = mat.shape
    if len(filter_size) == 2:
        if len(mat_size) == 3:
            channel = []
            for i in range(mat_size[-1]):
                pad_mat = np.pad(mat[:,:,i], ((padding[0], padding[1]), (padding[2], padding[3])), 'constant')
                temp = []
                for j in range(0,mat_size[0],strides[1]):
                    temp.append([])
                    for k in range(0,mat_size[1],strides[0]):
                        val = (filter*pad_mat[j:j+filter_size[0],k:k+filter_size[1]]).sum()
                        temp[-1].append(val)
                channel.append(np.array(temp))

            channel = tuple(channel)
            result = np.dstack(channel)
        elif len(mat_size) == 2:
            channel = []
            pad_mat = np.pad(mat, ((padding[0], padding[1]), (padding[2], padding[3])), 'constant')
            for j in range(0, mat_size[0], strides[1]):
                channel.append([])
                for k in range(0, mat_size[1], strides[0]):
                    val = (filter * pad_mat[j:j + filter_size[0],k:k + filter_size[1]]).sum()
                    channel[-1].append(val)


            result = np.array(channel)


    return result

def linear_convolve(filter,mat,padding=None,strides=[1,1]):
    '''
    :param filter:线性卷积核
    :param mat:图片
    :param padding:对齐
    :param strides:移动步长
    :return:返回卷积后的图片。(灰度图，彩图都适用) 若不是线性卷积核，返回None
    @author:bilibili-会飞的吴克
    '''
    result = None
    filter_size = filter.shape
    if len(filter_size) == 2 and 1 in filter_size:
        if padding == None or len(padding) < 2:
            if filter_size[1] == 1:
                padding = [filter_size[0]//2,filter_size[0]//2]
            elif filter_size[0] == 1:
                padding = [filter_size[1]//2,filter_size[1]//2]
        if filter_size[0] == 1:
            result = convolve(filter,mat,[0,0,padding[0],padding[1]],strides)
        elif filter_size[1] == 1:
            result = convolve(filter, mat, [padding[0],padding[1],0,0], strides)

    return result

def _2_dim_divided_convolve(filter,mat):
    '''
    :param filter: 线性卷积核,必须为二维(2 x 1也算二维) 否则返回None
    :param mat: 图片
    :return: 卷积后的图片,(灰度图，彩图都适用) 若不是线性卷积核，返回None
    '''
    result = None
    if 1 in filter.shape:
        result = linear_convolve(filter,mat)
        result = linear_convolve(filter.T,result)

    return result


def score_for_each_pixel(sq_img_gx,sq_img_gy,img_gx_gy,k):
    '''
    所传入的参数都只能有一个通道，且形状必须相同
    :param sq_img_gx: x方向上的梯度平方的图片
    :param sq_img_gy: y方向上的梯度平方的图片
    :param img_gx_gy: x,y方向上梯度乘积的图片
    :param k: 矩阵的迹前面的系数
    :return: 各点的得分
    '''
    result = []
    for i in range(sq_img_gx.shape[0]):
        result.append([])
        for j in range(sq_img_gx.shape[1]):
            M = np.array(
                [
                    [sq_img_gx[i,j],img_gx_gy[i,j]],
                    [img_gx_gy[i,j],sq_img_gy[i,j]]
                ]
            )
            result[-1].append(np.linalg.det(M)-k*(np.trace(M)**2))
    return np.array(result)

def Sign(img,score,area,decide_value=None,boder=[3,3,3,3]):
    '''
    :param img: 需要在角点处做标记的图片(可为多通道)
    :param score: 各个像素的角点得分
    :param area: 标记区域的大小(area[0] x area[1])
    :param decide_value: 决策是否为角点的阈值
    :param boder: 标记的边界宽度
    :return: 返回标记后的图片
    '''
    if decide_value == None:
        decide_value =  34*math.fabs(score.mean())  # 34这个参数可调
        print(decide_value)
    judger = score > decide_value
    final_img = img.copy()
    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            isLocalExtreme = score[i,j] >= score[max(i-(area[0]//2),0):min(i+(area[0]//2)+1,img.shape[0]),max(j-(area[1]//2),0):min(j+(area[1]//2)+1,img.shape[1])] #非极值抑制
            if judger[i,j] and isLocalExtreme.all():
                for k in range(min(boder[0],area[1]//6+1)):
                    final_img[max(i-(area[0]//2),0):min(i+(area[0]//2)+1,img.shape[0]),max(j-(area[1]//2),0)+k,:] = [255,0,0] #top
                for k in range(min(boder[1],area[1]//6+1)):
                    final_img[max(i-(area[0]//2),0):min(i+(area[0]//2)+1,img.shape[0]),min(j+(area[1]//2),img.shape[1]-1)-k,:] = [255,0,0] #bottom
                for k in range(min(boder[2],area[0]//6+1)):
                    final_img[max(i-(area[0]//2),0)+k,max(j-(area[1]//2),0):min(j+(area[1]//2)+1,img.shape[1]),:] = [255, 0, 0] #left
                for k in range(min(boder[3],area[0]//6+1)):
                    final_img[min(i+(area[0]//2),img.shape[0]-1)-k,max(j-(area[1]//2),0):min(j+(area[1]//2)+1,img.shape[1]),:] = [255,0,0]  # right
    return final_img

def OneDimensionStandardNormalDistribution(x,sigma):
    E = -0.5/(sigma*sigma)
    return 1/(math.sqrt(2*math.pi)*sigma)*math.exp(x*x*E)

if __name__ == '__main__':

    pic_path = './img/'
    pics = os.listdir(pic_path)

    window = 1.0/159*np.array([
        [2,4,5,4,2],
        [4,9,12,9,4],
        [5,12,15,12,5],
        [4,9,12,9,4],
        [2,4,5,4,2]
    ])   # window(5x5 Gaussisan kernel)
    linear_Gaussian_filter_5 = [2, 1, 0, 1, 2]
    sigma = 1.4
    linear_Gaussian_filter_5 = np.array([[OneDimensionStandardNormalDistribution(t, sigma) for t in linear_Gaussian_filter_5]])
    linear_Gaussian_filter_5 = linear_Gaussian_filter_5/linear_Gaussian_filter_5.sum()


    G_y = np.array(
        [
            [2, 2, 4, 2, 2],
            [1, 1, 2, 1, 1],
            [0 ,0 ,0 ,0 ,0],
            [-1,-1,-2,-1,-1],
            [-2,-2,-4,-2,-2]
        ]
    )
    G_x = np.array(
        [
            [2, 1, 0, -1, -2],
            [2, 1, 0, -1, -2],
            [4, 2, 0, -2, -4],
            [2, 1, 0, -1, -2],
            [2, 1, 0, -1, -2]
        ]
    ) #5x5 sobel kernel

    for i in pics:
        if i[-5:] == '.jpeg':
            orignal_img = plt.imread(pic_path+i)

            plt.imshow(orignal_img)
            plt.axis('off')
            plt.show()

            img = orignal_img.mean(axis=-1)

            img_gx = convolve(G_x,img,[2,2,2,2],[1,1])
            img_gy = convolve(G_y,img,[2,2,2,2],[1,1])

            sq_img_gx = img_gx * img_gx
            sq_img_gy = img_gy * img_gy
            img_gx_gy = img_gx * img_gy

            # sq_img_gx = convolve(window, sq_img_gx, [2, 2, 2, 2], [1, 1])
            # sq_img_gy = convolve(window, sq_img_gy, [2, 2, 2, 2], [1, 1])
            # img_gx_gy = convolve(window, img_gx_gy, [2, 2, 2, 2], [1, 1])

            sq_img_gx = _2_dim_divided_convolve(linear_Gaussian_filter_5,sq_img_gx)
            sq_img_gy = _2_dim_divided_convolve(linear_Gaussian_filter_5,sq_img_gy)
            img_gx_gy = _2_dim_divided_convolve(linear_Gaussian_filter_5,img_gx_gy)

            score = score_for_each_pixel(sq_img_gx,sq_img_gy,img_gx_gy,0.05)
            final_img = Sign(orignal_img,score,[12,12])

            plt.imshow(final_img.astype(np.uint8))
            plt.axis('off')
            plt.show()

参考

怎么深入了解 Harris 角点检测？ - 知乎

4.harris corner detection(角点检测)_哔哩哔哩_bilibili