目录
一. 方法介绍
1.1 光流法
1.1.1 光流法理论背景
光流(optical flow)是空间运动物体在观察成像平面上的像素运动的瞬时速度。光流法是利用图像序列中像素在时间域上的变化以及相邻帧之间的相关性来找到上一帧跟当前帧之间存在的对应关系,从而计算出相邻帧之间物体的运动信息的一种方法。通常将二维图像平面特定坐标点上的灰度瞬时变化率定义为光流矢量。
一言以概之:所谓光流就是瞬时速率,在时间间隔很小(比如视频的连续前后两帧之间)时,也等同于目标点的位移
1.1.2 光流的物理意义
一般而言,光流是由于场景中前景目标本身的移动、相机的运动,或者两者的共同运动所产生的。当人的眼睛观察运动物体时,物体的景象在人眼的视网膜上形成一系列连续变化的图像这一系列连续变化的信息不断“流过”视网膜(即图像平面),好像一种光的“流”,故称之为光流。光流表达了图像的变化,由于它包含了目标运动的信息,因此可被观察者用来确定目标的运动情况。
1.1.3 光流法的原理
基本原理:利用连续帧之间的像素亮度变化来估计场景中物体的运动。当物体在图像序列中移动时,它们在图像中的像素值也会随之变化,光流法通过分析这些像素值的变化来推断物体的运动方向和速度。
- 亮度恒定假设(Brightness Constancy Assumption):光流法基于亮度恒定假设,即在短时间内,同一物体的像素亮度保持不变。这意味着在相邻的图像帧中,同一点的像素值应该相似。
- 局部空间一致性假设(Spatial Coherence Assumption):这个假设认为相邻的像素点在运动上是相似的。通过分析相邻像素点之间的关系,可以估计物体的运动。
基于这些假设,光流法的计算步骤通常包括以下几个阶段:
- 特征检测和匹配:在两个相邻的图像帧中,需要检测并匹配一些特征点,这些特征点可以是角点、边缘等。通过这些特征点的匹配,可以建立两帧之间的对应关系。
- 光流计算:对于每个特征点,通过分析其在两帧之间的像素值变化,可以计算出该点的光流向量,即其在图像中的运动方向和速度。
- 光流约束:得到光流向量后,可以通过一些约束条件,如空间一致性和平滑性等,对光流场进行进一步的优化和过滤,以提高估计的准确性。
- 目标跟踪:通过对光流场的分析,可以获取目标的运动信息,从而实现目标的跟踪。
1.2 卡尔曼滤波法
1.2.1 卡尔曼滤波法的理论背景
卡尔曼滤波法的理论背景源自统计学和控制理论。该方法由数学家和控制理论家Rudolf E. Kálmán在20世纪60年代提出,并成为估计与控制领域中的重要工具。
-
马尔可夫过程: 卡尔曼滤波基于马尔可夫过程,即系统状态的未来只依赖于当前状态,而不依赖于过去的状态。这种特性使得卡尔曼滤波能够用递归方式对系统状态进行估计。
-
最小均方估计: 卡尔曼滤波采用了最小均方估计(Minimum Mean Square Estimation,MMSE)的原理。这意味着它通过最小化估计误差的平方来优化状态估计,考虑了系统模型的不确定性和观测数据的噪声。
-
动态系统建模: 卡尔曼滤波要求系统能够用线性动态方程和线性观测方程描述。因为卡尔曼滤波在对系统状态进行预测和更新时,使用了线性系统的性质,使得计算变得相对简单。
-
协方差矩阵: 卡尔曼滤波利用协方差矩阵来表示系统状态和观测数据的不确定性。通过有效地更新和传递协方差信息,卡尔曼滤波能够在不断的观测和预测中提供最优的状态估计。
-
卡尔曼增益: 卡尔曼滤波中的卡尔曼增益是一个关键的概念,它平衡了系统模型和观测数据在状态估计中的权衡。卡尔曼增益的计算考虑了系统不确定性和观测噪声,以提供最优的状态更新。
1.2.2 卡尔曼滤波法的基本原理
-
系统建模: 对待测系统进行建模。这包括定义系统的状态、状态转移方程和观测方程。系统的状态是描述系统行为的变量,状态转移方程描述状态如何随时间演变,观测方程描述系统状态如何被观测到。
-
预测状态: 在每个时间步,利用上一时刻的状态估计和状态转移方程来预测当前时刻的系统状态。这一步产生了预测状态和相应的不确定性。
-
测量更新: 使用传感器或测量设备获得的实际观测数据,通过观测方程将预测状态与观测数据进行比较。通过这一步,系统会更新状态的估计值,考虑了观测数据的影响。
-
卡尔曼增益计算: 卡尔曼增益是一个权重,表示预测状态和观测数据在状态估计中的相对重要性。增益的计算考虑了系统模型的不确定性和观测数据的不确定性。
-
状态更新: 利用卡尔曼增益,将预测状态与观测数据进行融合,得到更精确的状态估计。这一步考虑了系统模型和观测数据的权衡,使估计结果更加可靠。
-
循环迭代: 以上步骤随着时间的推移被循环迭代,以持续更新系统状态的估计值。
原静态视频内容:
二. 代码解决
2.1.1 光流法实现代码
while(cap.isOpened()):
ret, frame = cap.read()
if ret == True:
frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
# 计算光流
p1, _, _ = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)
# 如果找不到特征点,重新选择
if p1 is None:
p0 = cv2.goodFeaturesToTrack(old_gray, mask=None, maxCorners=6, qualityLevel=0.3, minDistance=7, blockSize=7)
continue
frame_centers.append(p1)
# 绘制光流的特征点
for i, (new, old) in enumerate(zip(p1, p0)):
a, b = np.int32(new.ravel()) # Ensure the coordinates are integers
c, d = np.int32(old.ravel())
cv2.circle(frame, (a, b), 5, (0, 0, 255), -1)
cv2.putText(frame, str(i + 1), (a, b), cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 255, 255), 2)
cv2.imshow('Frame', frame)
if cv2.waitKey(25) & 0xFF == ord('q'):
break
# 更新特征点
p0 = p1.reshape(-1, 1, 2)
old_gray = frame_gray.copy()
else:
break
效果展示:
![](https://i-blog.csdnimg.cn/blog_migrate/0999f9e88901279f83be570f2d6d0e1e.png)
可见效果很好,可以比较稳定的标记处光点的几何中心,并且使用红点进行标记,接下来我们计算峰峰值与标准差来查看其每一帧几何中心的变化,我们希望的是其几何中心的每一帧都可以相似或是极端的相等,这种情况对于我们来说是比较友好的,因为其可以稳定的在视频中找出光点并且进行标记,不会随每一帧的细微变化而变化。
# 计算标准差和峰峰值
std_dev_u = np.std(u)
std_dev_v = np.std(v)
peak_to_peak_u = np.max(u) - np.min(u)
peak_to_peak_v = np.max(v) - np.min(v)
print(f'Spot {i+1} U: std_dev={std_dev_u:.2f}, peak_to_peak={peak_to_peak_u}')
print(f'Spot {i+1} V: std_dev={std_dev_v:.2f}, peak_to_peak={peak_to_peak_v}')
运行结果:
Spot 1 U: std_dev=0.34, peak_to_peak=1.39764404296875
Spot 1 V: std_dev=0.28, peak_to_peak=1.35491943359375
Spot 2 U: std_dev=0.23, peak_to_peak=1.013427734375
Spot 2 V: std_dev=0.33, peak_to_peak=1.770263671875
Spot 3 U: std_dev=0.24, peak_to_peak=1.01568603515625
Spot 3 V: std_dev=0.27, peak_to_peak=1.3436279296875
Spot 4 U: std_dev=0.22, peak_to_peak=0.97344970703125
Spot 4 V: std_dev=0.27, peak_to_peak=1.058380126953125
Spot 5 U: std_dev=0.22, peak_to_peak=1.00396728515625
Spot 5 V: std_dev=0.29, peak_to_peak=1.235687255859375
Spot 6 U: std_dev=0.26, peak_to_peak=1.10943603515625
Spot 6 V: std_dev=0.37, peak_to_peak=1.6614990234375
绘图显示:
观察上图或者是运行结果,其峰峰值的变化范围很小,并且其标准差的值也可以稳定在0.5左右,这意味着本方法可以很好的进行该项目的目标追踪,接下来我们将具体情况具体分析,可以在每一帧中加入自己想要的操作进行分析,本方法的具体运用可参考2023年电赛红外光电追踪题目。
2.1.2 卡尔曼滤波法实现代码
kalman_filters = [cv2.KalmanFilter(4, 2) for _ in range(len(p0))]
for i, kf in enumerate(kalman_filters):
kf.measurementMatrix = np.array([[1, 0, 0, 0], [0, 1, 0, 0]], dtype=np.float32)
kf.transitionMatrix = np.array([[1, 0, 1, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]], dtype=np.float32)
kf.processNoiseCov = 1e-3 * np.eye(4, dtype=np.float32)
kf.measurementNoiseCov = 1e-1 * np.eye(2, dtype=np.float32)
while(cap.isOpened()):
ret, frame = cap.read()
if ret == True:
frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
p1, _, _ = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)
if p1 is None:
p0 = cv2.goodFeaturesToTrack(old_gray, mask=None, maxCorners=6, qualityLevel=0.3, minDistance=7, blockSize=7)
continue
frame_centers.append(p1)
for i, (new, old) in enumerate(zip(p1, p0)):
a, b = np.int32(new.ravel())
c, d = np.int32(old.ravel())
cv2.circle(frame, (a, b), 5, (0, 0, 255), -1)
cv2.putText(frame, str(i + 1), (a, b), cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 255, 255), 2)
# Kalman filtering
kalman = kalman_filters[i]
measurement = np.array([[a], [b]], dtype=np.float32)
kalman.correct(measurement)
prediction = kalman.predict()
cv2.circle(frame, (int(prediction[0]), int(prediction[1])), 5, (0, 255, 0), -1)
cv2.imshow('Frame', frame)
if cv2.waitKey(25) & 0xFF == ord('q'):
break
p0 = p1.reshape(-1, 1, 2)
old_gray = frame_gray.copy()
else:
break
效果展示:
该效果与上述光流法效果类似,但是该方法会有一个缓存前夕,在0.5s内会逐渐移动到目标光点集几何中心上,但最终结果是表现良好的。
对每一帧进行进行分析,并且绘制峰峰值与标准差图像:
for i in range(6):
tracked_points = [frame[i].ravel() for frame in frame_centers if len(frame) > i]
if not tracked_points:
continue
u, v = zip(*tracked_points)
if len(u) <= 1:
continue
plt.plot(u, label=f'Spot {i + 1} U')
plt.plot(v, label=f'Spot {i + 1} V')
std_dev_u = np.std(u)
std_dev_v = np.std(v)
peak_to_peak_u = np.max(u) - np.min(u)
peak_to_peak_v = np.max(v) - np.min(v)
print(f'Spot {i + 1} U: std_dev={std_dev_u:.2f}, peak_to_peak={peak_to_peak_u}')
print(f'Spot {i + 1} V: std_dev={std_dev_v:.2f}, peak_to_peak={peak_to_peak_v}')
plt.legend()
plt.show()
运行结果:
Spot 1 U: std_dev=0.32, peak_to_peak=1.39764404296875
Spot 1 V: std_dev=0.27, peak_to_peak=1.2662353515625
Spot 2 U: std_dev=0.24, peak_to_peak=1.013427734375
Spot 2 V: std_dev=0.34, peak_to_peak=1.770263671875
Spot 3 U: std_dev=0.23, peak_to_peak=1.01568603515625
Spot 3 V: std_dev=0.26, peak_to_peak=1.3436279296875
Spot 4 U: std_dev=0.22, peak_to_peak=0.97344970703125
Spot 4 V: std_dev=0.26, peak_to_peak=1.058380126953125
Spot 5 U: std_dev=0.22, peak_to_peak=1.00396728515625
Spot 5 V: std_dev=0.27, peak_to_peak=1.175537109375
Spot 6 U: std_dev=0.23, peak_to_peak=0.97100830078125
Spot 6 V: std_dev=0.37, peak_to_peak=1.6614990234375
绘图显示:
同样的,观察上图或者是运行结果,其峰峰值的变化范围很小,并且其标准差的值也可以稳定在0.4左右,这意味着本方法可以很好的进行该项目的目标追踪.
三. 完整代码
3.1 光流法
import cv2
import numpy as np
import matplotlib.pyplot as plt
import statistics
# 定义一个函数来获取白色的掩码
def get_white_mask(image):
lower_white = np.array([0, 0, 50])
upper_white = np.array([331, 10, 100])
hsv = cv2.cvtColor(image, cv2.COLOR_BGR2HSV)
mask = cv2.inRange(hsv, lower_white, upper_white)
return mask
cap = cv2.VideoCapture('静态视频.avi')
# 创建Lucas-Kanade光流法参数
lk_params = dict(winSize=(15, 15), maxLevel=2, criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))
# 选择初始特征点
ret, old_frame = cap.read()
old_gray = cv2.cvtColor(old_frame, cv2.COLOR_BGR2GRAY)
p0 = cv2.goodFeaturesToTrack(old_gray, mask=None, maxCorners=6, qualityLevel=0.3, minDistance=7, blockSize=7)
frame_centers = [] # 用于存储每一帧的光流法计算得到的光斑中心
while(cap.isOpened()):
ret, frame = cap.read()
if ret == True:
frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
# 计算光流
p1, _, _ = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)
# 如果找不到特征点,重新选择
if p1 is None:
p0 = cv2.goodFeaturesToTrack(old_gray, mask=None, maxCorners=6, qualityLevel=0.3, minDistance=7, blockSize=7)
continue
frame_centers.append(p1)
# 绘制光流的特征点
for i, (new, old) in enumerate(zip(p1, p0)):
a, b = np.int32(new.ravel()) # Ensure the coordinates are integers
c, d = np.int32(old.ravel())
cv2.circle(frame, (a, b), 5, (0, 0, 255), -1)
cv2.putText(frame, str(i + 1), (a, b), cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 255, 255), 2)
cv2.imshow('Frame', frame)
if cv2.waitKey(25) & 0xFF == ord('q'):
break
# 更新特征点
p0 = p1.reshape(-1, 1, 2)
old_gray = frame_gray.copy()
else:
break
cap.release()
cv2.destroyAllWindows()
# 绘制每一帧的光流法计算得到的光斑中心
for i in range(6): # 假设最多跟踪6个光斑
tracked_points = [frame[i].ravel() for frame in frame_centers if len(frame) > i]
if not tracked_points:
continue # 如果没有追踪到该光斑的点,则跳过
u, v = zip(*tracked_points)
if len(u) <= 1:
# 如果只有一个点,标准差未定义,跳过
continue
plt.plot(u, label=f'Spot {i+1} U')
plt.plot(v, label=f'Spot {i+1} V')
# 计算标准差和峰峰值
std_dev_u = np.std(u)
std_dev_v = np.std(v)
peak_to_peak_u = np.max(u) - np.min(u)
peak_to_peak_v = np.max(v) - np.min(v)
print(f'Spot {i+1} U: std_dev={std_dev_u:.2f}, peak_to_peak={peak_to_peak_u}')
print(f'Spot {i+1} V: std_dev={std_dev_v:.2f}, peak_to_peak={peak_to_peak_v}')
plt.legend()
plt.show()
3.2 卡尔曼滤波法
import cv2
import numpy as np
import matplotlib.pyplot as plt
def get_white_mask(image):
lower_white = np.array([0, 0, 50])
upper_white = np.array([331, 10, 100])
hsv = cv2.cvtColor(image, cv2.COLOR_BGR2HSV)
mask = cv2.inRange(hsv, lower_white, upper_white)
return mask
cap = cv2.VideoCapture('静态视频.avi')
lk_params = dict(winSize=(15, 15), maxLevel=2, criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0.03))
ret, old_frame = cap.read()
old_gray = cv2.cvtColor(old_frame, cv2.COLOR_BGR2GRAY)
p0 = cv2.goodFeaturesToTrack(old_gray, mask=None, maxCorners=6, qualityLevel=0.3, minDistance=7, blockSize=7)
frame_centers = []
kalman_filters = [cv2.KalmanFilter(4, 2) for _ in range(len(p0))]
for i, kf in enumerate(kalman_filters):
kf.measurementMatrix = np.array([[1, 0, 0, 0], [0, 1, 0, 0]], dtype=np.float32)
kf.transitionMatrix = np.array([[1, 0, 1, 0], [0, 1, 0, 1], [0, 0, 1, 0], [0, 0, 0, 1]], dtype=np.float32)
kf.processNoiseCov = 1e-3 * np.eye(4, dtype=np.float32)
kf.measurementNoiseCov = 1e-1 * np.eye(2, dtype=np.float32)
while(cap.isOpened()):
ret, frame = cap.read()
if ret == True:
frame_gray = cv2.cvtColor(frame, cv2.COLOR_BGR2GRAY)
p1, _, _ = cv2.calcOpticalFlowPyrLK(old_gray, frame_gray, p0, None, **lk_params)
if p1 is None:
p0 = cv2.goodFeaturesToTrack(old_gray, mask=None, maxCorners=6, qualityLevel=0.3, minDistance=7, blockSize=7)
continue
frame_centers.append(p1)
for i, (new, old) in enumerate(zip(p1, p0)):
a, b = np.int32(new.ravel())
c, d = np.int32(old.ravel())
cv2.circle(frame, (a, b), 5, (0, 0, 255), -1)
cv2.putText(frame, str(i + 1), (a, b), cv2.FONT_HERSHEY_SIMPLEX, 1, (255, 255, 255), 2)
kalman = kalman_filters[i]
measurement = np.array([[a], [b]], dtype=np.float32)
kalman.correct(measurement)
prediction = kalman.predict()
cv2.circle(frame, (int(prediction[0]), int(prediction[1])), 5, (0, 255, 0), -1)
cv2.imshow('Frame', frame)
if cv2.waitKey(25) & 0xFF == ord('q'):
break
p0 = p1.reshape(-1, 1, 2)
old_gray = frame_gray.copy()
else:
break
cap.release()
cv2.destroyAllWindows()
for i in range(6):
tracked_points = [frame[i].ravel() for frame in frame_centers if len(frame) > i]
if not tracked_points:
continue
u, v = zip(*tracked_points)
if len(u) <= 1:
continue
plt.plot(u, label=f'Spot {i + 1} U')
plt.plot(v, label=f'Spot {i + 1} V')
std_dev_u = np.std(u)
std_dev_v = np.std(v)
peak_to_peak_u = np.max(u) - np.min(u)
peak_to_peak_v = np.max(v) - np.min(v)
print(f'Spot {i + 1} U: std_dev={std_dev_u:.2f}, peak_to_peak={peak_to_peak_u}')
print(f'Spot {i + 1} V: std_dev={std_dev_v:.2f}, peak_to_peak={peak_to_peak_v}')
plt.legend()
plt.show()