% Anisotropic LPA-ICI Denoising of PSF estimate
%
% The files are prepared in Tampere University of Technology, Institute of
% Information Technology, 2006.
%
% Vladimir Katkovnik, Alessandro Foi, Dmitriy Paliy*
% *e-mail: dmitriy.paliy@tut.fi

%---------------------------------------------------------
% LPA WINDOWS PARAMETERS
%---------------------------------------------------------
ndirPSF=4; % number of directions
lenhPSF=length(h1_PSF); % number of scales in RI

%---------------------------------------------------------
% ICI threshold
%---------------------------------------------------------
% GammaParameterPSF = [1.4]; % this parameter always fixed and does not 
%                            % depend on the color channel

%---------------------------------------------------------
% Create LPA comvolution kernels
%---------------------------------------------------------
if sst-sst_filtering_v==1 & ccolor==1 & j==1,
    TYPE = 10;  window_type = 112; directional_resolution = ndirPSF; sig_winds=[ones(size(h1_PSF)); ones(size(h2_PSF))];    % Gaussian parameter

    [kernels0, kernels_higher_order0] = function_CreateLPAKernels([0 0],h1_PSF,h2_PSF,TYPE,window_type,directional_resolution,sig_winds,1);
    
    % save PSFkernels kernels0 kernels_higher_order0
    PSFKernels0 = kernels0;
    PSFKernelsHO0 = kernels_higher_order0;
else
    % load PSFkernels
    kernels0 = PSFKernels0;
    kernels_higher_order0 = PSFKernelsHO0;
end

lenh=lenhPSF;
directional_resolution=ndirPSF;

clear yh_RI stdh_RI var_opt_Q
[size_z_1,size_z_2] = size(vvv);
%---------------------------------------------------------
% LPA-ICI Denoising of PSF starts...
%---------------------------------------------------------
y_hat_RI=zeros(size_z_1,size_z_2);
var_inv=zeros(size_z_1,size_z_2);
%%%%%%%%%%%%%%%% DIRECTIONAL LPA %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for s1=1:ndirPSF  % cycle along the directions
    for s=1:lenhPSF,
        gh = kernels_higher_order0{s1,s,1}(:,:,1);   %gets single kernel from the cell array
        ghorigin(s1,s)=gh((end+1)/2,(end+1)/2);
        bound1=min([(find(sum(gh~=0,2)));abs(find(sum(gh~=0,2))-size(gh,1)-1)]); % removes unnecessary zeroes
        bound2=min([(find(sum(gh~=0,1))),abs(find(sum(gh~=0,1))-size(gh,2)-1)]); % removes unnecessary zeroes
        gh=gh(bound1:size(gh,1)-bound1+1,bound2:size(gh,2)-bound2+1);            % removes unnecessary zeroes

        %%%%%%%% LPA  %%%%%%%%%%
        % LPA Filtering (in spatial domain)        
        yh_RI(:,:,s)=conv2(vvv,gh,'same');
        % Standard deviation of LPA estimate
        stdh_RI(:,:,s)=repmat(sqrt(sum(sum(gh.^2)))*dev1,[size_z_1,size_z_2]);
    end %%%%%%%%%%%%%%% end for H  %%%%%%%%%%%%%%%%%%%%
    %%%%%% ICI %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [YICI_RIT,h_optRI,std_optRI1]=function_ICI(yh_RI,stdh_RI,GammaParameterPSF,2*(s1-1)*pi/directional_resolution);

    y_hat_Q_PSF(:,:,s1) = YICI_RIT; %% ADAPTIVE DIRECTIONAL ESTIMATES
    var_opt_Q(:,:,s1)   = (std_optRI1.^2+eps);  %%% VARIANCES OF THE ADAPTIVE DIRECTIONAL ESTIMATES
    h_opt_Q_PSF(:,:,s1) = h_optRI; %% STORES RESULTS OF DIRECTIONAL ADAPTIVE SCALES
    y_hat_RI            = y_hat_RI+y_hat_Q_PSF(:,:,s1)./var_opt_Q(:,:,s1);            %% FUSING WITH ADAPTIVE WEIGHTS %%%%%
    var_inv             = var_inv+1./var_opt_Q(:,:,s1);   %% SUM OF INVERSE VARIANCES DENOMINATOR FOR CONVEXIFICATION OF ADAPTIVE LINEAR COMBINATION
end   %%% END THETA LOOP

y_hat_PSF=y_hat_RI./var_inv; % final estimate
%%%%%%%%%%%%% END OF ANISOTROPIC LPA-ICI %%%%%%%%%%%%%%%%