function [output_name, statistiky, elapset_time] = stabil_CRS(filename,pocetjedincu,procento_obrazu,uhel_rotace,maxit,delta_nastav,hranice_stabilizace,kriterium_kvality)
%% ZADÁNÍ:
% Stabilizace videa využívající globální optimalizační algoritmy
%
% Naprogramujte aplikaci pro stabilizaci jednoduché videosekvence
% využívající optimalizační algoritmus CRS (Controlled Random Search).
%% NASTAVENÍ STABILIACE OBRAZU
% kriterium_kvality = 0.995; %kritérium kvality
% 
% maxit = 300; % max počet iterací
% delta_nastav = 0.0005; % rozdíl v kos. kritériu mezi jednotlivými iteracimi stabilizace
% %r = 0;  % pocet nejhorsich snimku, ktere se maji vynechat
% hranice_stabilizace = 1; % Hranice, pri ktere se ma obraz stabilizovat, mira shody predchoziho obrazku s novym, pri = 1 se bude stabilizovat vzdy.
% %% NASTAVENÍ CRS
% procento_obrazu = 0.025; % delta x a delta y
% uhel_rotace = 4; % [°]
% pocetjedincu = 40; % rozmer nahodnych cisel
Dimenze = 3; % musi byt 3, x, y, theta, zkos videa neuvazujeme
%% NAČTENÍ REFERENČNÍHO SNÍMKU
tic;
videoreader = VideoReader(filename); %kvuli parametrum snimky,vidH a vidW
srovnani = readFrame(videoreader);

video = vision.VideoFileReader(filename,'ImageColorSpace',... %kvuli pouziti fce step
    'intensity','VideoOutputDataType','double');


snimky = ceil(videoreader.Duration * videoreader.FrameRate);
progress_bar = waitbar(0,'Prosím čekejte');
%nacteni poctu pozadovanych snimku z videa
%pocet nactenych snimku v zavislosti na zvolene delce videa

%urceni rozmeru vystupniho obrazu
v_height_out = videoreader.Height;
v_width_out = videoreader.Width;
dx = round(v_height_out*procento_obrazu/100);
dy = round(v_width_out*procento_obrazu/100);

% vhodna prealokace proměnné movie pro spravne spusteni prikazem implay
movie(1:snimky) = struct('cdata', zeros(v_height_out, v_width_out),'colormap', []);

%referencni snimek - nacteni prvniho snimku
movie(1).cdata = step(video);
I1_in = movie(1).cdata;
I1n = I1_in(:);

%urceni rozmeru vstupniho obrazu
v_height_in = size(I1_in,1);
v_width_in = size(I1_in,2);

%generování náhodné populace pro CRS optimalizaci
uhel_rotace = deg2rad(uhel_rotace);
TransPopul = rand(Dimenze,pocetjedincu);
TransPopul(1,:) = -dy + (dy - -dy).*TransPopul(1,:); %pro posun x a y
TransPopul(2,:) = -dx + (dx - -dx).*TransPopul(2,:); %pro posun x a y
TransPopul(3,:) = -uhel_rotace + (uhel_rotace - -uhel_rotace).*TransPopul(3,:); %pro uhel

Ck_vektor_smycky(1) = 0.1;
vektor_posunu_x = zeros(1,snimky);
vektor_posunu_y = zeros(1,snimky);
Ck_porov_vektor = zeros(1,snimky);
Ck_vektor_transformovanych = zeros(1,snimky);
statistiky = zeros(3,snimky);
%% smycka pro vykonani stabilizace a vepsani do movie(k)
for k = 2:snimky
    
    movie(k).cdata = step(video);
    I2_in = movie(k).cdata;
    I2n = I2_in(:);
    
    %cosinove kriterium (bez registrace - pro porovnani s konecnym kriteriem)
    Ck_porov = (sum(I1n.*I2n))./(sqrt(sum(I1n.*I1n)).*sqrt(sum(I2n.*I2n)));
    
    Ck_porov_vektor(k) = Ck_porov;
    
    if Ck_porov < hranice_stabilizace %podminka pro rozhodnuti vykonani stabilizace
        Ck_v = zeros(1,pocetjedincu);
        %ohodnoceni kazdeho clena populace cosinovym kriteriem ve smycce
        for i = 1:pocetjedincu
            
            %vytvoreni transformacni matice R
            R = maketform('affine',[0+cos(TransPopul(3,i)) 0+sin(TransPopul(3,i)) 0;...
                0-sin(TransPopul(3,i)) 0+cos(TransPopul(3,i)) 0; TransPopul(1,i) TransPopul(2,i) 1]);
            
            %transformovany pohybovy snimek
            I3 = imtransform(I2_in,R,'XData',[1 v_width_in],'YData',[1 v_height_in]);
            %                 imshow(I3);
            %                 shg;
            
            %prevedeni I3 na vektor
            I3n = I3(:);
            bI3 = I3n~=0;
            
            %cosinove kriterium jako vysledny vektor vypoctenych kriterii ve smycce
            Ck_v(i) = (sum(I1n(bI3).*I3n(bI3)))./(sqrt(sum(I1n(bI3).*I1n(bI3))).*sqrt(sum(I3n(bI3).*I3n(bI3))));
        end
               
        delta_Ck = 1; %pocatecni deklarace aby mohlo while probehnout
        iterace = 0; %pocatecni deklarace aby mohlo while probehnout
        
        %iteracni smycka
        while ((max(Ck_v)<kriterium_kvality) || (delta_nastav < delta_Ck) || (delta_Ck == 0)) && (maxit>iterace)
            %zasadni podminka pro vykonavani CRS optimalizace,parametry lze nastavovat
            statistiky(1,k) = statistiky(1,k) + 1;
            [~,c] = min(Ck_v); %nalezeni minima a ulozeni jeho pozice do c
            
            Ck_gen = 0; %pocatecni deklarace aby mohlo while probehnout
            iterace = iterace+1; %inkrementace promenne iterace
            
            %generovani simplexu (viz CRS optimalizace)
            while Ck_gen < Ck_v(c)
                
                SimplexIndex = randperm(pocetjedincu);
                SimplexIndex = SimplexIndex(1:Dimenze+1);
                
                Simplex = TransPopul(:,SimplexIndex);
                Simplex_for_worst = TransPopul(:,SimplexIndex);
                Ck_vSimplex = Ck_v(:,SimplexIndex);
                
                %nalezeni a vyjmuti nejhorsiho bodu ze Simplexu
                [~,e] = min(Ck_vSimplex);
                Simplex(:,e) = [];
                
                worst = Simplex_for_worst(:,e);
                suma = (Simplex(:,1)+Simplex(:,2)+Simplex(:,3));
                
                %vypocet teziste
                g = (suma)./Dimenze;
                
                %generovany novy bod
                Xr = 2.*g-worst;
                
                while -procento_obrazu >=Xr(1,1)>= procento_obrazu && -procento_obrazu >=Xr(2,1)>= procento_obrazu && -uhel_rotace >=Xr(3,1)>= uhel_rotace
                    %pokud se novy bod nachazi mimo zvoleny prostor, generuje se znovu
                    
                    SimplexIndex = randperm(pocetjedincu);
                    SimplexIndex = SimplexIndex(1:Dimenze+1);
                    
                    Simplex = TransPopul(:,SimplexIndex);
                    Simplex_for_worst = TransPopul(:,SimplexIndex);
                    Ck_vSimplex = Ck_v(:,SimplexIndex);
                    
                    %nalezeni a vyjmuti nejhorsiho bodu ze Simplexu
                    [~,e] = min(Ck_vSimplex);
                    Simplex(:,e) = [];
                    
                    worst = Simplex_for_worst(:,e);
                    suma = (Simplex(:,1)+Simplex(:,2)+Simplex(:,3));
                    
                    %vypocet teziste
                    g = (suma)./Dimenze;
                    
                    %generovany novy bod
                    Xr = 2.*g-worst;
                end;
                
                R = maketform('affine',[0+cos(Xr(3,:)) 0+sin(Xr(3,:)) 0;...
                    0-sin(Xr(3,:)) 0+cos(Xr(3,:)) 0; Xr(1,:) Xr(2,:) 1]);
                
                I4 = imtransform(I2_in,R,'XData',[1 v_width_in],'YData',[1 v_height_in]);
                %                     imshow(I4);
                
                I4n = I4(:);
                bI4 = I4n~=0;
                
                %cosinove kriterium - kvalita noveho bodu
                Ck_gen = (sum(I1n(bI4).*I4n(bI4)))./(sqrt(sum(I1n(bI4).*I1n(bI4))).*sqrt(sum(I4n(bI4).*I4n(bI4))));
            end
                        
            TransPopul(:,c) = Xr; %nahrazeni nejhorsiho bodu populace bodem novym
            
            Ck_v(c)= Ck_gen; %nahrazeni kriteria
            
            Ck_vektor_smycky = [Ck_vektor_smycky max(Ck_v)]; %ukladani jednotlivych nejlepsich kriterii
            delta_Ck = Ck_vektor_smycky(end)-Ck_vektor_smycky(end-1);
            
        end
        
        Ck_vektor_transformovanych(k) = max(Ck_v); %finalni kriteria jednotlivych snimku
        
        %provedeni transformace
        [~,h]=max(Ck_v);
        Xl = TransPopul(:,h);
        R = maketform('affine',[0+cos(Xl(3,:)) 0+sin(Xl(3,:)) 0; 0-sin(Xl(3,:)) 0+cos(Xl(3,:)) 0; Xl(1,:) Xl(2,:) 1]);
        
        I4 = imtransform(movie(k).cdata,R,'XData',[1 v_width_out],'YData',[1 v_height_out]);
        %             imshow(I4);
        
        %nacteni transformovaneho obrazku do movie(k)
        movie(k).cdata = I4;
        A = double(R.tdata.T);
        [statistiky(2,k), statistiky(3,k)] =...
            mutualInfo(srovnani,imtransform(readFrame(videoreader),R,'XData',[1 v_width_out],'YData',[1 v_height_out]));
        
        %zjisteni maximalnich posunu pixelu pro orez obrazu
        vektor_posunu_x(k) = A(3,1);
        vektor_posunu_y(k) = A(3,2);
    else
        Ck_vektor_transformovanych(k) = Ck_porov;
        [statistiky(2,k), statistiky(3,k)] =...
            mutualInfo(srovnani,readFrame(videoreader));
        vektor_posunu_x(k) = 0;
        vektor_posunu_y(k) = 0;
    end
    waitbar(k/ snimky);
end
%% orez cernych pruhu
for l = 1:(snimky)
    orez_x = max(ceil(abs(vektor_posunu_x)));
    orez_y = max(ceil(abs(vektor_posunu_y)));
    if(orez_x == 0 && orez_y == 0)
        break;
    end
    if(orez_x == 0 )
        orez_x = 1;
    end
    if(orez_y == 0 )
        orez_y = 1;
    end
    movie(l).cdata = movie(l).cdata(orez_y:end-orez_y,orez_x:end-orez_x);
end
%% ulozeni stabilizovaneho videa
elapset_time = toc;
filename = strsplit(filename,'\');
output_name = [...
    'CRS_', '_PctJed',num2str(pocetjedincu),...
    '_KvKr',num2str(kriterium_kvality),...
    '_MaxIt',num2str(maxit),...
    '_Cas',num2str(elapset_time),...
    '_', filename{end}];
output = VideoWriter(output_name,'Uncompressed AVI');
open(output);

for m = 1:(snimky)
    writeVideo(output,movie(m).cdata);
end
close(progress_bar);
close(output);
end
