Harris Corner Detection
这个博客写的不详细，有兴趣的话可以看视频：https://www.bilibili.com/video/BV1Hv411e7CW。
使用一个固定窗口在图像上进行任意方向上的滑动，比较滑动前与滑动后两种情况，窗口中的像素灰度变化程度，如果存在任意方向上的滑动，
都有着较大灰度变化，那么我们就可以认为该窗口中存在角点。举个例子，下图中 $A$ 是角点，$B$ 不是角点，如果用一个小窗口(蓝色)在点 $A$
处滑动，无论朝哪个方向走，窗口内的像素值都会发生明显的变化，比如窗口往左，那白色像素就越多，红色越少；但是对于 $B$ 点，无论朝哪
个方向滑动，窗口内的像素基本不变。
进一步地
基本数学公式如下：
$$E(u,v) = \sum_{(x,y) \in W}^{}W(x,y)\left [ I(x+u,y+v) - I(x,y) \right ]^{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) \right ]^{2}
= \sum_{(x,y) \in W}^{}W(x,y) \left (u^2I_x^2 + 2uvI_xIy + v^2I_y^2 \right ) \\
= \sum_{(x,y) \in W}^{}W(x,y)\begin{bmatrix}
u & v
\end{bmatrix}\begin{bmatrix}
I_x^2 & IxIy\\
IxIy & I_y^2
\end{bmatrix}\begin{bmatrix}
u \\
v
\end{bmatrix} \\
= \begin{bmatrix}
u & v
\end{bmatrix}\left (\sum_{(x,y) \in W}^{}W(x,y)\begin{bmatrix}
I_x^2 & IxIy\\
IxIy & I_y^2
\end{bmatrix} \right ) \begin{bmatrix}
u \\
v
\end{bmatrix} \\
=
\begin{bmatrix}
u & v
\end{bmatrix}M \begin{bmatrix}
u \\
v
\end{bmatrix}$$
每个实对称矩阵累加之后还是一个实对称矩阵，所以得到的 $M$ 是一个实对称矩阵。进一步有：
$$E(u,v) \approx \begin{bmatrix}
u & v
\end{bmatrix}P\begin{bmatrix}
\lambda_1 & \\
& \lambda_2
\end{bmatrix}P^T \begin{bmatrix}
u \\
v
\end{bmatrix} \\
= \begin{bmatrix}
u & v
\end{bmatrix}P\begin{bmatrix}
\lambda_1 & \\
& \lambda_2
\end{bmatrix} \left (\begin{bmatrix}
u & v
\end{bmatrix}P \right )^T \\
= \begin{bmatrix}
u^{'} & v^{'}
\end{bmatrix}\begin{bmatrix}
\lambda_1 & \\
& \lambda_2
\end{bmatrix} \begin{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_1\lambda_2 - \left ( \lambda_1 + \lambda_2 \right )$$
代码实现过程主要在于求得矩阵 $M$，即
$$M = \left (\sum_{(x,y) \in W}^{}w(x,y)\begin{bmatrix} I_x^2 & IxIy\\ IxIy & I_y^2 \end{bmatrix} \right )$$
假设窗口大小为 $5 \times 5$，即
$$W = \begin{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 = \begin{bmatrix}
w_{11}I_x^2 & w_{11}I_xI_y \\
w_{11}I_xI_y & w_{11}I_y^2
\end{bmatrix} + \begin{bmatrix}
w_{12}I_x^2 & w_{12}I_xI_y \\
w_{12}I_xI_y & w_{12}I_y^2
\end{bmatrix} + \cdots + \begin{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 \times 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")
