Get the original image from google drive. This movie was already stabilized in iMovie. This removed the high frequency shaking, but left the slow drift motion. This notebook detects the sun in every frame, shifts each frame such that the sun is always in the same position, and writes the motion-corrected frames to disk, and finally stitches them together into a movie.
There are a lot of ways to accomplish this task. This algorithm is very simple -- matching a sun template -- but not very flexible. An alternative is to directly match consecutive frames, i.e., register the images.
In [1]:
### Get the data
import os
import requests
# Download original movie file
src_filename = 'eclipse.mp4'
if not os.path.exists(src_filename):
!wget "https://docs.google.com/uc?export=download&id=0B-ayIMQgiR1dUGx3MEVFRzVxYzQ" -O $src_filename
else:
print ('Source movie file already downloaded.')
# Prepare temporary output directory
output_frame_dir = 'out_frames'
if not os.path.exists(output_frame_dir):
os.makedirs(output_frame_dir)
In [2]:
import numpy as np
import skvideo.io
videodata = skvideo.io.vread(src_filename)
gray_frames = videodata.sum(axis=-1)
print(videodata.shape, gray_frames.shape)
In [3]:
### Register each frame to a reference frame
from skimage.feature import register_translation
shifts, errs = [], []
ref_frame = gray_frames[200]
for frame_num, frame in enumerate(gray_frames):
if frame_num % 50 == 0:
print('Frame %i of %i' % (frame_num, len(videodata)))
shift, err, _ = register_translation(ref_frame, frame)
shifts.append(shift)
errs.append(err)
shifts = np.round(np.array(shifts)).astype(int)
Note that with a frame/template of this size, naive convolution is too slow, and fftconvolve is necessary. It would be much more efficient to FFT the template only once, compute the spatial FFT of the video data, and then just multiply them. This is fast enough for a demo!
In [4]:
### Smooth an plot the "trajectory" of the sun
from math import factorial
import matplotlib.pyplot as plt
%matplotlib inline
# from here: http://scipy.github.io/old-wiki/pages/Cookbook/SavitzkyGolay
def savitzky_golay(y, window_size, order, deriv=0, rate=1):
try:
window_size = np.abs(np.int(window_size))
order = np.abs(np.int(order))
except ValueError:
raise ValueError("window_size and order have to be of type int")
if window_size % 2 != 1 or window_size < 1:
raise TypeError("window_size size must be a positive odd number")
if window_size < order + 2:
raise TypeError("window_size is too small for the polynomials order")
order_range = range(order+1)
half_window = (window_size -1) // 2
# precompute coefficients
b = np.mat([[k**i for i in order_range] for k in range(-half_window, half_window+1)])
m = np.linalg.pinv(b).A[deriv] * rate**deriv * factorial(deriv)
# pad the signal at the extremes with
# values taken from the signal itself
firstvals = y[0] - np.abs( y[1:half_window+1][::-1] - y[0] )
lastvals = y[-1] + np.abs(y[-half_window-1:-1][::-1] - y[-1])
y = np.concatenate((firstvals, y, lastvals))
return np.convolve( m[::-1], y, mode='valid')
t = np.arange(len(shifts))
win_size = 11
degree = 5
row_shifts_smooth = savitzky_golay(shifts[:,0], win_size, degree)
col_shifts_smooth = savitzky_golay(shifts[:,1], win_size, degree)
shifts_smooth = np.vstack([row_shifts_smooth, col_shifts_smooth]).T
plt.figure(figsize=(9,6))
plt.plot(shifts_smooth[:,1], shifts_smooth[:,0], 'b-')
plt.plot(shifts[:,1], shifts[:,0], 'rx')
Out[4]:
In [5]:
import numpy as np
from skimage.util import crop
crop_size = int(round(np.abs(shifts).max()))
shifts_smooth = np.round(shifts_smooth).astype(int)
for frame_num, (frame, shift) in enumerate(zip(videodata, shifts_smooth)):
if frame_num % 50 == 0:
print('Frame %i of %i' % (frame_num, len(videodata)))
# Shift the frame
shifted_frame = np.roll(frame, shift=(shift[0], shift[1], 0), axis=(0,1,2))
# Zero-fill the wrapped values
out_frame = shifted_frame.copy()
if shift[0] != 0:
if shift[0] > 0:
out_frame[:shift[0]] = 0
else:
out_frame[shift[0]:] = 0
if shift[1] != 0:
if shift[1] > 0:
out_frame[:, :shift[1]] = 0
else:
out_frame[:, shift[1]:] = 0
# Save motion-compensated frame
out_fn = os.path.join(output_frame_dir, 'frame%03d.png' % frame_num)
plt.imsave(out_fn, out_frame)
In [6]:
### Convert frames to a mp4 movie
!ffmpeg -r 30 -i $output_frame_dir/frame%03d.png -vcodec libx264 -crf 25 -pix_fmt yuv420p eclipse_corrected.mp4
In [7]:
ring_frame = videodata[130]
plt.figure(figsize=(16,10))
plt.imshow(ring_frame)
plt.imsave('diamond_ring.jpg', ring_frame)
In [ ]: