角点检测(2) - harris算子 - 理论与Python代码

数字图像，图像=矩阵，[m*n]从[0,255]的灰度值

角点检测：物体边缘的拐点

->应用：图像匹配与检索、图像物体形变恢复（摄像机标定）、三维重建

Harris角点检测（早期，原理简单，视频跟踪，快速检测）

• 1988年-Harris+Stepehns
• 图像中一个局部窗口在不同方向少量平移后窗口内图像灰度值的平均变化

• 动机：特征点具有局部差异性
• 以每个点为中心取一个窗口，例如，5×5/7×7的像素，描述特征点周围环境
• 此点具有差异性->窗口往任意方向移动，则周围环境变化较大->具有局部差异性

分类：

• 光滑 - flat
• 边缘 - edge
• 角点 - corner

一个方向上的灰度变化

E(u,v)=\sum_{x,y} w(x,y)[I(x+u,y+v)-I(x,y)]^{2} 越大越可能为角点

• 竖直和水平方向的偏移，
• w(x,y)窗口中心，高斯滤波/均值滤波
• x+u,y+v中心加偏移，
• x,y中心灰度值

改进——多个反向上的像素变化：

• 泰勒展开公式f(x+u,y+v)\approx f(x,y)+uf_x(x,y)+vf_y(x,y)

所以：

\\ \begin{equation} \begin{split} E(u,v)&=\sum_{(x,y)} w(x,y)[I(x+u,y+v)-I(x,y)]^2 \\ &\approx \sum_{(x,y)} w(x,y)[I(x,y)+\frac{ \partial I }{ \partial x}(x,y)u+\frac{ \partial I }{ \partial y}(x,y)v-I(x,y)]^2 (一阶泰勒展开)\\ &\approx \sum_{(x,y)} w(x,y)[\frac{ \partial I }{ \partial x}(x,y)u+\frac{ \partial I }{ \partial y}(x,y)v ]^2 （消除重复项）\\ &=\sum_{x,y} w(x,y)(u^2 f_x^2(x,y)+2uvf_x(x,y)f_y(x,y)+v^2f_y^2(x,y)) (二阶泰勒展开)\\&=\sum_{x,y} w(x,y)(u^2I_x ^2+2uv I_x I_y+v^2I_y ^2) （简化） \\ &=\left[ \begin{array}{ccc}u & v\\\end{array}\right ](\sum w(x,y)\left[ \begin{array}{ccc}I_x^2 & I_xI_y \\ I_xI_y & I_y^2\end{array}\right ])\left[ \begin{array}{ccc}u\\v\end{array}\right ] \\ &= \left[\begin{matrix} u \\ v \\ \end{matrix} \right] ^T H \left[\begin{matrix} u \\ v \\ \end{matrix} \right] \\ \end{split} \end{equation}

• 加和符号：表示窗口内每个像素
• w：表示权重，权值1或者以点为中心的高斯权重（离点越近权重越大）
• I：表示像素，RGB/灰度
• u，v：窗口移动的方向
• H：harris矩阵，由两个方向上的梯度构建而成
• 图像梯度： \Delta I(x,y)=(\frac{\partial I}{\partial x}(x,y),\frac{\partial I}{\partial y}(x,y))
• Harris矩阵： H= \left[ \begin{matrix} \sum_{(x,y)}w(x,y)(\frac{\partial I}{\partial x}(x,y))^2 & \sum_{(x,y)}w(x,y)(\frac{\partial I}{\partial x}(x,y))\frac{\partial I}{\partial y}(x,y)) \\ \sum_{(x,y)}w(x,y)(\frac{\partial I}{\partial x}(x,y))\frac{\partial I}{\partial y}(x,y)) & \sum_{(x,y)}w(x,y)(\frac{\partial I}{\partial y}(x,y))^2\\ \end{matrix} \right]
• 图像水平梯度、图像垂直梯度

I_x =\frac{ \partial I(x+u,y+v) }{\partial x} I_y =\frac{ \partial I(x+u,y+v) }{\partial y}

Harris矩阵H 的特征值分析

• 两个特征值反映相互垂直方向上的变化情况，分别代表变化最快和最慢的方向，特征值大变化快，特征值小变化慢
• SVD(H)=U\sum V, (\lambda_1, \lambda_2), \ \ \lambda_1>\lambda_2
• λ1 ≈ λ2 ≈ 0， 两个方向上变化都很小，兴趣点位于光滑区域
• λ1 > 0 , λ2 ≈ 0 ，一个方向变化快，一个方向变化慢，兴趣点位于边缘区域
• λ1 , λ2 > 0 ， 两个方向变化都很快，兴趣点位于角点区域（容易判断）

Harris角点准则代替矩阵分解： C=det(H)-ktrace(H)^2=\lambda_1\lambda_2-k(\lambda_1+\lambda_2)^2, k=0.04

• 反映特征值情况，trace为迹
• k的值越小，检测子越敏感
• 只有当λ1和λ2同时取得最大值时，C才能取得较大值
• 避免了特征值分解，提高检测计算效率
• 非极大值抑制(Non-maximal Suppression) 选取局部响应最大值，避免重复的检测

算法流程：

• 0）滤波、平滑，避免出现阶跃函数
• 1）计算图像水平和垂直方向的梯度 \frac{\partial L(x,y,\sigma _D)}{\partial x}=\frac{\partial G(x,y,\sigma _D)}{\partial x} * I(x,y), \frac{\partial L(x,y,\sigma _D)}{\partial y}=\frac{\partial G(x,y,\sigma _D)}{\partial y} * I(x,y)
• 2）计算每个像素位置的Harris矩阵 H= G(x,y, \sigma_1)*\left[ \begin{matrix} \sum_{(x,y)}w(x,y)(\frac{\partial L}{\partial x}(x,y))^2 & \sum_{(x,y)}w(x,y)(\frac{\partial L}{\partial x}(x,y))\frac{\partial L}{\partial y}(x,y)) \\ \sum_{(x,y)}w(x,y)(\frac{\partial L}{\partial x}(x,y))\frac{\partial L}{\partial y}(x,y)) & \sum_{(x,y)}w(x,y)(\frac{\partial L}{\partial y}(x,y))^2\\ \end{matrix} \right]
• 3）计算每个像素位置的Harris角点响应值 C=det(H)-ktrace(H)^2=\lambda_1\lambda_2-k(\lambda_1+\lambda_2)^2
• 3+）非极大值抑制
• 4）找到Harris角点响应值大于给定阈值且局部最大的位置作为特征点

检测结果：

求导：使用Sobel梯度算子

和sobel算子做卷积计算，3*3与 3*3卷积改写中心值

Harris确定角点

判断：

Python代码

import numpy as np
import cv2

def cornerHarris(img, blocksize=2, ksize=3, k=0.04):
# ksize是sobel的大小，blocksize是滑动窗口的大小
    def _clacHarris(cov,k):
        result = np.zeros([cov.shape[0], cov.shape[1]], dtype=np.float32)
        for i in range(cov.shape[0]):
            for j in range(cov.shape[1]):
                a = cov[i,j,0]
                b = cov[i,j,1]
                c = cov[i,j,2]
                result[i,j] = a*c - b*b - k*(a+c)*(a+c)
        return result

    Dx = cv2.Sobel(img,cv2.CV_32F,1,0,ksize=ksize)
    Dy = cv2.Sobel(img,cv2.CV_32F,0,1,ksize=ksize)

    cov = np.zeros([img.shape[0], img.shape[1], 3], dtype=np.float32)

    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            cov[i,j,0] = Dx[i,j] * Dx[i,j]
            cov[i,j,1] = Dx[i,j] * Dx[i,j]
            cov[i,j,2] = Dy[i,j] * Dy[i,j]
    #算块内的梯度和，
    #和[1,1 
    #   1,1]相乘
    cov = cv2.boxFilter(cov, -1, (blocksize, blocksize), normalize = False)
    return _clacHarris(cov,k)

if __name__ == '__main__':
    img = cv2.imread('test.jpg')
    gray_img = cv2.cvtColor(img, cv2.COLOR_BGR2GRAY)
    result = cornerHarris(gray_img, 2, 3, 0.04)
    pos = cv2.goodFeaturesToTrack(result, 0, 0.01, 10)
    for i in range(len(pos)):
        cv2.circle(img, (pos[i][0][0], pos[i][0][1]), 5, [255,0,0], thickness=3)
    cv2.imshow('harris',img)
    cv2.waitKey(0)