%Funkce pro optimální nalezení parametrů geometrické transformace rojem           %
% včel BA                                                                         %
%vystup - parametry transformace, kvalita - hodnota kriteriální funkce            %
%vstup - obraz1, obraz2, [parametry modelu]                                       %
%příklad spuštění:[vystup,kvalita] = BA_MSE(ref_obraz,transf_obraz,20,10,5,3,1,20)%
% ______________________________________________________________________          %

function [vystup,kvalita] = BA_MSE(ref_obraz,transf_obraz,n_pruzkumnic,m_best_pruzkumnic,elita,vyck,ngh,maxit)
%% výběr 2 obrázků 
obr1=ref_obraz;
obr2=transf_obraz;
%okno=[100,200,90,330];
% výběr obdélníkového okna
imshow(ref_obraz)
title('Tažením myši vyberte obdélníkovou oblast pro registraci')
h=imrect()
pos = h.getPosition();
pos_rect = round(pos);
okno=[pos_rect(2),pos_rect(2)+pos_rect(4),pos_rect(1),pos_rect(1)+pos_rect(3)];
%% zadana uloha
fcn=@(x,y) immse(x,y);  %definice kriteriální fce
xmin=-10;ymin=-10;zmin=-10;  %hranice
xmax=10;ymax=10;zmax=10;
%% inicializace
n=n_pruzkumnic;  %počet včel průzkumnic 
m=m_best_pruzkumnic;   %počet nejkvalitnějších průzkumnic 
e=elita;    % počet elitních průzkumnic
% vyck=3; 
% ngh=2;  
% maxit=100;
dim=3;
min_xy=[-10, -10]; % minima pro x a y
max_xy=[10, 10];
%% pocatecni populace
populace = cell(1,m);
populace2=zeros(1,(dim));
pop=zeros(1,length(dim));
rychlost = cell(1,m);
krit=zeros(1,n); % hodnota kvality

for i=1:n
    for j=1:dim
pop(j)=min_xy(1,1)+rand().*(max_xy(1,1)-min_xy(1,1));

    end
    populace{i} = pop;
    populace2(i,:)=pop;
    
end
%% samotny algoritmus
nove_vcely={};
ite=1;
h = waitbar(0,'Please wait...');
while ite<=maxit
    % 1. faze - vyber vcel
    for j=1:n
        
        pom=imtranslate(obr2,[populace2(j,1) populace2(j,2)]);
        pom=imrotate(pom,populace2(j,3),'bilinear','crop');
        krit(j)=-fcn(obr1(okno(1):okno(2),okno(3):okno(4)),pom(okno(1):okno(2),okno(3):okno(4)));
    end
   
    [kvalita_pred,indexy]=sort(krit,'descend');
    populace_ser=populace2(indexy,:);
    pracovni=populace_ser(1:m,:); % pouze stanoveny pocet nejlepsich vcel
    
        for i=1:m
            if i<=e
                nah_vcel=2*vyck;
            else
              nah_vcel=vyck;
            end
        
     
        delta=2*ngh*(rand(nah_vcel,3)-0.5);
        nove_vcel=pracovni(i,:)+delta;
        
            krit2=zeros(nah_vcel,1);
            for k=1:nah_vcel 
                
                
                pom2=imtranslate(obr2,[nove_vcel(k,1) nove_vcel(k,2)]);
                pom2=imrotate(pom2,nove_vcel(k,3),'bilinear','crop');
                krit2(k)=-fcn(obr1(okno(1):okno(2),okno(3):okno(4)),pom2(okno(1):okno(2),okno(3):okno(4)));
                
            end
         
            if any(krit2>kvalita_pred(i))
                [maxi, pos]=max(krit2);
                pracovni(i,:)=nove_vcel(pos,:);       
                 
            end
         end
     
     populace2(1:m,:)=pracovni;
       
    populace2(m+1:n,:)=[xmin,ymin,zmin]+rand(n-m,3).*([xmax,ymax,zmax]-[xmin,ymin,zmin]);
   for j=1:n
        
        pom3=imtranslate(obr2,[ populace2(j,1) populace2(j,2)]);
        pom3=imrotate(pom3,populace2(j,3),'bilinear','crop');
        krit3(j)=-fcn(obr1(okno(1):okno(2),okno(3):okno(4)),pom3(okno(1):okno(2),okno(3):okno(4)));
        
    end
       kv_mean(ite)=mean(krit3);
      [kv_max(ite),pos]=max(krit3);
    
    
    ite=ite+1;
    pause(.1);
    
    waitbar(ite / maxit)
end
close(h)
vystup=populace2(pos,:);
kvalita=-kv_max(ite-1)  ;      
end