function [xs,ys,radius,xp,yp,im] = radial_symmetry (im)
%% metoda rychle radialni symetrie
% skript prevzat od ..., ktery jsem dostala k dispozici od doc. Kolare
% im - vstupni sedotonovy snimek
% xs,ys  - souradnice stredu detekovane kruznice
% radius - polomer detekovane kruznice
% xp, yp - souradnice detekovane kruznice
%--------------------------------------------------------------------------
    % vektror prohledavanych polomeru duhovky
    radii = 20:65; % 
    alpha = 3;% 3 nejvice prisny z 1-3
    % velikost snimku
    [rows, cols] = size(im);

    % Gradient snimku - pouziti gradientni funkce
    [imgx,imgy] = gradient (double(im),2); % 2 je interval mezi body v kazdem smeru 
    norm = sqrt(imgx.^2 + imgy.^2)+eps; % eps vrátí vzdálenost od 1,0 do dalšího největší cisla s dvojitou přesností
    imgx = imgx./norm;   
    imgy = imgy./norm;

    S = zeros(rows,cols);  % nadefinovani symetricke matice

    [x,y] = meshgrid(1:cols, 1:rows); % transformuje doménu specifikovanou vektory x a y do polí X a Y, které mohou být použity k vyhodnocení funkce dvou proměnných a trojrozměrných ok sítě / povrch plots.
    
    MIN = zeros (1,length(radii)); % promenna pro jednotlive polomery
    r = zeros(length(radii),2);
    k = 1; % preddefinovani promenne pro ukladani matice S

    %% postupne prochazeni vsech polomeru
    for n = radii
        M = zeros(rows,cols);  % matice mohutnosti, inicializace nulová
        O = zeros(rows,cols);  % matice orientace

        % Souradnice pozitivne/negativne ovlivnenych pixelu
        posx = x + round(n*imgx); 
        posy = y + round(n*imgy);
        
        negx = x - round(n*imgx);
        negy = y - round(n*imgy);
        
        % Ummisteni vsech dotcenych bodu do matice
        % zaporne prevedeny na 1
        % kladne bud zustanou, nebo upraveny podle poctu radku/sloupcu
        % snimku
        posx( posx<1 )    = 1;    % postx - pozitivne ovlivnene pixely
        posx( posx>cols ) = cols; % cols - celkovy pocet sloupcu puvodniho snimku
        posy( posy<1 )    = 1;
        posy( posy>rows ) = rows;  % rows - celkovy pocet radku puvodniho snimku
        
        negx( negx<1 )    = 1;     % negx - negativne ovlivnene pixely
        negx( negx>cols ) = cols;
        negy( negy<1 )    = 1;
        negy( negy>rows ) = rows;
        
         % indexy k logickemu indexovani 
         pos = sub2ind([rows cols], posy, posx); % příkaz určuje jeden index odpovídající sadě hodnot dolnich indexu.
         neg = sub2ind([rows cols], negy, negx); 

         % matrix to vector
         pos = pos (:);
         neg = neg (:);


        % projekce orientace & mohutnosti
        % narustani matic pro vsechny pozitivni a negaticni pixely
        O(:) = full(sparse([pos,neg], 1, [ones(rows,cols),-ones(rows,cols)])); % sparse - řidká matice
        % sparse - vytesnuje nulove prvky 
        M(:) = full(sparse([pos,neg], 1, [norm,-norm]));

        % detekce pixelu, ktere jsou mimo vstupni snimek
        O(1,:) = 0;
        O(rows,:) = 0;
        O(:,1) = 0;
        O(:,cols) = 0;

        M(1,:) = 0;
        M(rows,:) = 0;
        M(:,1) = 0;
        M(:,cols) = 0;
     
    if max([pos,neg]) ~= length(pos)
        
        O((max([pos,neg])+1):length(pos)) = 0;
        M((max([pos,neg])+1):length(pos)) = 0;
        
    end
        %         
        % Set the orientation matrix into the of -kappa +kappa according to
        % the article suggestion

        if n == 1 
           kappa = 8; 
        
        else
           kappa = 9.9;
        end
        
        O(O >  kappa) =  kappa;  
        O(O < -kappa) = -kappa;  
       
        % vypocet F
        F = M./kappa .* (abs(O)/kappa).^alpha;
        
        % vytvoreni vazeneho gaussian kernel A
        A = fspecial('gaussian',[n n], 0.25*n) * n; % n = radii  
        
        % konvoluce, ktera vypocita symetrickou matici S
        S1 = imfilter(F,A); % konvoluce mezi F a A
        S(:,:,k) = S1;
        [C,I] = min(S1(:)); % hledani minima - protoze je tmava duhovka       
        MIN(k) = C; % hodnota minima ze S1 pro dany polomer
        [r(k,1) r(k,2)] = ind2sub(size(S1),I); % souradnice maxima z S1
        k = k + 1;
    end  

%% Zobrazeni
% nalezeni minima ze symetricke mapy a lokalizace stredu a vnitrni hranice
% duhovky vstupniho obraxu

    [C,I] = min (MIN); % najiti minima ze vsech polomeru, min protoze je tmava duhovka
    S2 = S (:,:,I); % vyber ze symetricke matice, I je pozice minima
    radius = radii(I); % vysledny polomer, I-ta pozice z radii

    [C,I] = min (S2(:)); % zjisteni indexu minima
    [ys xs] = ind2sub(size(S2),I); % zjisteni souradnic stredu 

    Fkrok = 0.5; % krok pro stupne
    % 
    i = 0 : Fkrok : 360; %2*pi+0.1;
    xp = radius*cos((i*pi)/180); % souradnice vnitriho polomeru
    yp = radius*sin((i*pi)/180);



end