% Anisotropic LPA-ICI Denoising of signal 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
%---------------------------------------------------------
ndirRI=8; % number of directions
lenhRI=length(h1_Im); % number of scales in RI

%---------------------------------------------------------
% ICI threshold
%---------------------------------------------------------
GammaParameterRI  = GammaSignal(ccolor);

%---------------------------------------------------------
% Create LPA comvolution kernels
%---------------------------------------------------------
if sst-sst_filtering_y==1 & ccolor==1,
    TYPE = 10;  window_type = 112; directional_resolution = ndirRI; sig_winds=[ones(size(h1_Im)); ones(size(h2_Im))];    % Gaussian parameter

    [kernels0, kernels_higher_order0] = function_CreateLPAKernels([0 0],h1_Im,h2_Im,TYPE,window_type,directional_resolution,sig_winds,1);

    % save SignalKernels kernels0 kernels_higher_order0
    SignalKernels0 = kernels0;
    SignalKernelsHO0 = kernels_higher_order0;
else
    % load SignalKernels
    kernels0 = SignalKernels0;
    kernels_higher_order0 = SignalKernelsHO0;
end

lenh=lenhRI;
directional_resolution=ndirRI;

clear yh_RI stdh_RI var_opt_Q

y_hat_RI = zeros(size_z_1,size_z_2);
var_inv  = zeros(size_z_1,size_z_2);

%---------------------------------------------------------
% LPA-ICI Signal Denoising starts...
%---------------------------------------------------------
for s1=1:ndirRI  % cycle along the directions
    for s=1:lenhRI,
        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(y_est-100000,gh,'same')+100000;   
        % Standard deviation of LPA estimate
        stdh_RI(:,:,s) = repmat(sqrt(sum(sum(gh.^2)))*dev2,[size_z_1,size_z_2]);
    end %%%%%%%%%%%%%%% end for H  %%%%%%%%%%%%%%%%%%%%

    %%%%%% ICI %%%%%%%%%%%%%%%%%%%%%%%%%%%%
    [YICI_RIT,h_optRI,std_optRI1]=function_ICI(yh_RI,stdh_RI,GammaParameterRI,2*(s1-1)*pi/directional_resolution);
    % YICI_RIT          = max(0,min(1,YICI_RIT)); %% impose [0,1] constraint
    y_hat_Q_y(:,:,s1) = YICI_RIT; %% ADAPTIVE DIRECTIONAL ESTIMATES
    var_opt_Q(:,:,s1) = (std_optRI1.^2+eps);  %%% VARIANCES OF THE ADAPTIVE DIRECTIONAL ESTIMATES
    h_opt_Q_y(:,:,s1) = h_optRI; %% STORES RESULTS OF DIRECTIONAL ADAPTIVE SCALES
    y_hat_RI          = y_hat_RI+y_hat_Q_y(:,:,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_est_y = y_hat_RI./var_inv; % final estimate
%%%%%%%%%%%%% END OF ANISOTROPIC LPA-ICI %%%%%%%%%%%