%**************************************************************************
% główna funkcja do testowania metod nieliniowej aproksymacji w wielu
% bazach
%**************************************************************************

clear all; close all;

measurements = [ 1000 5000:5000:65000];
GraphError =   [0.01   0.02  0.05  0.1   0.2   0.5   1    3   5    10  20  50  100];

X = double(imread('lena.tif'));

%% fantom 1
% fantom_rozdzielczosc;
% X1 = double(xL);
% a = zeros(112,400);
% X = [X1; a]; 
% b = zeros(512,112);
% X = [X b];
% X = double(X);
% % X = X(1:512,1:512);
% % imshow(X,[]);

%% fantom 2
% fantom_real_small_jasnosc;
%% obraz źródłowy
imshow(X,[]);
x = X(:);
n = size(X,1);
N = n*n;
%imshow(X,[]); title('source image');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%liczba uwzględnionych współczynników
numa = 6500; hh=0;%1300; %numa = int(0.1*N);
%hh = 100; numa=6500+hh;
disp('LiczbaWSP='); disp(numa);
%energy = 1000;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Metody
FOURIER = 1; DCT = 1; BDCT = 1; WAV = 1; CONT = 1;
UWAV = 1; CURV = 1; CWT = 1; SHEAR = 1; SURF = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

if FOURIER
%aproksymacja DCT
    A = fft2(X);

    a = A(:); 
    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(abs(a),'descend');
  %  for i=1
  %   if sum(absa(1:i)) > energy
    
    OMSORT = ind(1+hh:numa);
    c(OMSORT) = a(OMSORT);
    C = reshape(c,n,n);

    Y = real(ifft2(C));
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('Fourier approximation');
     disp('FourierAPP='); disp(psnr(X,YY)); disp(rozm/N);
     imwrite(uint8(YY),sprintf('atom_fur.pgm'),'pgm');
end

if DCT
%aproksymacja DCT
    A = dct2(X);
    a = A(:); 
    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(abs(a),'descend'); %nonlinear
    OMSORT = ind(1+hh:numa); %nonlinear
    c(OMSORT) = a(OMSORT); %nonlinear
    C = reshape(c,n,n); %nonlinear
%     lporder = bdct_linapprox_ordering(n, n); %linear
%     C(lporder(1+hh:numa)) = A(lporder(1+hh:numa)); %linear
%     c(1+hh:numa) = a(1+hh:numa);%linear

    Y = idct2(C);
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('DCT approximation');
     disp('DctAPP='); disp(psnr(X,YY)); disp(rozm/N);
     imwrite(uint8(YY),sprintf('atom_dct.pgm'),'pgm');
end

if BDCT
%aproksymacja DCT w blokach
    w = 32; %wielkość bloku
    bsiz = uint8(n/w);
    fL = zeros(n,n);

    for i=1:bsiz
        for j=1:bsiz
            seli = (i-1)*w+1:i*w;
            selj = (j-1)*w+1:j*w;
            fL(seli,selj) = dct2( X(seli,selj) );
        end
    end
    %A = dct2(X);

    a = fL(:); 
    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(abs(a),'descend');
    OMSORT = ind(1+hh:numa);
    c(OMSORT) = a(OMSORT); %nonlinear
    %c(1+hh:numa) = a(1+hh:numa);%linear
    C = reshape(c,n,n);
    Y = zeros(n,n);
    for i=1:bsiz
        for j=1:bsiz
            seli = (i-1)*w+1:i*w;
            selj = (j-1)*w+1:j*w;
            Y(seli,selj) = idct2( C(seli,selj) );
        end
    end
    
    YY = clamp(Y,0,255);
    
    figure, imshow(YY,[]); title('BDCT approximation');
     disp('BDctAPP='); disp(psnr(X,YY)); disp(rozm/N);
     imwrite(uint8(YY),sprintf('atom_bdct.pgm'),'pgm');
end

if WAV
    %aproksymacja falkowa
    wkind = 'ort';
    
    if strcmp(wkind,'ort')
        qmf = MakeONFilter('Coiflet',5);%'Battle',5);%'Vaidyanathan');%'Symmlet',10);%'Beylkin');%'Coiflet',5); %ortogonaly
        A = FWT2_PO(X,3,qmf);
    else
        [qmf,dqmf] = MakeBSFilter('CDF',[4 4]);%'Villasenor',5); %biortogonaly
        A = FWT2_SBS(X,3,qmf,dqmf);
    end

    a = A(:); 
    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(abs(a),'descend');
    OMSORT = ind(1+hh:numa);
    c(OMSORT) = a(OMSORT);
    C = reshape(c,n,n);

    if strcmp(wkind,'ort')
        Y = IWT2_PO(C,3,qmf);
    else
        Y = IWT2_SBS(C,3,qmf,dqmf);
    end
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('Wavelet approximation');
    disp('WaveletAPP='); disp(psnr(X,YY)); disp(rozm/N);
    disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(X(70:136,167:250)))- var(var(X(70:136,17:100))))/(var(var(X(70:136,167:250)))+ var(var(X(70:136,17:100)))),abs(var(var(YY(70:136,17:100)))- var(var(YY(70:136,167:250))))/(var(var(YY(70:136,167:250)))+ var(var(YY(70:136,17:100))))));
    imwrite(uint8(YY),sprintf('atom_fal.pgm'),'pgm');
end

if UWAV,
    %aproksymacja falkowa undecimated
    wkind = 'sym2';%'sym1';%'coif2';%'bior4.4';
    levels = 4;
    A = swt2(X,levels,wkind);
    
    a = A(:); 
    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(abs(a),'descend');
    OMSORT = ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
    c(OMSORT) = a(OMSORT);
    C = reshape(c,n,n,1+3*levels);

    Y = iswt2(C,wkind);
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('UnWavelet approximation');
    disp('UnWaveletAPP='); disp(psnr(X,YY)); disp(rozm/N);

    %mapa aproksymacji
% % %     A = swt2(X,levels,wkind);
% % %     
% % %     a = A(:); 
% % %     rozm = size(a,1); c = zeros(rozm,1);
% % %     [absa,ind] = sort(abs(a),'descend');
% % %     
% % %     R = zeros(n); M = zeros(n); S = X; P = zeros(n); K = zeros(n);  
% % %     pr1 = zeros(n); pr2 = zeros(n);
% % %     krok = round(rozm/50);
% % %     for numa = 100000:100000:rozm,
% % %         OMSORT = ind(1:numa); %OMSORT = ind(1:numa);
% % %         c(OMSORT) = a(OMSORT);
% % %         C = reshape(c,n,n,1+3*levels);
% % %         Y = iswt2(C,wkind);
% % %         YY = clamp(Y,0,255);
% % %         R = abs(S-YY); %T = abs(R-S);
% % %         M = max(M,R); S = YY; %macierz maksymalnych gradientów
% % %         pom = abs(X-YY) < 0.9*X; pom = pom | pr1;
% % %         P(xor(pom,pr1)) = numa; pr1 = pom;
% % %         
% % %         pom = abs(X-YY) < 0.1*X; pom = pom | pr2;
% % %         K(xor(pom,pr2)) = numa; pr2 = pom;       
% % % %         if R < 0.9*X && P == 0,
% % % %             P = numa;
% % % %         elseif R < 0.1*X && K == 0,
% % % %             K = numa;
% % % %         end
% % %     end
% % %     M = M./(abs(K-P)+1);
% % %     MN = M./max(max(M));
% % %     figure, imshow(M,[]); title('UnWavelet approximation map1');
% % %     figure, imshow(abs(K-P),[]); title('UnWavelet approximation map2');
% % %     figure, imshow(MN,[]); title('UnWavelet approximation map3');
%*****************
end
if CONT
    %aproksymacja conturletowa
    % Parameters
    pfilt = '9-7';
    dfilt = 'pkva';
    nlevs = [0, 0, 0, 3, 4]; %[0, 0, 4, 4, 5]

    A = pdfbdec(X, pfilt, dfilt, nlevs);
    [a, s] = pdfb2vec(A);

    rozm = size(a,2); c = zeros(1,rozm);
    [absa,ind] = sort(abs(a),'descend');
    OMSORT = ind(1+hh:numa);%ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
    c(OMSORT) = a(OMSORT);

    C = vec2pdfb(c, s);
    Y = pdfbrec(C, pfilt, dfilt);
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('Contourlet approximation');
    disp('ContourletAPP='); disp(psnr(X,YY)); disp(rozm/N);
    disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(X(70:136,167:250)))- var(var(X(70:136,17:100))))/(var(var(X(70:136,167:250)))+ var(var(X(70:136,17:100)))),abs(var(var(YY(70:136,17:100)))- var(var(YY(70:136,167:250))))/(var(var(YY(70:136,167:250)))+ var(var(YY(70:136,17:100))))));
    imwrite(uint8(YY),sprintf('atom_cont.pgm'),'pgm');
end

if CURV
    %aproksymacja curveletowa
    % Parameters

    C = fdct_wrapping(X,0,2,3,16); 

    nbscales = length(C);
    nbangles_coarse = length(C{2});
    nbangles = [1, nbangles_coarse .* 2.^(ceil((nbscales-(nbscales:-1:2))/2))];

    a = C{1}{1}(:);
    rozma = size(a,1); na = size(C{1}{1},1);
    for s = 2:length(C)
      for w = 1:length(C{s})
        b = C{s}{w}(:);
        a = [a;b(:)];
      end
    end

    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(abs(a),'descend');
    OMSORT = ind(1:numa);%ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
    c(OMSORT) = a(OMSORT);

    C{1}{1} = reshape(c(1:rozma),na,na); 
    rozmx = rozma;
    for s = 2:length(C)
      for w = 1:length(C{s})
        blokx = size(C{s}{w},1); bloky = size(C{s}{w},2);
        pocz = rozmx + 1;  kon = rozmx + blokx*bloky; 
        C{s}{w} = reshape(c(pocz:kon),blokx,bloky);
        rozmx = rozmx + blokx*bloky;
      end
    end

    nbscales = length(C);
    nbangles_coarse = length(C{2});
    nbangles = [1, nbangles_coarse .* 2.^(ceil((nbscales-(nbscales:-1:2))/2))];

    Y = real(ifdct_wrapping(C,0)); 
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('Curvelet approximation');
    disp('CurveletAPP='); disp(psnr(X,YY)); disp(rozm/N);
    disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(X(70:136,167:250)))- var(var(X(70:136,17:100))))/(var(var(X(70:136,167:250)))+ var(var(X(70:136,17:100)))),abs(var(var(YY(70:136,17:100)))- var(var(YY(70:136,167:250))))/(var(var(YY(70:136,167:250)))+ var(var(YY(70:136,17:100))))));
end

if CWT
    %aproksymacja falek zespolonych
    % Parameters
    biort = 'near_sym_b';%'legall';%'antonini'; %'near_sym_b'; % biort = 'antonini';
    qshift = 'qshift_06';%'qshift_d'; 
    gm = ones(6,4); gm(:,1)= 1;gm(:,2)= 1;gm(:,3)=1;gm(:,4)=1;
    wzmoc = 1.8;%1.13;
    levels = 8;%9;

    [Yl,Yh] = dtwavexfm2(X,levels,biort,qshift);
    %[Yl,Yh] = dtwavexfm2ap(X,4,biort,qshift); %znormalizowane filtry qshift

% % %     ar = Ylr(:); %aa =a;
% % %     %rozma = size(a,1); na = size(Yl,1);
% % %     for j=1:4
% % %         %wzm = wzmoc^j;
% % %         for k=1:6
% % %             b = Yhr{j,1}(:,:,k);
% % %             ar = [ar;b(:)];
% % %             %aa = [aa;b(:) * wzm];
% % %         end
% % %     end
    a = Yl(:); %aa =a;%
    rozma = size(a,1); na = size(Yl,1);
    for j=1:levels
        %wzm = wzmoc^j;
        for k=1:6
            b = Yh{j,1}(:,:,k);
            a = [a;b(:)];
            %aa = [aa;b(:) * wzm];
        end
    end

    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(a,'descend'); %a  sortowanie po zwiększonych współczynnikach
    OMSORT = ind(1+hh:numa);%ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
    c(OMSORT) = a(OMSORT);

    Yl = reshape(c(1:rozma),na,na);
    rozmx = rozma; blok = n;
    for j=1:levels
        blok = blok/2;
        for k=1:6
            pocz = rozmx + 1;  kon = rozmx + blok*blok; 
            b = reshape(c(pocz:kon),blok,blok);
            Yh{j,1}(:,:,k) = b;
            rozmx = rozmx + blok*blok;
        end
    end

    Y = dtwaveifm2(Yl,Yh,biort,qshift);%,gm);
    %YY = (Y-min(min(Y)))/(max(max(Y))-min(min(Y)))*255;
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('Complex Wavelet approximation');
    disp('DT-CWTAPP='); disp(psnr(X,YY)); disp(rozm/N*2); %2 bo zespolone
    disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(X(70:136,167:250)))- var(var(X(70:136,17:100))))/(var(var(X(70:136,167:250)))+ var(var(X(70:136,17:100)))),abs(var(var(YY(70:136,17:100)))- var(var(YY(70:136,167:250))))/(var(var(YY(70:136,167:250)))+ var(var(YY(70:136,17:100))))));
    imwrite(uint8(YY),sprintf('atom_cwt.pgm'),'pgm');
    %Z = histeq(uint8(Y));
    %figure, imshow(Z,[]); 


% % % %mapa aproksymacji
% % %     [Yl,Yh] = dtwavexfm2(X,levels,biort,qshift);
% % % 
% % %     a = Yl(:); 
% % %     rozma = size(a,1); na = size(Yl,1);
% % %     for j=1:levels
% % %         for k=1:6
% % %             b = Yh{j,1}(:,:,k);
% % %             a = [a;b(:)];
% % %         end
% % %     end
% % % 
% % %     rozm = size(a,1); c = zeros(rozm,1);
% % %     [absa,ind] = sort(a,'descend'); %a  sortowanie po zwiększonych współczynnikach
% % %     
% % %     R = zeros(n); M = zeros(n); S = X; P = zeros(n); K = zeros(n);  
% % %     pr1 = zeros(n); pr2 = zeros(n);
% % %     krok = round(rozm/1000); zasieg = 0.1*rozm;
% % %     for numa = krok:krok:zasieg,
% % %         OMSORT = ind(1:numa); %OMSORT = ind(1:numa);
% % %         c(OMSORT) = a(OMSORT);
% % %         Yl = reshape(c(1:rozma),na,na);
% % %         rozmx = rozma; blok = n;
% % %         for j=1:levels
% % %             blok = blok/2;
% % %             for k=1:6
% % %                 pocz = rozmx + 1;  kon = rozmx + blok*blok; 
% % %                 b = reshape(c(pocz:kon),blok,blok);
% % %                 Yh{j,1}(:,:,k) = b;
% % %                 rozmx = rozmx + blok*blok;
% % %             end
% % %         end
% % %         Y = dtwaveifm2(Yl,Yh,biort,qshift);%,gm);
% % %         YY = clamp(Y,0,255);
% % %         R = abs(S-YY); %T = abs(R-S);
% % %         M = max(M,R);% (M+R)/2;
% % %         S = YY; %macierz maksymalnych gradientów
% % %         pom = abs(X-YY)<0.95*X; pom = pom | pr1;
% % %         P(xor(pom,pr1)) = numa; pr1 = pom;
% % %         
% % %         pom = abs(X-YY) < 0.5*X; pom = pom | pr2;
% % %         K(xor(pom,pr2)) = numa; pr2 = pom;       
% % %     end
% % %     M = M./(abs(K-P)+1);%próba z mnożeniem
% % %     MN = M./max(max(M))*255;
% % %     figure, imshow(adapthisteq(uint8(M)),[]); title('UnWavelet approximation map1');
% % %     figure, imshow(adapthisteq(uint8(abs(K-P))),[]); title('UnWavelet approximation map2');
% % %     figure, imshow(adapthisteq(uint8(MN)),[]); title('UnWavelet approximation map3');
% % % %      Z = adapthisteq(uint8(M));
% % % %     figure, imshow(Z,[]); 
% % % %***********koniec mapy
end

if SHEAR
    %aproksymacja shearlets
    % Parameters

    [ST, Psi] = shearletTransformSpect(X);

    a = ST(:); 
    sz = size(ST);
    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(abs(a),'descend');
    OMSORT = ind(1:numa);%ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
    c(OMSORT) = a(OMSORT);
    C = reshape(c,sz(1),sz(2),sz(3));

    Y = inverseShearletTransformSpect(C,Psi);
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('Shearlet approximation');
    disp('ShearletAPP='); disp(psnr(X,YY)); disp(rozm/N);
    disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(X(70:136,167:250)))- var(var(X(70:136,17:100))))/(var(var(X(70:136,167:250)))+ var(var(X(70:136,17:100)))),abs(var(var(YY(70:136,17:100)))- var(var(YY(70:136,167:250))))/(var(var(YY(70:136,167:250)))+ var(var(YY(70:136,17:100))))));
end

if SURF
    %aproksymacja surfacelets
    % Parameters
    L = 2; % We will have 2^(L+1) different directional subbands at each scale
    % The order of the checkerboard filter banks. The higher bo is, the
    % sharper the frequency response is.
    bo = 12; 

%     % Different configurations for the directional decomposition
%     Level_64 = [-1 3 3; 3 -1 3; 3 3 -1]; % 3 * 64 directions
%     Level_16 = [-1 2 2; 2 -1 2; 2 2 -1]; % 3 * 64 directions
%     Level_4 =  [-1 1 1; 1 -1 1; 1 1 -1]; % 3 * 64 directions
%     Level_1 =  [-1 0 0; 0 -1 0; 0 0 -1]; % 3 directions (i.e. the hourglass decomposition)
% 
%     Pyr_mode = 1;%2;
%     %Lev_array = {Level_64, Level_16, 1, 1};
%     Lev_array = {Level_64, Level_64, Level_16, Level_4};
    
    for i = 1 : 2
        Lev_array{i} = [-1 0; 0 -1];
    end
        Lev_array{3} = [-1 L; L -1];

    Pyr_mode = 1;
    
    [C, Recinfo] = Surfdec(X, Pyr_mode, Lev_array, 'ritf', 'bo', bo);
    G = C;
    %licz = 0; %a=0;
    a = C{1}{1}{1}(:); r =  size( C{1}{1}{1},1); r = [r;size( C{1}{1}{1},2)];
    C{1}{1}{1}(:) = 0;
    for n = 1 : length(C) - 1
    for dd = 1 : length(C{n})
        for sd = 1 : length(C{n}{dd})
            if n ~= 1 || dd ~= 1 || sd ~= 1,
                b = C{n}{dd}{sd};  C{n}{dd}{sd}(:) = 0;
                r = [r;size(b,1)];r = [r;size(b,2)];
                a = [a; b(:)];
            end
        end
    end
    end

    rozm = size(a,1); c = zeros(rozm,1);
    [absa,ind] = sort(abs(a),'descend');
    OMSORT = ind(1:numa);%ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
    c(OMSORT) = a(OMSORT);
    %c = a; %dodatek
    licz = 1; pocz = 1; kon = 0;
    for n = 1 : length(C) - 1
    for dd = 1 : length(C{n})
        for sd = 1 : length(C{n}{dd})
            rozx = r(licz); rozy = r(licz+1); kon =  kon + rozx*rozy; 
            b = reshape(c(pocz:kon),rozx,rozy);
            C{n}{dd}{sd} = b;
            pocz =kon+1; licz = licz+2;
        end
    end
    end   
    
     
    Y = Surfrec(C, Recinfo);
    YY = clamp(Y,0,255);
    figure, imshow(YY,[]); title('Surface approximation');
    disp('SurfaceletAPP='); disp(psnr(X,YY)); disp(rozm/N);
    disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(X(70:136,167:250)))- var(var(X(70:136,17:100))))/(var(var(X(70:136,167:250)))+ var(var(X(70:136,17:100)))),abs(var(var(YY(70:136,17:100)))- var(var(YY(70:136,167:250))))/(var(var(YY(70:136,167:250)))+ var(var(YY(70:136,17:100))))));
    
%filtracja wyjścia 

% % % f4=[
% % % 0     0    -1    -2    -2    -2    -1     0     0
% % % 0    -2    -3    -5    -5    -5    -3    -2     0
% % % -1    -3    -5    -3     0    -3    -5    -3    -1
% % % -2    -5    -3    11    21    11    -3    -5    -2
% % % -2    -5     0    21    40    21     0    -5    -2
% % % -2    -5    -3    11    21    11    -3    -5    -2
% % % -1    -3    -5    -3     0    -3    -5    -3    -1
% % % 0    -2    -3    -5    -5    -5    -3    -2     0
% % % 0     0    -1    -2    -2    -2    -1     0     0
% % % ];%/170;
% % % 
% % % %Y4=uint8(conv2(X,f4,'valid'));
% % % %Y4=conv2(YY,f4,'valid');
% % % Y4 = ncanny(YY,'canny');
% % % figure, imshow(adapthisteq(uint8(Y4)),[]); title('Filtered Surface approximation');    
    
end
% % % %%%% ocena szybkości opadania
% % % 
% % % % linear
% % % v_lin = cumsum( abs(a(end:-1:1)).^2 );
% % % v_lin = v_lin(end:-1:1);
% % % % non-linear
% % % v = sort( abs(a) );
% % % v_nl = cumsum( v.^2 );
% % % v_nl = v_nl(end:-1:1);
% % % % log/log plot of the errors
% % % clf; loglog( 1:length(v_lin), v_lin, 1:length(v_lin), v_nl );
% % % legend('linear', 'non-linear');
% % % axis([1 length(v_lin) 1e-5 max(v_lin) ]);