用opencv实现视频稳定防抖video stabilization（来自learnopencv）

代码来自https://www.learnopencv.com/video-stabilization-using-point-feature-matching-in-opencv/

其中涉及视频光流、图像特征点提取、仿射变换（二维移动、旋转，缩放）、投影变换（三维投影）、用卷积实现移动平均等，是挺好的opencv综合练习材料。

原文的参考资料 https://abhitronix.github.io/2018/11/30/humanoid-AEAM-3/ 有一些视频解释，有助于理解代码。

import numpy as np
import cv2

cap = cv2.VideoCapture('video.mp4');
n_frames = int(cap.get(cv2.CAP_PROP_FRAME_COUNT)) 
w = int(cap.get(cv2.CAP_PROP_FRAME_WIDTH)) 
h = int(cap.get(cv2.CAP_PROP_FRAME_HEIGHT))

# https://www.fourcc.org/codecs.php    https://www.cnblogs.com/zhangzhihui/p/12503152.html 关于用frame构造新视频写入时压缩格式和视频大小的说明
fourcc = cv2.VideoWriter_fourcc('X','V','I','D')
fps = cap.get(cv2.CAP_PROP_FPS);

out = cv2.VideoWriter('video_out.avi', fourcc, fps, (2*w, h))
# Read first frame
_, prev = cap.read()
# Convert frame to grayscale
prev_gray = cv2.cvtColor(prev, cv2.COLOR_BGR2GRAY) 


# Pre-define transformation-store array
transforms = np.zeros((n_frames-1, 3), np.float32) 

for i in range(n_frames-2):
  # Detect feature points in previous frame
  prev_pts = cv2.goodFeaturesToTrack(prev_gray,
                                     maxCorners=200,
                                     qualityLevel=0.01,
                                     minDistance=30,
                                     blockSize=3)
   
  # Read next frame
  success, curr = cap.read() 
  if not success: 
    break 

  # Convert to grayscale
  curr_gray = cv2.cvtColor(curr, cv2.COLOR_BGR2GRAY) 

  # Calculate optical flow (i.e. track feature points)
  curr_pts, status, err = cv2.calcOpticalFlowPyrLK(prev_gray, curr_gray, prev_pts, None) 

  # Sanity check
  assert prev_pts.shape == curr_pts.shape 

  # Filter only valid points
  idx = np.where(status==1)[0]
  prev_pts = prev_pts[idx]
  curr_pts = curr_pts[idx]

  #Find transformation matrix
  m, _ = cv2.estimateAffinePartial2D(prev_pts, curr_pts)
   
  # Extract traslation
  dx = m[0,2]
  dy = m[1,2]

  # Extract rotation angle
  da = np.arctan2(m[1,0], m[0,0])
   
  # Store transformation
  transforms[i] = [dx,dy,da]
   
  # Move to next frame
  prev_gray = curr_gray

  print("Frame: " + str(i) +  "/" + str(n_frames) + " -  Tracked points : " + str(len(prev_pts)))

# https://xiaosongshine.blog.csdn.net/article/details/90258331
# https://www.cnblogs.com/xiaosongshine/p/10874644.html
# https://www.cnblogs.com/yuxiangyang/p/11161789.html 有卷积的计算方法
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 

SMOOTHING_RADIUS = 50
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

# 通过cumsum累计将frame每次变动的相对偏移建立同样的基准。每个trajectory值都变为相对原点的变动。
# Compute trajectory using cumulative sum of transformations
trajectory = np.cumsum(transforms, axis=0) 

# Smooth trajectory using moving average filter
smoothed_trajectory = smooth(trajectory); 

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

# 还原到相对量用于相邻两图的仿射变换
# Calculate newer transformation array
transforms_smooth = transforms + difference

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

# 通过仿射变换来放大图像，遮住因图片移动、旋转、缩放导致不能覆盖住frame的长宽的问题
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

# 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]
  
  # https://www.cnblogs.com/silence-cho/p/10926248.html 解释仿射变换和投影的计算，下面的cos/sin/...都是按照固定公式进行的计算
  # 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)

cap.release()
out.release()
cv2.destroyAllWindows()