See https://www.reddit.com/r/spacex/comments/528b1k/looped_video_70_frames_prefireball_with_a_couple/d7i45o9.

sx.m

% 1. download files from http://alumni.soe.ucsc.edu/%7Exzhu/doc/turbulence.html to current folder
%    and http://tmp.hejnoah.com/frames_Sothh.zip to subfolder "frames_Sothh"
%
% 2. run
%    $ youtube-dl -o sx.mp4 https://www.youtube.com/watch?v=_BgJEXQkjNQ
%    to download the video and
%    $ ffmpeg -ss 00:01:05 -i sx.mp4 -filter:v fps=fps=60/1 frames/ffmpeg_%3d.png
%    to extract the frames - quit this once it reaches frame 450 or so
%
% 3. then run this script
%
% 4 to generate video, run
%   $ cd frames_allthree
%   $ mogrify -resize 1280x720 *.png
%   $ ffmpeg -framerate 10 -i frame_%03d.png -c:v libx264 -vf "scale=trunc(iw/2)*2:trunc(ih/2)*2" -r 30 -pix_fmt yuv420p frames.mp4

% TODO possible improvement: look at color channels separately

if exist('frames_reg','var') == 0
    
    % read images (if not already read)
    if exist('imgs','var') == 0
        imgPath = 'frames/';
        imgType = '*.png';
        images  = dir([imgPath imgType]);
        for idx = 1:length(images)
            images(idx).name
            img = imread([imgPath images(idx).name]);
            img = img(300:460,480:610,:);  % area of interest
            img = rgb2gray(img);
            imgs(:,:,idx) = img;
        end
    end

    % setup variables
    frames = imgs;
    Ref = mean(imgs, 3);
    imshow(Ref);

    % run the turbulence removal, note that this takes ~30s per frame
    % TODO improve settings (see Example.m)
    Example
end

% add additional frame to skip handling mmultiple cases during diff later on
frames = cat(3, frames, frames(:,:,size(frames, 3)));
frames_reg = cat(3, frames_reg, frames_reg(:,:,size(frames_reg, 3)));

% perform blind deconvolution
for idx = 1:size(frames_reg, 3)
    idx
    img = frames_reg(:,:,idx);
    [J,PSF] = deconvblind(img, ones(size(img)));
    frames_deconvblind(:,:,idx) = J;
end

% write results
for idx = 1:size(frames_reg, 3) - 1
    
    % for comparison: write area of interest from original
    mkdir('frames_orig');
    filename = sprintf('frames_orig/frame_%03d.png', idx);
    data0raw = frames(:,:,idx);
    data0 = uint8(repmat(data0raw,[1 1 3]));
    imwrite(data0,filename,'png');
    
    % for comparison: write result from turbulence removal
    mkdir('frames_reg');
    filename = sprintf('frames_reg/frame_%03d.png', idx);
    data1raw = frames_reg(E+1:end-E,E+1:end-E,idx);
    data1 = uint8(repmat(data1raw,[1 1 3]));
    imwrite(data1,filename,'png');
    
    mkdir('frames_deconvblind');
    filename = sprintf('frames_deconvblind/frame_%03d.png', idx);
    data2raw = frames_deconvblind(E+1:end-E,E+1:end-E,idx);
    data2 = uint8(repmat(data2raw,[1 1 3]));
    imwrite(data2,filename,'png');
    
    mkdir('frames_diffs1');
    filename = sprintf('frames_diffs1/frame_%03d.png', idx);
    data3raw = frames(:,:,idx);
    data3raw1 = frames(:,:,idx+1);
    data3raw = abs(data3raw - data3raw1) * 50;
    data3 = uint8(repmat(data3raw,[1 1 3]));
    imwrite(data3,filename,'png');
    
    mkdir('frames_diffs2');
    filename = sprintf('frames_diffs2/frame_%03d.png', idx);
    data4raw = frames_reg(E+1:end-E,E+1:end-E,idx);
    data4raw1 = frames_reg(E+1:end-E,E+1:end-E,idx+1);
    data4raw = abs(data4raw - data4raw1) * 50;
    data4 = uint8(repmat(data4raw,[1 1 3]));
    imwrite(data4,filename,'png');
    
    mkdir('frames_diffs1');
    filename = sprintf('frames_diffs1/frame_%03d.png', idx);
    data5raw = frames_deconvblind(E+1:end-E,E+1:end-E,idx);
    data5raw1 = frames_deconvblind(E+1:end-E,E+1:end-E,idx+1);
    data5raw = abs(data5raw - data5raw1) * 50;
    data5 = uint8(repmat(data5raw,[1 1 3]));
    imwrite(data5,filename,'png');
    
    filenameSothh = sprintf('frames_Sothh/frame_%03d.png', idx);
    imgSothh = imread(filenameSothh);
    imgSothh = rgb2gray(imgSothh);
    filenameSothh1 = sprintf('frames_Sothh/frame_%03d.png', idx+1);
    imgSothh1 = imread(filenameSothh1);
    imgSothh1 = rgb2gray(imgSothh1);
    diffSothh = abs(imgSothh - imgSothh1) * 50;
    
    % final (for concat with ffmpeg)
    mkdir('frames_allthree');
    filename = sprintf('frames_allthree/frame_%03d.png', idx);
    data0raw = insertText(data0raw,[0 0],'raw video','BoxColor','white');
    data0raw = data0raw(:,:,1);
    data1raw = insertText(uint8(data1raw),[0 0],'w/o turbulences','BoxColor','white');
    data1raw = data1raw(:,:,1);
    data2raw = insertText(uint8(data2raw),[0 0],'deconvolved','BoxColor','white');
    data2raw = data2raw(:,:,1);
    data3raw = insertText(uint8(data3raw),[0 0],'diff w/ next frame','BoxColor','white');
    data3raw = data3raw(:,:,1);
    imgSothh = insertText(uint8(imgSothh),[0 0],'deconv by /u/Sothh','BoxColor','white');
    imgSothh = imgSothh(:,:,1);
    dataraw = [data0raw data1raw data2raw imgSothh; data3raw data4raw data5raw diffSothh];
    data = uint8(repmat(dataraw,[1 1 3]));
    imwrite(data,filename,'png');
end

% remove additional frame
frames = frames(:,:,1:size(frames, 3)-1);
frames_reg = frames_reg(:,:,1:size(frames_reg, 3)-1);