function [EllipseParamsAccurate] = StarOperator(Img,EllipseData,Star)

%Ucelem funkce je vysoce presny vypocet souradnic stredu elips v obraze,
%vstupem je cernobily obraz, data o elipsach a pocet radialnich paprsku. 
%Pro vylepseni reseni je mozne v budoucnu pouzit polarni souradnice misto 
%kartezskych.

%VYGENEROVANI PARSKU OKOLO OCEKAVANE HRANY

for i=1:size(EllipseData,1);                                                %pro vsechny elipticke objekty
    A = EllipseData(i,1);                                                   %nachystani hodnot pro nasledne rovnice - parametry elipsy, pozn. - u x0 a y0 je hodota posunuta tak, aby byl pocatek s.s. v levem hornim rohu snimku
    B = EllipseData(i,2);
    C = EllipseData(i,3);
    x0 = EllipseData(i,4)+0.5;
    y0 = EllipseData(i,5)+0.5;
    
    for j=0:Star-1;                                                         %vytvoreni paprsku                                                       
        Fi = pi/Star*j;                                                     %vytvoreni uhlu
        vx(i,j+1) = cos(Fi);                                                %smerove vektory
        vy(i,j+1) = sin(Fi);
        t = sqrt(1/(A*vx(i,j+1)^2+C*vy(i,j+1)^2+B*vx(i,j+1)*vy(i,j+1)));    %vypocet parametru t primky - urci hranici elipsy, vznika dosazenim parametricke rovnice primky do rovnice elipsy
                                
        x1 = x0+t*vx(i,j+1);                                                %vytvoreni pruseciku s elipsou - reseni se rozpada na dve poloviny - "kladnou" a "zapornou", kazda symbolizuje polovinu elipsy
        x2 = x0-t*vx(i,j+1);
        y1 = y0+t*vy(i,j+1);
        y2 = y0-t*vy(i,j+1);
               
        for k=-5:5                                                          %interpolace bodu na primce - 10 kroku radialne dovnitr i ven okolo predpokladane hrany pro obe poloviny elipsy
            
            sx1(i,j+1) = x1-vx(i,j+1)*k;                                    %vytvoreni souradnic interpolovynych bodu
            sy1(i,j+1) = y1-k*vy(i,j+1);
            sx2(i,j+1) = x2+vx(i,j+1)*k;                                              
            sy2(i,j+1) = y2+k*vy(i,j+1);
            
            interpol1(j+1,k+6,i) = BiCubicConvolution(Img, sx1(i,j+1)+0.5, sy1(i,j+1)+0.5);       
                                                                            %interpolace hodnot jasu pixelu na bodech na paprscich
            
            interpol2(j+1,k+6,i) = BiCubicConvolution(Img, sx2(i,j+1)+0.5, sy2(i,j+1)+0.5);
            
             if k==-5                                                       %ulozeni pocatku paprsku
                sxs1(i,j+1) = sx1(i,j+1);
                sys1(i,j+1) = sy1(i,j+1);
                sxs2(i,j+1) = sx2(i,j+1);
                sys2(i,j+1) = sy2(i,j+1);
             end
        end
    end
end

%DETEKCE HRAN - MOMENT PRESERVATION

for i=1:size(EllipseData,1)                                                                                                                                     
    for j=1:Star
        Index1(i,j) = MomentPreserv(interpol1(j,:,i));                      %V tomto indexu (subindexu) vektoru je nalezena hrana 
        Index2(i,j) = MomentPreserv(interpol2(j,:,i));
        hx1(i,j)=sxs1(i,j)-Index1(i,j)*vx(i,j);                             %Souradnice bodu lezicich na hranach elipsy
        hx2(i,j)=sxs2(i,j)+Index2(i,j)*vx(i,j);
        hy1(i,j)=sys1(i,j)-Index1(i,j)*vy(i,j);
        hy2(i,j)=sys2(i,j)+Index2(i,j)*vy(i,j);      
    end
end

hx=horzcat(hx1,hx2);                                                        %Slouceni obou polovin reseni dohromady - vystup detekce - body lezici na hranach elipsy
hy=horzcat(hy1,hy2);

%FITTING ELIPS DO DAT ZISKANYCH DETEKCI HRAN

for i=1:size(hx,1)
    
    a0 = [EllipseData(i,1) EllipseData(i,2) EllipseData(i,3) EllipseData(i,4) EllipseData(i,5)];                    
                                                                            %Startovni hodnoty
    f = @(a) (a(1).*hx(i,:).^2+(a(2).*hx(i,:).*hy(i,:))+(a(3).*((hy(i,:).^2)))-((2.*a(1).*a(4)+a(2).*a(5)).*(hx(i,:)))-((2.*a(3).*a(5)+a(2).*a(4)).*hy(i,:))+(a(1).*(a(4).^2))+(a(2).*a(4).*a(5))+(a(3).*(a(5).^2))-1);
                                                                            %Zavedeni rovnice elipsy
    options = optimset('Display','off');                                    %Nastaveni iteracniho reseni - default
                                                                            
    [EllipseParamsAccurate(i,:)] = lsqnonlin(f, a0, [], [], options);       %Reseni metodou nelinearnich nejmensich ctvercu   

end

end

