Harris Corner Detection

这个博客写的不详细，有兴趣的话可以看视频：https://www.bilibili.com/video/BV1Yo4y1Z7sD。

声音可能不太好听，是用 TTS 生成的。

使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，

都有着较大灰度变化，那么我们就可以认为该窗口中存在角点。举个例子，下图中 $A$ 是角点，$B$ 不是角点，如果用一个小窗口(蓝色)在点 $A$

处滑动，无论朝哪个方向走，窗口内的像素值都会发生明显的变化，比如窗口往左，那白色像素就越多，红色越少；但是对于 $B$ 点，无论朝哪

个方向滑动，窗口内的像素基本不变。

    

进一步地

      

基本数学公式如下：

$$E(u,v) = sum_{(x,y) in W}^{}W(x,y)left [ I(x+u,y+v) - I(x,y) ight ]^{2}$$

其中 $w(x, y)$ 表示移动窗口，$I(x, y)$ 表示像素灰度值强度，范围为 $0-255$，$x,y$ 表示窗口内的像素坐标，$u,v$ 表示方向。

函数 $E(u,v)$ 是和图像一个窗口区域一一对应的，即首先需要用 $W(x,y)$ 确定一个窗口，然后再计算这个图像区域对应的 $E(u,v)$

函数来判断该窗口内是否存在角点。根据二元泰勒展开式：

$$I(x+u,y+u) = I(x,y) + uI_x(x,y) + vI_y(x,y) + R_n$$

去掉高阶项，原式近似为

$$E(u,v) approx sum_{(x,y) in W}^{}W(x,y)left [ uI_x(x,y) + vI_y(x,y) ight ]^{2} 
= sum_{(x,y) in W}^{}W(x,y) left (u^2I_x^2 + 2uvI_xIy + v^2I_y^2 ight ) \
= sum_{(x,y) in W}^{}W(x,y)egin{bmatrix}
u & v 
end{bmatrix}egin{bmatrix}
I_x^2 & IxIy\ 
IxIy & I_y^2
end{bmatrix}egin{bmatrix}
u \ 
v
end{bmatrix} \ 
= egin{bmatrix}
u & v 
end{bmatrix}left (sum_{(x,y) in W}^{}W(x,y)egin{bmatrix}
I_x^2 & IxIy\ 
IxIy & I_y^2
end{bmatrix} ight ) egin{bmatrix}
u \ 
v
end{bmatrix} \
= 
egin{bmatrix}
u & v 
end{bmatrix}M egin{bmatrix}
u \ 
v
end{bmatrix}$$

每个实对称矩阵累加之后还是一个实对称矩阵，所以得到的 $M$ 是一个实对称矩阵。进一步有：

$$E(u,v) approx egin{bmatrix}
u & v 
end{bmatrix}Pegin{bmatrix}
lambda_1 & \ 
& lambda_2
end{bmatrix}P^T egin{bmatrix}
u \ 
v
end{bmatrix} \
= egin{bmatrix}
u & v 
end{bmatrix}Pegin{bmatrix}
lambda_1 & \ 
& lambda_2
end{bmatrix} left (egin{bmatrix}
u & v
end{bmatrix}P ight )^T \
= egin{bmatrix}
u^{'} & v^{'}
end{bmatrix}egin{bmatrix}
lambda_1 & \ 
& lambda_2
end{bmatrix} egin{bmatrix}
u^{'} \ v^{'}
end{bmatrix} \
= lambda_1u^{'2} + lambda_2v^{'2} \
= frac{u^{'2}}{1/lambda_1} + frac{v^{'2}}{1/lambda_2}$$

可以看出函数 $E(x,y)$ 在 $xoy$ 平面的等高线就是椭圆。

我们的目标是：无论 $u,v$ 取什么值，或者说窗口 $w$ 无论沿哪个方向运动，$E(u,v)$ 都要尽可能地大。这样才符合角点的性质。

也就是说如果一个区域存在角点，那么相比与其它区域，相同的高度对应的等高线(椭圆)越小。

所以分母要尽可能地小，即特征值 $lambda_1,lambda_2$ 要尽可能地大。

    

当 $lambda_1 >> lambda_2$ 或 $lambda_2 >> lambda_1$ 时，检测到的就是边；当 $lambda_1,lambda_2$ 都很小，$E$ 接近常数，检查测到的是平坦区域；两个都很大，

并且近似一样大时，检测到的就是角点。引入一个量化的指标来衡量这种关系：

$$R = det M - k(trace ; M) = lambda_1lambda_2 - left ( lambda_1 + lambda_2 ight )$$

代码实现过程主要在于求得矩阵 $M$，即

$$M = left (sum_{(x,y) in W}^{}w(x,y)egin{bmatrix} I_x^2 & IxIy\  IxIy & I_y^2 end{bmatrix} ight )$$

假设窗口大小为 $5 imes 5$，即

$$W = egin{bmatrix}
w_{11} & w_{12} & w_{13} & w_{14} & w_{15} \ 
w_{21} & w_{22} & w_{23} & w_{24} & w_{25} \ 
w_{31} & w_{32} & w_{33} & w_{34} & w_{35} \ 
w_{41} & w_{42} & w_{43} & w_{44} & w_{45} \ 
w_{51} & w_{52} & w_{53} & w_{54} & w_{55}
end{bmatrix}$$

那么矩阵 $M$ 可以表示为

$$M = egin{bmatrix}
w_{11}I_x^2 & w_{11}I_xI_y \ 
w_{11}I_xI_y & w_{11}I_y^2
end{bmatrix} + egin{bmatrix}
w_{12}I_x^2 & w_{12}I_xI_y \ 
w_{12}I_xI_y & w_{12}I_y^2
end{bmatrix} + cdots + egin{bmatrix}
w_{55}I_x^2 & w_{55}I_xI_y \ 
w_{55}I_xI_y & w_{55}I_y^2
end{bmatrix}$$

因为不好表示，所以这里每个矩阵内都写成了 $I_x,I_y$，但它们是不同的，表示图像对应窗口元素的梯度值。

就看 $M[0][0]$ 这个元素，发现其实它就是矩阵 $W$ 和矩阵 $I_x$ 在相同窗口位置的卷积值。用一个 $5 imes 5$ 的高斯核作为 $W$ 矩阵。

代码如下：

#!/usr/bin/env python3
# -*- coding:utf8 -*-

import math
import cv2
import numpy as np

def Conv2D(kernel, img, padding, strides=(1, 1)):
    """
    单通道图像卷积
    """
    padImg = np.pad(img, padding, "edge")
    result = []
    for i in range(0, padImg.shape[0] - kernel.shape[0] + 1, strides[0]):
        result.append([])
        for j in range(0, padImg.shape[1] - kernel.shape[1] + 1, strides[1]):
            val = (kernel * padImg[i:i + kernel.shape[0], j:j + kernel.shape[1]]).sum()
            result[-1].append(val)
    return np.array(result)

def Map2Gray(data):
    """
    将数组中的元素映射到[0,255]区间
    """
    maxVal = np.max(data)
    minVal = np.min(data)
    interval = maxVal - minVal
    data = np.floor(255 / interval * (data - minVal))
    return data

def CalImgGrad(f):
    """
    用 5x5 sobel kernel 计算每个像素点的梯度
    """
    sobelx = 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]
    ])
    sobely = 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]
    ])

    fx = Conv2D(sobelx, f, ((2, 2), (2, 2)))
    fy = Conv2D(sobely, f, ((2, 2), (2, 2)))
    fx2 = fx * fx
    fy2 = fy * fy
    fxfy = fx * fy

    return fx2, fy2, fxfy

def CalRvals(fx2, fy2, fxfy):
    """
    计算每个像素对应的M矩阵和R值
    """
    W = 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]
    ])
    k = 0.05

    fx2 = Conv2D(W, fx2, ((2, 2), (2, 2)))
    fy2 = Conv2D(W, fy2, ((2, 2), (2, 2)))
    fxfy = Conv2D(W, fxfy, ((2, 2), (2, 2)))

    result = []
    maxEValue = []   # 存储最大特征值
    minEValue = []   # 存储最小特征值
    for i in range(fx2.shape[0]):
        result.append([])
        maxEValue.append([])
        minEValue.append([])
        for j in range(fx2.shape[1]):
            M = np.array([
                [fx2[i, j], fxfy[i, j]],
                [fxfy[i, j], fy2[i, j]]
            ])
            a = np.linalg.det(M)
            b = np.trace(M)
            """
            M是实对称矩阵，按理必然存在实际特征值，可能由于精度缺失，
            方程 x^2 - bx + a = 0 未必有解，所以 b ** 2 - 4 * a 取了个绝对值，
            这样算出来的值就不对。
            """
            delta = math.sqrt(math.fabs(b ** 2 - 4 * a))
            e1 = 0.5 * (b - delta)
            e2 = 0.5 * (b + delta)
            r = a - k * b
            result[-1].append(r)
            maxEValue[-1].append(max(e1, e2))
            minEValue[-1].append(min(e1, e2))

    g1 = Map2Gray(np.array(maxEValue))         # 最大特征值图
    g2 = Map2Gray(np.array(minEValue))         # 最小特征值图
    g3 = Map2Gray(np.array(result).copy())     # R图

    cv2.imshow('eMax', g1)
    cv2.waitKey(0)
    cv2.imshow('eMin', g2)
    cv2.waitKey(0)
    cv2.imshow('R', g3)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

    return np.array(result), g1, g2, g3

def Sign(img, rVals):
    threshold = 25 * math.fabs(rVals.mean())
    judge = rVals >= threshold
    rectangle = (12, 12)    # 标记矩形大小
    thickness = 2           # 标记矩形边的厚度

    for i in range(img.shape[0]):
        for j in range(img.shape[1]):
            up     = max(i - rectangle[0] // 2, 0)
            bottom = min(i + rectangle[0] // 2, img.shape[0] - 1)
            left   = max(j - rectangle[1] // 2, 0)
            right  = min(j + rectangle[1] // 2, img.shape[0] - 1)
            # 为了使矩形框不发生重叠需要进行非极值抑制
            isLocalExtreme = rVals[i,j] >= rVals[up:bottom + 1, left:right + 1]
            if judge[i, j] and isLocalExtreme.all():
                cv2.rectangle(img, (left, up), (right, bottom), (0, 0, 255), thickness)
    return img

def HarrisCornerDetection(img):
    grayImg = img.mean(axis=-1)  # 变成单通道灰度图像
    fx2, fy2, fxfy = CalImgGrad(grayImg)
    rVals, g1, g2, g3 = CalRvals(fx2, fy2, fxfy)
    finalImg = Sign(img, rVals)
    return finalImg, g1, g2, g3

# =============================================================================

def RecordVideo(fullFileName):
    fps = 25
    size = (640, 480)
    duration = 10  # 修改录制时长，单位s
    frameCount = duration * fps

    cap = cv2.VideoCapture(0, cv2.CAP_DSHOW)  # 打开摄像头
    videoWriter = cv2.VideoWriter(fullFileName, cv2.VideoWriter_fourcc('M', 'J', 'P', 'G'), fps, size)

    print("Begin to record %ss video" % str(duration))

    success, frame = cap.read()
    while success and frameCount > 0:
        videoWriter.write(frame)
        success, frame = cap.read()
        frameCount -= 1

    cap.release()
    print("End of recording")

def TestVideo():
    fullFileName = "video.avi"
    RecordVideo(fullFileName)

    triggerCount = 0
    capVideo = cv2.VideoCapture(fullFileName)
    fps = capVideo.get(cv2.CAP_PROP_FPS)
    frameInterval = int(1000 // fps)   # ms

    print()
    print("Begin to playback video")

    ret, videoFrame = capVideo.read()
    while ret:
        cv2.imshow('image', videoFrame)
        flag = cv2.waitKey(frameInterval)
        if flag == ord(' '):
            triggerCount += 1
            print("[%s]st trigger detection" % str(triggerCount))
            cv2.destroyAllWindows()

            finalImg, g1, g2, g3 = HarrisCornerDetection(videoFrame)
            cv2.imshow('image', finalImg)
            cv2.waitKey(0)

            cv2.imwrite("[%s]-e-max.jpg" % str(triggerCount), g1)
            cv2.imwrite("[%s]-e-min.jpg" % str(triggerCount), g2)
            cv2.imwrite("[%s]-R.jpg" % str(triggerCount), g3)
            cv2.imwrite("[%s]-final.jpg" % str(triggerCount), finalImg)
        ret, videoFrame = capVideo.read()

    capVideo.release()
    cv2.destroyAllWindows()
    print("End of playing")

def TestImage(fullFileName):
    img = cv2.imread(fullFileName)
    finalImg, g1, g2, g3 = HarrisCornerDetection(img)
    cv2.imshow('image', finalImg)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

TestVideo()
# TestImage("2.jpg")