Code covered by the BSD License  

Highlights from
GUI for denoising video signals with Kalman filter

image thumbnail

GUI for denoising video signals with Kalman filter

by

 

filter grayscale video signals(avi || Webcam) with 1st Order Estimation & Kalman filter .

[X,P,K]=KALMANFILTER_gui(video,option)
function [X,P,K]=KALMANFILTER_gui(video,option)
% function [X,P,K] = KALMANFILTER_gui(video,option)
% Usage :
% X : the estimated state
% P : the estimated error covariance
% K : Kalman gain  0<= K <= 1.
%
% The system is defined as follows :
%
% X(n+1) =  A(n) .* X(n) + W(n) ,  with W(n) ~ N(0, Q) .
% Y(n+1) =  H(n) .* X(n) + V(n) ,  with V(n) ~ N(0, R) .
%
% The dynamic matrix A is specific for a every model  and relates the two 
% states X(n+1) and X(n). 
% The algorithm :
%
% K(n+1)     = P_est(n+1) ./ (P_est(n+1)+ Q) .
% P_est(n+1) = (I - K).* P_est(n) .
% X_est(n+1) = X_est(n) + K .*( Z(n) - H. *X_est(n)) .  
% option specifies how to estimate the mean state of pixels 
%
% option =1 means that the mean X is calculed on cube [3*3*3]
% option =2 means that the mean X is calculated on cube [3*3*2] 
% Generally the second option preserve edges better
% Class of the input Argument "video" : [double] .
%
% copyright : KHMOU Youssef
% 2011/2012
%
% Example:
%        load mri
%        D=squeeze(D);
%        D=im2double(D);
%        D=imnoise(D,'Gaussian',0,0.06);
%        [X,P,K]= KALMANFILTER(D,1);
%
% In case of Stack of size >= 640*480*20 , use the following functions
% to decompose the SIGNAL into 4 parts :  Decomposition & fusion
%
% [Bloc1, Bloc2 Bloc3, Bloc4] =  Decomposition( video);
 
%Bloc1=Estimation(Bloc1);
%Bloc2=Estimation(Bloc2);
%Bloc3=Estimation(Bloc3);
%Bloc4=Estimation(Bloc4);
 
% Bloc1=Perform_Filter(Bloc1);
%Bloc2=Perform_Filter(Bloc2);
%Bloc3=Perform_Filter(Bloc3);
%Bloc4=Perform_Filter(Bloc4);
 
%Y= fusion(Bloc1, Bloc2, Bloc3, Bloc4 );
 tic
video(:,:,2:end)    = Estimation(video(:,:,2:end),option);
[X,P,K]             = Perform_Filter(video);
 toc
%figure,
%plot(P,'LineWidth',2)
%xlabel('Frames');
%ylabel('Mean Matrix Covariance ');
%title('Mean Matrix Covariance P over filtered sequence ');
%figure,
%plot(K,'LineWidth',2)
%xlabel('Frames');
%ylabel('Mean of Gain');
%title(' Mean of Kalman Gain');
%================================Perform_Filter============================
function [X_est,PP,KK] = Perform_Filter(video)
         
% Variables
line             = size(video,1);
column           = size(video,2);
time             = size(video,3);
PP               = zeros(time,1);
KK               = zeros(time,1);
X_est            = double(zeros(size(video)));
I                = ones(line,column);

% Initial conditions .
%X1              = wiener2(video(:,:,1),[3 3]);
X1               = Initial_Estimation(video(:,:,1));
X_est(:,:,1)     = X1;
E                = X1-video(:,:,1);
Q                = rand(1)/10*ones(line,column);
R                = rand(1)/10*ones(line,column);
P_est            = cov(E);
if (line>column)
    delta        = line-column;
    portion      = P_est(1:delta,:);
    P_est        = [P_est;portion];
end
if(line<column)
    P_est        = P_est(1:line,:);
end
K                = P_est./(P_est+R);
PP(1)            = mean(P_est(:));
KK(1)            = mean(K(:));


% The  ALGORITHM
for i= 2 : time
    
    %A                          =  video(:,:,i)./video(:,:,i-1);
    %A(isnan(A))                =  0.5;
    %A(isinf(A))                =  0.5;
    X_pred                     =  X_est(:,:,i-1);
    X_pred(isnan(X_pred))      =  0.5;
    X_pred(isinf(X_pred))      =  0.5;
    %P_pred                     =  (A.*P_est.*A) + Q ; 
    P_pred                    =  P_est + Q ; 
    K                          =  P_pred./(P_pred + R);
    Correction                 =  K .*(video(:,:,i)-X_pred);
    X_est(:,:,i)               =  X_pred+Correction;
    % X_est(:,:,i)             =  imadjust(X_est(:,:,i));
    X_est(isnan(X_est(:,:,i))) =  0.5;
    X_est(isinf(X_est(:,:,i))) =  0.5;   
    P_est                      =  (I-K).*P_pred;
    PP(i)                      =  mean(P_pred(:));
    KK(i)                      =  mean(K(:));
  
end
% Elimnating the Nan and Inf values from the estimated State
%for i=2:line-1
 %   for j=2:column-1
 %       for k=1:time
 %            if(isnan(X_est(i,j,k)) || isinf(X_est(i,j,k)))

 %               X_est(i,j,k)=( X_est(i-1,j,k)+X_est(i+1,j,k)+...
 %                   X_est(i,j-1,k)+X_est(i,j+1,k) ) ./4;
 %           end
 %       end
 %   end
%end
% eliminating the INF AND NAN in the edges
X_est(1,:,:)     = X_est(2,:,:);
X_est(:,1,:)     = X_est(:,2,:);
X_est(line,:,:)  = X_est(line-1,:,:);
X_est(:,column,:)= X_est(:,column-1,:);
return
%====================================Estimation============================
function Y = Estimation(X,option)

% Normalize the data : covert data to Double
if ~isa(X,'double') 
    X=double(X)./255;
end
if max(X(:)) > 1.00
    X=X./255.00;
end

n = size(X,1);
p = size(X,2);
m = size(X,3);
Y = double(zeros(size(X)));
% Estimation
% Intial anf Final Estimations : X(t_0) and X(t_final)
for i=2:n-1
    for j=2:p-1
      %  Y(i,j,1)   = (  X(i-1,j,1)+X(i+1,j,1)+X(i,j-1,1)+X(i,j+1,1) ) ./ 4;
        Y(i,j,end) = (  X(i-1,j,m)+X(i+1,j,m)+X(i,j-1,m)+X(i,j+1,m) ) ./ 4;
    end
end

if option==1
        
for i=2:n-1
    for j=2:p-1
        for k=2:m-1
            Y(i,j,k)= (( X(i+1,j-1,k) + X(i+1,j,k)+ X(i+1,j+1,k)+X(i,j-1,k)+...
             X(i,j+1,k)+X(i-1,j-1,k)+X(i-1,j,k)+...
                X(i-1,j+1,k) ) + (X(i+1,j-1,k-1) + X(i+1,j,k-1)+ X(i+1,j+1,k-1)+X(i,j-1,k-1)+...
             X(i,j+1,k-1)+X(i-1,j-1,k-1)+X(i-1,j,k-1)+...
                X(i-1,j+1,k-1) ) + (X(i+1,j-1,k+1) + X(i+1,j,k+1)+ X(i+1,j+1,k+1)+X(i,j-1,k+1)+...
             X(i,j+1,k+1)+X(i-1,j-1,k+1)+X(i-1,j,k+1)+...
                X(i-1,j+1,k+1) )+X(i,j,k) )./27;
        end
    end
end


elseif option==2

    for i=2:n-1
        for j=2:p-1
            for k=2:m
                Y(i,j,k)=(( X(i+1,j-1,k) + X(i+1,j,k)+ X(i+1,j+1,k)+X(i,j-1,k)+...
             X(i,j+1,k)+X(i-1,j-1,k)+X(i-1,j,k)+...
                X(i-1,j+1,k) ) + (X(i+1,j-1,k-1) + X(i+1,j,k-1)+ X(i+1,j+1,k-1)+X(i,j-1,k-1)+...
             X(i,j+1,k-1)+X(i-1,j-1,k-1)+X(i-1,j,k-1)+...
                X(i-1,j+1,k-1) +X(i,j,k)) )./18;
            end
        end
    end
end

% Now all pixels of Y are estimated by 3*3*2 pixels of noisy X
% Temprary Stop
 
if (option~=1 &&option~=2)
    error(' Invalid Option for type of  estimation')
end
return
%=========================================================================
function Z = Initial_Estimation(X)

[n p]=size(X);
Z=double(zeros(size(X)));
for i=2:n-1
    for j=2:p-1
        Z(i,j)=X(i-1,j)+X(i+1,j)+X(i,j-1)+X(i,j+1);
        Z(i,j)=Z(i,j)./4;
    end
end
return

%function [Bloc1, Bloc2 Bloc3, Bloc4] = Decomposition( Matrix)

%[n p k]= size(Matrix);

%t1=n/2;
%t2=p/2;

%Bloc1=squeeze(Matrix(1:t1,1:t2,:));
%Bloc2=squeeze(Matrix(1:t1,t2-3:p,:));
%Bloc3=squeeze(Matrix(t1-3:n,1:t2 ,:));
%Bloc4=squeeze(Matrix(t1-3:n, t2-3:p,:));

%return

%==========================================================================
%function Y= fusion(Bloc1, Bloc2, Bloc3, Bloc4 )

%[n1 p1 k]=size(Bloc1);
%[n2 p2 k]=size(Bloc2);
%[n3 p3 k]=size(Bloc3);
%[n4 p4 k]=size(Bloc4);

%Y=double(zeros(2*n1, 2*p1, k));
%Y(1 : n1, 1 : p1, :)             = Bloc1;
%Y(1 : n1, p1-3 : 2*p1, :)       = Bloc2;
%Y(n1-3 : 2*n1, 1 : p1, :)        = Bloc3;
%Y(n1-3 : 2*n1, p1-3 : 2*p1,:) =Bloc4;
%return






    
    
    

Contact us