雪花飘飘

前言

去雪

本文主要是借用卡尔曼滤波器的特性来追踪雪花，从而达到去除静态背景图中雪花的目的。
关于卡尔曼滤波器的原理，参考这里。
这里直接上代码，考虑到矩阵操作的便利性，所以选择了matlab进行试验，工程性选择还是应该首选opencv和c++

matlab代码

% demo.m

filename = 'snow1.avi';
[F1,E1] = Kalman(filename);
N = size(F1,2);
WriterObj = VideoWriter([filename(1:end-4) '_kalman_removal.avi']);
open(WriterObj);

for i = 1:N-1
    // 展示结果
    subplot(1,2,1);imshow(F1{i});
    subplot(1,2,2);imshow(E1{i});
    pause(0.001); 
    // 原始图片和结果图片上下拼接   
    R = cat(2,F1{i},E1{i});
    // 回写文件
    writeVideo(WriterObj,R);
end
close(WriterObj);

卡尔曼滤波

% Kalman.m

function [Frames,Estimates] = Kalman(filename)

    % 读取视频文件对象
    ReaderObj = VideoReader(filename); 
    
    % 获取视频帧数
    N = get(ReaderObj, 'NumberOfFrames'); 
    % 尾帧裁剪
    N = N-10;
    
    % 视频数据，一共4个维度：宽度、高度、RGB、帧数
    VideoFrames = read(ReaderObj);
    
    % 获取视频帧数
    sz = size(VideoFrames(:,:,:,1));
    
    % 抽取亮度数据
    light = zeros(sz(1),sz(2),N);
    for i = 1:N
        % 转换颜色空间：RGB->HSI
        HSI = rgb2hsv(VideoFrames(:,:,:,i));
        % 这里只关心第三维的亮度，色调、饱和度不关心
        light(:,:,i) = HSI(:,:,3); 
    end

    % 过程噪声
    W = randn(sz(1),sz(2),N) * 0.2;
    % 过程噪声的标准差和协方差
    R = std(W,0,3);Rww = R.^2;
    
    % 测量噪声
    V = randn(sz(1),sz(2),N) * 5;
    % 测量噪声的标准差和协方差
    Q = std(V,0,3);Rvv = Q.^2;
    
    % 系统参数
    A = 1;
    % 测量系统参数
    H = 1;
    
    % X代表有过程噪声的亮度数据，Y代表原来的亮度数据
    X{1} = light(:,:,1) * 255;
    Y{1} = X{1};
    for k = 2:N
        X{k} = A*X{k-1} + W(:,:,k);
        Y{k} = light(:,:,k)*255;
    end

    % 协方差初始值
    prior = ones(sz(1),sz(2)) * 2;

    % 首帧数据直接取原数据
    s{1} = X{1};

    for t = 2:N
        % 计算t时刻的协方差，有系统噪声
        p = A.^2 * prior + Rww;
        
        % 计算卡尔曼增益
        Kg = H*p ./ ((H.^2)*p + Rvv);
        
        % 更新t时刻的结果
        s{t} = A*s{t-1} + Kg .* (Y{t} - A*H*s{t-1});
        
        % 更新协方差
        prior = p - H*Kg.*p;
    end
    
    % 拼接色调和饱和度，转换颜色空间，并输出
    for t = 1:N
        HSI = rgb2hsi(VideoFrames(:,:,:,i));
        s{t} = hsv2rgb(HSI(:,:,1),HSI(:,:,2),im2double(s{t}/255));
        Y{t} = hsv2rgb(HSI(:,:,1),HSI(:,:,2),im2double(Y{t}/255));
    end
    Frames = Y;
    Estimates = s;
end

效果图

总结

先写到这里，后续有时间的话，会增添opencv代码