声音可能不太好听，是用 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) \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")