function [appwalki] = image_approx(kolorimage, klasa, metoda,factx)
%kolorimage - obraz wejściowy; klasa - grupa metod; metoda - konkretna; metoda,factx - współczynnik aproksymacji(opcjonalny)
kolorimageo = kolorimage;
rescaling = 0;
if klasa == 1 || klasa == 5
    calkowita = log2(size(kolorimage,1));
    if floor(calkowita) ~= calkowita %** dopasowanie rozmiaru obrazów do potęgi 2
        if metoda == 4 || metoda == 5 || metoda == 6 || metoda == 9 || metoda == 10 || klasa == 5
            kolorimageo = kolorimage;
            sizev = size(kolorimage,1); sizeh = size(kolorimage,2);
            if sizev > 2048 || sizeh > 2048 
                a = padarray(kolorimage,[(4096-sizev)/2 (4096-sizeh)/2],0,'both'); kolorimage = a;
                % a = kolorimage((size(kolorimage,1)-1024)/2:(size(kolorimage,1)-1024)/2+1023,:); kolorimage = a;
            elseif sizev > 1024 || sizeh > 1024 
                a = padarray(kolorimage,[(2048-sizev)/2 (2048-sizeh)/2],0,'both'); kolorimage = a;
                % a = kolorimage((size(kolorimage,1)-1024)/2:(size(kolorimage,1)-1024)/2+1023,:); kolorimage = a;
            elseif sizev > 512 || sizeh > 512
                a = padarray(kolorimage,[(1024-sizev)/2 (1024-sizeh)/2],0,'both'); kolorimage = a;
                %a = kolorimage((size(kolorimage,1)-512)/2:(size(kolorimage,1)-512)/2+511,:); kolorimage = a;
            elseif sizev > 256 || sizeh > 256
                a = padarray(kolorimage,[(512-sizev)/2 (512-sizeh)/2],0,'both'); kolorimage = a;   
            end
%             if size(kolorimage,2) > 1024, a = kolorimage(:,(size(a,2)-1024)/2:(size(a,2)-1024)/2+1023); kolorimage = a;
%             elseif size(kolorimage,2) > 512, a = kolorimage(:,(size(kolorimage,2)-512)/2:(size(kolorimage,2)-512)/2+511); kolorimage = a;
          %  kolorimage = a;
                 %metoda = 0; 
            fprintf('Rozmiar obrazu skorygowano do potęgi dwójki\n'); 
            rescaling = 1; 
       end 
    end      
end
if klasa == 0
    appwalki = kolorimage;
elseif klasa == 1
% funkcja aproksymacji/ekstrakcji wałków 
    FOURIER = 0; DCT = 0; BDCT = 0; WAV = 0; CONT = 0; UWAV = 0; CURV = 0; CWT = 0; SHEAR = 0; SURF = 0;

    if metoda == 0, YY = kolorimage; fact = 1; %if PAR3==1, numth = 4; etyk = 2; end
    elseif metoda == 1, FOURIER = 1; if factx > 0, fact = factx; else fact=0.0025; end 
    elseif metoda == 2, DCT = 1; if factx > 0, fact = factx; else fact=0.001; end
    elseif metoda == 3, BDCT = 1; if factx > 0, fact = factx; else fact=0.99; end  
    elseif metoda == 4, WAV = 1; if factx > 0, fact = factx; else fact=0.001; end
    elseif metoda == 5, CONT = 1; if factx > 0, fact = factx; else fact=0.001; end
    elseif metoda == 6, UWAV= 1; if factx > 0, fact = factx; else fact=0.007; end
    elseif metoda == 7, CURV = 1; if factx > 0, fact = factx; else fact=0.005; end
    elseif metoda == 8, CWT = 1; if factx > 0, fact = factx; else fact=.0003; end 
    elseif metoda == 9, SHEAR = 1; if factx > 0, fact = factx; else fact=0.04; end 
    elseif metoda == 10, SURF = 1; if factx > 0, fact = factx; else fact=0.02; end 
  
    end
    n1 = size(kolorimage,1); n2 = size(kolorimage,2); N = n1*n2;
    numa = round(fact*N); hh=0; %liczba uwzglednionych wspolczynników okreslonej reprezentacji

    %for k=1:size(kolorimage,3)
    
    %%*********** aproksymacja z bazach *********
    if FOURIER %aproksymacja DCT
        A = fft2(kolorimage); a = A(:); 
        rozm = size(a,1); c = zeros(rozm,1);
        [fabsa,ind] = sort(abs(a),'descend');

        OMSORT = ind(1+hh:numa);
        c(OMSORT) = a(OMSORT); C = reshape(c,n1,n2);

        Y = real(ifft2(C)); 
%         figure, imshow(YY,[]); title('Fourier approximation');
%         disp('FourierAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N);
%         imwrite(uint8(YY),sprintf('atom_fur.pgm'),'pgm');
%         figure, imshow(histeq(uint8(YY)),[]); title('Hist Fourier approx'); 
    elseif DCT %aproksymacja DCT
        A = dct2(kolorimage); a = A(:); 
        rozm = size(a,1); c = zeros(rozm,1);
        [dabsa,ind] = sort(abs(a),'descend'); %nonlinear

        OMSORT = ind(1+hh:numa); %nonlinear
        c(OMSORT) = a(OMSORT); %nonlinear
        C = reshape(c,n1,n2); %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(kolorimage,YY)); disp(rozm/N);
%         %imwrite(uint8(YY),sprintf('atom_dct.pgm'),'pgm');
%         figure, imshow(histeq(uint8(YY)),[]); title('Hist DCT approx'); 
    elseif BDCT %aproksymacja DCT w blokach
        w = 32; %wielkość bloku
        fL = zeros(n1,n2);

        for i=1:n1/w
            for j=1:n2/w
                seli = (i-1)*w+1:i*w; selj = (j-1)*w+1:j*w;
                fL(seli,selj) = dct2( kolorimage(seli,selj) );
            end
        end
        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,n1,n2);
        Y = zeros(n1,n2);
        for i=1:n1/w 
            for j=1:n2/w
                seli = (i-1)*w+1:i*w; selj = (j-1)*w+1:j*w;
                Y(seli,selj) = idct2( C(seli,selj) );
            end
        end
%           figure, imshow(YY,[]); title('BDCT approximation');
%          disp('BDctAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N);
%          imwrite(uint8(YY),sprintf('atom_bdct.pgm'),'pgm');
%         figure, imshow(histeq(uint8(YY)),[]); title('Hist BDCT approx'); 
    elseif WAV
     %korekta rozmiaru obrazu  kolorimage = kolorimage(51:562,71:582); n1=512;n2=512;  
        %aproksymacja falkowa
        wkind = 'ort';

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

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

        if strcmp(wkind,'ort'), Y = IWT2_PO(C,3,qmf);
        else Y = IWT2_SBS(C,3,qmf,dqmf); end
%        figure, imshow(YY,[]); title('Wavelet approximation');
%         disp('WaveletAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N);
%         disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(kolorimage(70:136,167:250)))- var(var(kolorimage(70:136,17:100))))/(var(var(kolorimage(70:136,167:250)))+ var(var(kolorimage(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');
%             figure, imshow(histeq(uint8(YY)),[]); title('Hist Wave approx'); 
    elseif UWAV
        %aproksymacja falkowa undecimated
        wkind = 'sym2';%'sym1';%'coif2';%'bior4.4';
        levels = 4;
        if mod(n1,4)~=0, n1 = n1+mod(n1,4); end 
        if mod(n2,4)~=0, n2=n2+mod(n2,4); end
        A = swt2(kolorimage,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,n1,n2,1+3*levels);

        Y = iswt2(C,wkind);
%         figure, imshow(YY,[]); title('UnWavelet approximation');
%         disp('UnWaveletAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N);
    elseif CONT
        %aproksymacja conturletowa
        % Parameters
        pfilt = '9-7'; dfilt = 'pkva';
        nlevs = [2, 2, 3];%[0, 0, 0, 3, 4]; %[0, 0, 4, 4, 5]1
        %korekta rozmiaru obrazu
    %     if mod(n1,4)~=0, n1 = n1+mod(n1,4); end 
    %     if mod(n2,4)~=0, n2=n2+mod(n2,4); end
    %     kolorimage = imresize(kolorimage, [n1 n2]);
        A = pdfbdec(kolorimage, pfilt, dfilt, nlevs);
        [a,s] = pdfb2vec(A);

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

        C = vec2pdfb(c, s);
        Y = pdfbrec(C, pfilt, dfilt); 
%         figure, imshow(YY,[]); title('Contourlet approximation');
%         disp('ContourletAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N);
%         disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(kolorimage(70:136,167:250)))- var(var(kolorimage(70:136,17:100))))/(var(var(kolorimage(70:136,167:250)))+ var(var(kolorimage(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');
%         figure, imshow(histeq(uint8(YY)),[]); title('Hist CONT approx'); 
    elseif CURV
        %aproksymacja curveletowa
        C = fdct_wrapping(kolorimage,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}(:); ac = size(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);
        [cuabsa,ind] = sort(abs(a),'descend');
        OMSORT = ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
        c(OMSORT) = a(OMSORT);

        C{1}{1} = reshape(c(1:rozma),ac(1),ac(2)); 
        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)); 
%         figure, imshow(YY,[]); title('Curvelet approximation');
%         disp('CurveletAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N);
%         disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(kolorimage(70:136,167:250)))- var(var(kolorimage(70:136,17:100))))/(var(var(kolorimage(70:136,167:250)))+ var(var(kolorimage(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))))));
%         figure, imshow(histeq(uint8(YY)),[]); title('Hist CURV approx'); 
    elseif 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(kolorimage,levels,biort,qshift);
        %[Yl,Yh] = dtwavexfm2ap(kolorimage,4,biort,qshift); %znormalizowane filtry qshift

        a = Yl(:); %aa =a;%
        rozma = size(a,1); na = size(Yl);
        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);
        [cwtabsa,ind] = sort(abs(a),'descend'); %a  sortowanie po zwiększonych współczynnikach
        OMSORT = ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
        c(OMSORT) = a(OMSORT);

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

        Y = dtwaveifm2(Yl,Yh,biort,qshift);%,gm);
        %YY = (Y-min(min(Y)))/(max(max(Y))-min(min(Y)))*255;
%         figure, imshow(YY,[]); title('Complex Wavelet approximation');
%         disp('DT-CWTAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N*2); %2 bo zespolone
%         %disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(kolorimage(70:136,167:250)))- var(var(kolorimage(70:136,17:100))))/(var(var(kolorimage(70:136,167:250)))+ var(var(kolorimage(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');
%         figure, imshow(histeq(uint8(YY)),[]); title('Hist CW approx'); 
     %   figure, imshow(adapthisteq(uint8(YY)),[]); title('Adapt Hist CW approx'); 

    elseif SHEAR
        %aproksymacja shearlets
        [ST, Psi] = shearletTransformSpect(kolorimage);
        %[ST, Psi] = shear_trans(kolorimage);

        a = ST(:); 
        sz = size(ST);
        rozm = size(a,1); c = zeros(rozm,1);
        [shabsa,ind] = sort(abs(a),'descend');
        OMSORT = 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);
        %Y = inverse_shear(C,Psi);
%         figure, imshow(YY,[]); title('Shearlet approximation');
%         disp('ShearletAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N);
%         disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(kolorimage(70:136,167:250)))- var(var(kolorimage(70:136,17:100))))/(var(var(kolorimage(70:136,167:250)))+ var(var(kolorimage(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))))));
%             figure, imshow(histeq(uint8(YY)),[]); title('Hist sherlets approx'); 
    elseif 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;
                                            %Lev_array
        [G, Recinfo] = Surfdec(kolorimage, Pyr_mode, Lev_array, 'ritf', 'bo', bo);
        C = G;
       % licz = 0; %a=0;
        %a = C{1}{1}{1}(:); r =  size( C{1}{1}{1},1); r = [r;size( C{1}{1}{1},2)];
        a = C{4,1}(:); r =  size( C{4,1},1); r = [r; size(C{4,1},2)];
        C{4,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);
        [suabsa,ind] = sort(abs(a),'descend');
        OMSORT = ind(1:round(numa*rozm/N)); %OMSORT = ind(1:numa);
        c(OMSORT) = a(OMSORT);
         disp('liczba zer');
         nnz(c)
        %c = a; %dodatek
        rozx = r(1); rozy = r(2); kon =  rozx*rozy; 
        b = reshape(c(1:kon),rozx,rozy); 
        C{4,1}(:) = b;
        licz = 3; pocz = kon+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);
                    %max(max(b))
                    C{n}{dd}{sd} = b;
                    pocz =kon+1; licz = licz+2;
                end
            end
        end   

        Y = Surfrec(C, Recinfo); 
%         figure, imshow(YY,[]); title('Surface approximation');
%         disp('SurfaceletAPP='); disp(psnr(kolorimage,YY)); disp(rozm/N);
%         disp(sprintf('Wzmocnienie3): oryginal=%5.5f process=%5.5f', abs(var(var(kolorimage(70:136,167:250)))- var(var(kolorimage(70:136,17:100))))/(var(var(kolorimage(70:136,167:250)))+ var(var(kolorimage(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))))));
%         figure, imshow(histeq(uint8(YY)),[]); title('Hist sufrace approx'); 
    else
        Y = kolorimage;    
    end
    appwalki = rescale(Y,0,255); %clamp(Y,0,255);
elseif klasa == 2 %PCA
    I = double(kolorimage);
    %X = reshape(I,size(I,1)*size(I,2),3);
    X = reshape(I,size(I,1)*size(I,2),1);
    coeff = pca(X); coeff(abs(coeff) < 2*median(abs(coeff))) = 0; %2.5
    Itransformed = X.*coeff;
    Ipc1 = reshape(Itransformed,size(I,1),size(I,2));
%    Ipc1 = reshape(Itransformed(:,1),size(I,1),size(I,2));
%     Ipc2 = reshape(Itransformed(:,2),size(I,1),size(I,2));
%     Ipc3 = reshape(Itransformed(:,3),size(I,1),size(I,2));
%     figure, imshow(Ipc1,[]);
%     figure, imshow(Ipc2,[]);
%     figure, imshow(Ipc3,[]);
    appwalki = rescale(Ipc1,0,255); %clamp(Y,0,255);
elseif klasa == 3 %ICA
    ica
    appwalki = rescale(Ipc1,0,255); %clamp(Y,0,255);
elseif klasa == 4 %KSVD
    if metoda == 1, Y = ksvddenwal(kolorimage); 
    else
       Y = kolorimage; 
    end
    appwalki = rescale(Y,0,255); %clamp(Y,0,255);
elseif klasa == 5 %Split Bregman method for TV denoising 
    mu = 3;%150;
   Y = SB_ATV(double(kolorimage),mu);
   % Y = SB_ITV(double(kolorimage),mu);
    Y = reshape(Y,size(kolorimage,1),size(kolorimage,2));
    appwalki = rescale(Y,0,255); %clamp(Y,0,255);
end

if rescaling == 1
    % przywrócenie rozmiaru oryginalnego obrazka
    targetSize = [sizev sizeh]; r = centerCropWindow2d(size(appwalki),targetSize); a = imcrop(appwalki,r);
    appwalki = a;  %b = imcrop(kolorimage,rect);
elseif sum(size(kolorimageo) ~= size(appwalki)) > 0
    targetSize = size(kolorimageo); r = centerCropWindow2d(size(appwalki),targetSize); a = imcrop(appwalki,r);
    appwalki = a;    
end

ssimval = ssim(uint8(appwalki),uint8(kolorimageo));
disp(ssimval);