Code covered by the BSD License

# GUI for denoising video signals with Kalman filter

### Youssef Khmou (view profile)

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] .
%
% 2011/2012
%
% Example:
%        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(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

```