function FilteredI = EIGAveF(ImageR,Delta,D)
%function FilteredI = EIGAveF(ImageR,Delta,D)

%Code for
%Latorre-Carmona, P., Miñana, JJ. & Morillas, S. 
%Colour image denoising by eigenvector analysis of neighbourhood colour samples. 
%SIViP 14, 483–490 (2020). https://doi.org/10.1007/s11760-019-01575-5
%ImageR is the noisy Image
%Delta fixes the filterin window of size (Delta+1)x(Delta+1), Delta=1,2,...

[a b c]= size(ImageR);
FilteredI = ImageR;
NDatosWindow = (2*Delta+1)^2;
PosCenter = ceil(NDatosWindow/2);
%loop through the whole image
for i=(Delta+1):1:(a-Delta)
    for j=(Delta+1):1:(b-Delta)
     
        %Noisy data
        WindowRR = ImageR(i-Delta:i+Delta,j-Delta:j+Delta,1);
        WindowRG = ImageR(i-Delta:i+Delta,j-Delta:j+Delta,2);
        WindowRB = ImageR(i-Delta:i+Delta,j-Delta:j+Delta,3);
        
        %build data matrix to perform later eigenvector analysis
        VectorRR = double(Matrix2Vector(WindowRR));
        VectorRG = double(Matrix2Vector(WindowRG));
        VectorRB = double(Matrix2Vector(WindowRB));
        %means
        MeanR = mean(VectorRR);
        MeanG = mean(VectorRG);
        MeanB = mean(VectorRB);
        %Matrix
        DatosR = double([VectorRR; VectorRG; VectorRB ]');
        
        %Color vector processed in this iteration
        ColorVectorIJ = [DatosR(PosCenter,1) DatosR(PosCenter,2) DatosR(PosCenter,3)];        
        
        %Center data around means
        DatosR = double([(VectorRR - MeanR); (VectorRG - MeanG); (VectorRB - MeanB) ]');
                
        %Move to EIG space (reference system)
        %EIG: eigenvectors and eigenvalues
        %Solves Eqs (1-3) in the paper
        [EIGVEC EIGVAL] = eig(cov(DatosR));
        
        %Rotate data to EIGVEC reference system
        %Eq (4) in the paper
        DatosNuevosR = double(DatosR)*EIGVEC;

        %Vector to filter in EIGVEC space
        VectorSVD = [ DatosNuevosR(PosCenter,1) DatosNuevosR(PosCenter,2) DatosNuevosR(PosCenter,3)];

        %Compute standard deviation of each variable in EIGVEC space 
        %In the paper these are sigma(U^{i,j,k}) 
        stdDatosNuevosR = std(DatosNuevosR);
        
        %Normalize standard deviations: Eq (7) in the paper.       
        sumstdDatosNuevosR = sum(stdDatosNuevosR);
        if ( sumstdDatosNuevosR == 0) %it is unlikely but just in case...
            stdNormDatosNuevosR = [1/3 1/3 1/3];
        else
            stdNormDatosNuevosR = stdDatosNuevosR ./ sumstdDatosNuevosR; %Eq (7)
        end
        %add tiny amount just in case one of them is zero
        stdNormDatosNuevosR = stdNormDatosNuevosR + 0.0000000001; 
        
        
        
        %Processing in EIGVEC space Eq. (6)
        %First compute weights in Eq (8)
        %Weights, Numerator and denominator for each channel
        NumD1 = 0;
        DenD1 = 0;
        NumD2 = 0;
        DenD2 = 0;
        NumD3 = 0;
        DenD3 = 0;
        WeightsD1 = [];
        WeightsD2 = [];
        WeightsD3 = [];
        %loop through the data
        for k=1:1:NDatosWindow
            Weight = 1/exp(abs(VectorSVD(1)-DatosNuevosR(k,1))/(D/stdNormDatosNuevosR(1)));
            WeightsD1 = [WeightsD1 Weight];
            NumD1 =  NumD1 + Weight.*DatosNuevosR(k,1);
            DenD1 = DenD1 + Weight; 

            Weight = 1/exp(abs(VectorSVD(2)-DatosNuevosR(k,2))/(D/stdNormDatosNuevosR(2)));
            WeightsD2 = [WeightsD2 Weight];
            NumD2 =  NumD2 + Weight.*DatosNuevosR(k,2);
            DenD2 = DenD2 + Weight; 

            Weight = 1/exp(abs(VectorSVD(3)-DatosNuevosR(k,3))/(D/stdNormDatosNuevosR(3)));
            WeightsD3 = [WeightsD3 Weight];
            NumD3 =  NumD3 + Weight.*DatosNuevosR(k,3);
            DenD3 = DenD3 + Weight; 
            
        end
        %Cocient
        VectorSVDFiltrado = [ NumD1/DenD1 NumD2/DenD2 NumD3/DenD3];
        
        %Back to RGB from EIGVEC space Eq. (5): inverse rotation plus add
        %means 
        ColorVectorIJFiltrado = VectorSVDFiltrado*inv(EIGVEC) + [MeanR MeanG MeanB];
        
        %Assign result to current filtered image position
        FilteredI(i,j,:) = ColorVectorIJFiltrado; 
        
%         %You can always see results in every pixel and pause to see how
%         it is working 
%         DatosR' 
%         DatosNuevosR'
%         [WeightsD1;WeightsD2;WeightsD3]
%         
%         [stdNormDatosNuevosR;VectorSVD; VectorSVDFiltrado; ColorVectorIJ; ColorVectorIJFiltrado]'
%         
%         %Show data
%         subplot(1,3,1),imshow(ImageR(i-Delta:i+Delta,j-Delta:j+Delta,:))
%         subplot(1,3,2),imshow(FilteredI(i-Delta:i+Delta,j-Delta:j+Delta,:))
%         subplot(1,3,3),plot3(DatosR(:,1),DatosR(:,2),DatosR(:,3),'b.'),axis([-50 50 -50 50 -50 50]), grid on, hold on
%         %center
%         subplot(1,3,3),plot3(DatosR(PosCenter,1),DatosR(PosCenter,2),DatosR(PosCenter,3),'ko'),axis([-50 50 -50 50 -50 50])     
%         %center filtered
%         subplot(1,3,3),plot3(ColorVectorIJFiltrado(1),ColorVectorIJFiltrado(2),ColorVectorIJFiltrado(3),'bo'),axis([-50 50 -50 50 -50 50])             
%         %noisy rotated
%         subplot(1,3,3),plot3(DatosNuevosR(:,1),DatosNuevosR(:,2),DatosNuevosR(:,3),'r.'),axis([-50 50 -50 50 -50 50])
%         %center
%         subplot(1,3,3),plot3(DatosNuevosR(PosCenter,1),DatosNuevosR(PosCenter,2),DatosNuevosR(PosCenter,3),'ko'),axis([-50 50 -50 50 -50 50])        
%         %center filtered
%         subplot(1,3,3),plot3(VectorSVDFiltrado(1),VectorSVDFiltrado(2),VectorSVDFiltrado(3),'ro'),axis([-50 50 -50 50 -50 50])             
%         hold off
% 
%         pause
    end
end