// Pre-define transformation-store array
vector <TransformParam> transforms;

//
Mat last_T;

for(int i = 1; i < n_frames-1; i++)
{
// Vector from previous and current feature points
vector <Point2f> prev_pts, curr_pts;

// Detect features in previous frame
goodFeaturesToTrack(prev_gray, prev_pts, 200, 0.01, 30);

// Read next frame
bool success = cap.read(curr);
if(!success) break;

// Convert to grayscale
cvtColor(curr, curr_gray, COLOR_BGR2GRAY);

// Calculate optical flow (i.e. track feature points)
vector <uchar> status;
vector <float> err;
calcOpticalFlowPyrLK(prev_gray, curr_gray, prev_pts, curr_pts, status, err);

// Filter only valid points
auto prev_it = prev_pts.begin();
auto curr_it = curr_pts.begin();
for(size_t k = 0; k < status.size(); k++)
{
if(status[k])
{
prev_it++;
curr_it++;
}
else
{
prev_it = prev_pts.erase(prev_it);
curr_it = curr_pts.erase(curr_it);
}
}

// Find transformation matrix
Mat T = estimateRigidTransform(prev_pts, curr_pts, false);

// In rare cases no transform is found.
// We'll just use the last known good transform.
if(T.data == NULL) last_T.copyTo(T);
T.copyTo(last_T);

// Extract traslation
double dx = T.at<double>(0,2);
double dy = T.at<double>(1,2);

// Extract rotation angle
double da = atan2(T.at<double>(1,0), T.at<double>(0,0));

// Store transformation
transforms.push_back(TransformParam(dx, dy, da));

// Move to next frame
curr_gray.copyTo(prev_gray);

cout << "Frame: " << i << "/" << n_frames << " -  Tracked points : " << prev_pts.size() << endl;
}

在前面的步骤中，我们估计帧之间的运动并将它们存储在一个数组中。我们现在需要通过叠加上一步估计的微分运动来找到运动轨迹。

步骤4.1:轨迹计算

在这一步，我们将增加运动之间的帧来计算轨迹。我们的最终目标是平滑这条轨迹。

Python 在Python中，可以很容易地使用numpy中的cumsum(累计和)来实现。

# Compute trajectory using cumulative sum of transformations
trajectory = np.cumsum(transforms, axis=0

C++

在c++中，我们定义了一个名为Trajectory的结构体来存储转换参数的累积和。

struct Trajectory
{
Trajectory() {}
Trajectory(double _x, double _y, double _a) {
x = _x;
y = _y;
a = _a;
}

double x;
double y;
double a; // angle
};

vector<Trajectory> cumsum(vector<TransformParam> &transforms)
{
vector <Trajectory> trajectory; // trajectory at all frames
// Accumulated frame to frame transform
double a = 0;
double x = 0;
double y = 0;

for(size_t i=0; i < transforms.size(); i++)
{
x += transforms[i].dx;
y += transforms[i].dy;
a += transforms[i].da;

trajectory.push_back(Trajectory(x,y,a));

}

return trajectory;
}

在上一步中，我们计算了运动轨迹。所以我们有三条曲线来显示运动(x, y，和角度)如何随时间变化。

在这一步，我们将展示如何平滑这三条曲线。

平滑任何曲线最简单的方法是使用移动平均滤波器（moving average filter）。顾名思义，移动平均过滤器将函数在某一点上的值替换为由窗口定义的其相邻函数的平均值。让我们看一个例子。

假设我们在数组c中存储了一条曲线，那么曲线上的点是c[0]…c[n-1]。设f是我们通过宽度为5的移动平均滤波器过滤c得到的平滑曲线。

该曲线的k^{th}元素使用

如您所见，平滑曲线的值是噪声曲线在一个小窗口上的平均值。下图显示了左边的噪点曲线的例子，使用右边的尺度为5 滤波器进行平滑。

Python

在Python实现中，我们定义了一个移动平均滤波器，它接受任何曲线(即1-D的数字)作为输入，并返回曲线的平滑版本。

def movingAverage(curve, radius):
window_size = 2 * radius + 1
# Define the filter
f = np.ones(window_size)/window_size
# Add padding to the boundaries
curve_pad = np.lib.pad(curve, (radius, radius), 'edge')
# Apply convolution
curve_smoothed = np.convolve(curve_pad, f, mode='same')
# Remove padding
curve_smoothed = curve_smoothed[radius:-radius]
# return smoothed curve
return curve_smoothed

def smooth(trajectory):
smoothed_trajectory = np.copy(trajectory)
# Filter the x, y and angle curves
for i in range(3):
smoothed_trajectory[:,i] = movingAverage(trajectory[:,i], radius=SMOOTHING_RADIUS)

return smoothed_trajectory

# Compute trajectory using cumulative sum of transformations
trajectory = np.cumsum(transforms, axis=0)

在c++版本中，我们定义了一个名为smooth的函数，用于计算平滑移动平均轨迹。

vector <Trajectory> smooth(vector <Trajectory>& trajectory, int radius)
{
vector <Trajectory> smoothed_trajectory;
for(size_t i=0; i < trajectory.size(); i++) {
double sum_x = 0;
double sum_y = 0;
double sum_a = 0;
int count = 0;

for(int j=-radius; j <= radius; j++) { if(i+j >= 0 && i+j < trajectory.size()) {
sum_x += trajectory[i+j].x;
sum_y += trajectory[i+j].y;
sum_a += trajectory[i+j].a;

count++;
}
}

double avg_a = sum_a / count;
double avg_x = sum_x / count;
double avg_y = sum_y / count;

smoothed_trajectory.push_back(Trajectory(avg_x, avg_y, avg_a));
}

return smoothed_trajectory;
}

// Smooth trajectory using moving average filter
vector <Trajectory> smoothed_trajectory = smooth(trajectory, SMOOTHING_RADIUS);

这是通过找到平滑轨迹和原始轨迹之间的差异，并将这些差异加回到原始的变换中来完成的。

Python

# Calculate difference in smoothed_trajectory and trajectory
difference = smoothed_trajectory - trajectory

# Calculate newer transformation array
transforms_smooth = transforms + difference

C++

vector <TransformParam> transforms_smooth;
for(size_t i=0; i < transforms.size(); i++)
{
// Calculate difference in smoothed_trajectory and trajectory
double diff_x = smoothed_trajectory[i].x - trajectory[i].x;
double diff_y = smoothed_trajectory[i].y - trajectory[i].y;
double diff_a = smoothed_trajectory[i].a - trajectory[i].a;
// Calculate newer transformation array
double dx = transforms[i].dx + diff_x;
double dy = transforms[i].dy + diff_y;
double da = transforms[i].da + diff_a;
transforms_smooth.push_back(TransformParam(dx, dy, da));
}

第五步:将平滑的摄像机运动应用到帧中

差不多做完了。现在我们所需要做的就是循环帧并应用我们刚刚计算的变换。

如果我们有一个指定为(x, y, \theta),的运动，对应的变换矩阵是

请阅读代码中的注释以进行后续操作。

Python

# Reset stream to first frame
cap.set(cv2.CAP_PROP_POS_FRAMES, 0)

# Write n_frames-1 transformed frames
for i in range(n_frames-2):
# Read next frame
success, frame = cap.read()
if not success:
break

# Extract transformations from the new transformation array
dx = transforms_smooth[i,0]
dy = transforms_smooth[i,1]
da = transforms_smooth[i,2]

# Reconstruct transformation matrix accordingly to new values
m = np.zeros((2,3), np.float32)
m[0,0] = np.cos(da)
m[0,1] = -np.sin(da)
m[1,0] = np.sin(da)
m[1,1] = np.cos(da)
m[0,2] = dx
m[1,2] = dy

# Apply affine wrapping to the given frame
frame_stabilized = cv2.warpAffine(frame, m, (w,h))

# Fix border artifacts
frame_stabilized = fixBorder(frame_stabilized)

# Write the frame to the file
frame_out = cv2.hconcat([frame, frame_stabilized])

# If the image is too big, resize it.
if(frame_out.shape[1] > 1920):
frame_out = cv2.resize(frame_out, (frame_out.shape[1]/2, frame_out.shape[0]/2));

cv2.imshow("Before and After", frame_out)
cv2.waitKey(10)
out.write(frame_out)

C++

cap.set(CV_CAP_PROP_POS_FRAMES, 1);
Mat T(2,3,CV_64F);
Mat frame, frame_stabilized, frame_out;

for( int i = 0; i < n_frames-1; i++) { bool success = cap.read(frame); if(!success) break; // Extract transform from translation and rotation angle. transforms_smooth[i].getTransform(T); // Apply affine wrapping to the given frame warpAffine(frame, frame_stabilized, T, frame.size()); // Scale image to remove black border artifact fixBorder(frame_stabilized); // Now draw the original and stabilised side by side for coolness hconcat(frame, frame_stabilized, frame_out); // If the image is too big, resize it. if(frame_out.cols > 1920)
{
resize(frame_out, frame_out, Size(frame_out.cols/2, frame_out.rows/2));
}

imshow("Before and After", frame_out);
out.write(frame_out);
waitKey(10);
}

步骤5.1:修复边界伪影
当我们稳定一个视频，我们可能会看到一些黑色的边界伪影。这是意料之中的，因为为了稳定视频，帧可能不得不缩小大小。

我们可以通过将视频的中心缩小一小部分(例如4%)来缓解这个问题。

下面的fixBorder函数显示了实现。我们使用getRotationMatrix2D，因为它在不移动图像中心的情况下缩放和旋转图像。我们所需要做的就是调用这个函数时，旋转为0，缩放为1.04(也就是提升4%)。

Python

def fixBorder(frame):
s = frame.shape
# Scale the image 4% without moving the center
T = cv2.getRotationMatrix2D((s[1]/2, s[0]/2), 0, 1.04)
frame = cv2.warpAffine(frame, T, (s[1], s[0]))
return frame

void fixBorder(Mat &frame_stabilized)
{
Mat T = getRotationvoid fixBorder(Mat &frame_stabilized)
{
Mat T = getRotationMatrix2D(Point2f(frame_stabilized.cols/2, frame_stabilized.rows/2), 0, 1.04);
warpAffine(frame_stabilized, frame_stabilized, T, frame_stabilized.size());
}Matrix2D(Point2f(frame_stabilized.cols/2, frame_stabilized.rows/2), 0, 1.04);
warpAffine(frame_stabilized, frame_stabilized, T, frame_stabilized.size());
}

我们分享的视频防抖代码的结果如上所示。我们的目标是显著减少运动，但不是完全消除它。

我们留给读者去思考如何修改代码来完全消除帧之间的移动。如果你试图消除所有的相机运动，会有什么副作用?

目前的方法只适用于固定长度的视频，而不适用于实时feed。我们不得不对这个方法进行大量修改，以获得实时视频输出，这超出了本文的范围，但这是可以实现的，更多的信息可以在这里找到。

https://abhitronix.github.io/2018/11/30/humanoid-AEAM-3/

优点和缺点

优点

这种方法对低频运动(较慢的振动)具有良好的稳定性。这种方法内存消耗低，因此非常适合嵌入式设备(如树莓派)。这种方法对视频缩放抖动有很好的效果。

缺点

这种方法对高频扰动的抵抗效果很差。如果有一个严重的运动模糊，特征跟踪将失败，结果将不是最佳的。这种方法也不适用于滚动快门失真。

References:

1. Example video and Code reference from Nghia Ho’s post

   http://nghiaho.com/uploads/code/videostab.cpp

2. Various References, data, and image from my website

   https://abhitronix.github.io/

3. https://www.learnopencv.com/video-stabilization-using-point-feature-matching-in-opencv/

1. 关注【3D视觉工坊】公众号

2. 在【3D视觉工坊】公众号后台回复 【视频防抖】 即可获取

本文仅做学术分享，如有侵权，请联系删文。下载1在「3D视觉工坊」公众号后台回复：3D视觉，即可下载 3D视觉相关资料干货，涉及相机标定、三维重建、立体视觉、SLAM、深度学习、点云后处理、多视图几何等方向。下载2在「3D视觉工坊」公众号后台回复：3D视觉github资源汇总，即可下载包括结构光、标定源码、缺陷检测源码、深度估计与深度补全源码、点云处理相关源码、立体匹配源码、单目、双目3D检测、基于点云的3D检测、6D姿态估计源码汇总等。下载3在「3D视觉工坊」公众号后台回复：相机标定，即可下载独家相机标定学习课件与视频网址；后台回复：立体匹配，即可下载独家立体匹配学习课件与视频网址。
重磅！3DCVer-学术论文写作投稿 交流群已成立
扫码添加小助手微信，可申请加入3D视觉工坊-学术论文写作与投稿 微信交流群，旨在交流顶会、顶刊、SCI、EI等写作与投稿事宜。

同时也可申请加入我们的细分方向交流群，目前主要有3D视觉、CV&深度学习、SLAM、三维重建、点云后处理、自动驾驶、CV入门、三维测量、VR/AR、3D人脸识别、医疗影像、缺陷检测、行人重识别、目标跟踪、视觉产品落地、视觉竞赛、车牌识别、硬件选型、学术交流、求职交流等微信群。
一定要备注：研究方向+学校/公司+昵称，例如：”3D视觉 + 上海交大 + 静静“。请按照格式备注，可快速被通过且邀请进群。原创投稿也请联系。
▲长按加微信群或投稿▲长按关注公众号
3D视觉从入门到精通知识星球：针对3D视觉领域的知识点汇总、入门进阶学习路线、最新paper分享、疑问解答四个方面进行深耕，更有各类大厂的算法工程人员进行技术指导。与此同时，星球将联合知名企业发布3D视觉相关算法开发岗位以及项目对接信息，打造成集技术与就业为一体的铁杆粉丝聚集区，近2000星球成员为创造更好的AI世界共同进步，知识星球入口：
学习3D视觉核心技术，扫描查看介绍，3天内无条件退款

 圈里有高质量教程资料、可答疑解惑、助你高效解决问题
觉得有用，麻烦给个赞和在看~