function [zh, zt] = SAIF(y,z,sigma,kernel,estimator)
% Spatially Adaptive Iterative Filtering %
% H. Talebi, X. Zhu and P. Milanfar, "How to SAIF-ly Boost Denoising Performance", 
% IEEE Transactions on Image Processing, vol 22, No. 4, pp. 1470-1485, April 2013.
%
% Input: y: noisy image
%        z: clean image
%        sigma: noise standard deviation
%        kernel: baseline kernel and also the pre-filter:
%                Non-loal Means: 'NLM'
%                Locally adaptive kernel regression: 'LRK'
%                Bilateral: 'BLT'
%        estimator: Plug-in: 'PLGN' 
%                   Stein's Unbiased Estimator: 'SURE'
    
[M,N] = size(z);

%%%%%%%%%%% parameters %%%%%%%%%%%%%%

R = 5; % Patch Radius
L = 2*R+1; % Patch Size = L x L
ksize = L;
skip=3; % Patch Skip Size
maxk=100; % Maximum Iteration Number
stepk=0.05; % Iteration Step Size


%%%%%%%%%%%%% padding %%%%%%%%%%%%%%

zz = padarray(z,[R R],'symmetric');
yy = padarray(y,[R R],'symmetric');

%%%%%%%%% Estimator Selection %%%%%%%%%%%%%%%

if (exist('estimator') ~= 1),
    kernel               = 'PLGN'; % default kernel is NLM
end


%%%%%%%%% Kernel Selection %%%%%%%%%%%%%%%
tic
if (exist('kernel') ~= 1),
    kernel               = 'NLM'; % default kernel is NLM
end

if kernel=='NLM'
disp(sprintf('Prefiltering by Non-Local Means...'))
    
%NLM
h0 = 0.32*sigma;
h1 =0.43*sigma;
ksz = 7;
wrad = 10; % patch radius
[K_rec,zt] = NLM(yy,wrad,h0,ksz);
clear K_rec;
[K_rec] = NLM(zt,R,h1,ksz);

elseif kernel=='LRK'
disp(sprintf('Prefiltering by LARK...'))
   
%LARK
wsize = 7;
h0 = 0.18*sigma;
h1 = 0.25*sigma;
ksize0 = 21;
[K_rec,zt] = LARK(yy,ksize0,wsize,h0);
clear K_rec;
[K_rec] = LARK(zt,ksize,wsize,h1);


elseif kernel=='BLT'
disp(sprintf('Prefiltering by Bilateral...'))
 
%BL
hs = 2;
hr = 2*sigma;
[K_rec,zt] = BLT(yy,R,hs,hr);
clear K_rec;
[K_rec] = BLT(zt,R,hs,hr);

else 
    
        error('SAIF accepts LRK, NLM or BLT as the input Kernel');

end
time = toc;
disp(sprintf('Elapsed Time = %.2f sec', time))

zt(zt>255) = 255; zt(zt<0) = 0;
ztt = zt(R+1:end-R,R+1:end-R);
et = ztt - z;
MSE_PreFilter = mean(et(:).^2);
PSNR_PreFilter = 10*log10((255^2)/MSE_PreFilter);
disp(sprintf('PreFilter PSNR = %.2f dB', PSNR_PreFilter))

rr = [1:skip:M];
cc = [1:skip:N];
  
%  matlabpool(2)
% matlabpool close

disp(sprintf('Filtering by SAIF...'))

    tic
for m=1:length(rr)   
%      tic   
    for n=1:length(cc)        
         
   loc_p = [rr(m)+R , cc(n)+R]; % patch location  
   [W] = WM(loc_p,ksize,K_rec,R,kernel); % patch-based filter   

   [U,s,V]=svd(W); % svd decomposition 
   lambda = diag(s); % 

   yp = yy(loc_p(1)-R :loc_p(1)+R , loc_p(2)-R :loc_p(2)+R);
   yv = yp(:);

   ztp = zt(loc_p(1)-R :loc_p(1)+R , loc_p(2)-R :loc_p(2)+R);
   ztv = ztp(:);

% pre-filtered projected patches
b0e = V'*ztv;

% noisy pojected patches
b0h = V'*yv;


if estimator == 'PLGN'
    
% Plug-in estimator %
[k_plug_diff(m,n),minMSE_diff] = plug_MSE_diff(lambda,b0e,sigma,maxk,stepk); % diffusion
[k_plug_twic(m,n),minMSE_twic] = plug_MSE_twic(lambda,b0e,sigma,maxk,stepk); % twicing

if minMSE_diff < minMSE_twic  
%diffusion    
s_plug = diag(lambda.^k_plug_diff(m,n)); 
W_plug = real(V*s_plug*V');
else
%twicing
s_plug = diag(1-(1-lambda).^(k_plug_twic(m,n)+1));
W_plug = real(V*s_plug*V');
end

% weight aggregation (plug-in) %
T=0.1*sigma^2;
zh_plug = W_plug*yv;
zhp_plug = reshape(zh_plug, L, L);

cov_zh_plug = (sigma^2)*(W_plug^2);
var_plug = diag(cov_zh_plug); % variance of estimation for each pixel
bias2_plug = ((W_plug-eye(L^2))*ztv).^2; % bias2 of estimation for each pixel
risk = real(bias2_plug + var_plug);
weight = exp(-risk/T);

patch_weight(:,:,m,n) = reshape(weight, L, L);
patch(:,:,m,n) = patch_weight(:,:,m,n).*zhp_plug;

else
% SURE estimator%
[k_sure_diff(m,n),minMSE_diff] = SURE_MSE_diff(lambda,b0h,sigma,maxk,stepk);
[k_sure_twic(m,n),minMSE_twic] = SURE_MSE_twic(lambda,b0h,sigma,maxk,stepk);

if minMSE_diff < minMSE_twic  
%diffusion    
s_sure = diag(lambda.^k_sure_diff(m,n)); 
W_sure = real(V*s_sure*V');
else
%twicing
s_sure = diag(1-(1-lambda).^(k_sure_twic(m,n)+1));
W_sure = real(V*s_sure*V');
end

% weight aggregation (sure) %
P = 0.05*sigma^2;
cov_zh_sure = (sigma^2)*(W_sure^2);
var_sure = diag(cov_zh_sure);
% weight_sure = 1./var_sure;
weight = exp(-var_sure/P);
zh_sure = W_sure*yv;
zhp_sure = reshape(zh_sure, L, L);
patch_weight(:,:,m,n) = reshape(weight, L, L);
patch(:,:,m,n) = patch_weight(:,:,m,n).*zhp_sure;

end

    end
%    m
%    time=toc
      
end
time = toc;
disp(sprintf('Elapsed Time = %.2f sec', time))
out = aggreg(patch,patch_weight,M,N,R,skip);
clear patch;
clear patch_weight;
clear K_rec;

zt = zt(R+1:end-R,R+1:end-R);

out(out>255) = 255; out(out<0) = 0;
zh = real(out);

end
 