function [threshold,curve,curvefit] = cortical_bones_thresh(histogram,type,optim_method)
%Funkce pro zjištění globálního prahu pro segmentaci kortikálních částí
%kostí. Vstupem je histogram dat, typ pravděpobnostního rozdělení -
%Gaussovo nebo Pearsonovo typu VII a také typ optimalizační metody.

switch type   %vybrané rozdělení
    case 'gauss' %normalní (Gaussovo) rozdělení
        mygauss = @(x,xdata) x(1)*exp(-((xdata-x(2)).^2/(2*x(3).^2))); %Gaussova funkce
        % x(1) = výška píku (a)
        % x(2) = střední hodnota (u)
        % x(3) = směrodatná odchylka (sigma)
    
        [a u] = max(histogram); %amplituda a střední hodnota
        u = (u-1)/1023; % -1 index kvuli matlabovské indexaci od 1
        sigma = 0.015; %std počáteční sigma
        x = (0:1023)/1023; %osa x
        curve = mygauss([a u sigma],x); %počáteční odhad funkce
        
        switch optim_method %vybraná metoda optimalizace
            case 'exhaustive_search' %úplné prohledávání
                N = length(x); %délka osy x
                sig = Inf; %počáteční kriteriální hodnota
                s = 0.0001:0.0001:sigma; %vektor prohledávaných hodnot
                for i=1:length(s)
                    pom = mygauss([a u s(i)],x); %výpočet funkce
                    krit = sum((histogram'-pom).^2)/N; %výpočet kriteriální hodnoty
                    if krit<sig %pokud je kriteriální hodnota menši než doposud nejmenši
                       sig = krit; %přepsání
                       opt_s = s(i); %uložení hodnoty sigma
                    end
                end
                curvefit = mygauss([a u opt_s],x);
                threshold = u+12*opt_s; %stanovení prahu
                
            case 'random_search' %náhodné prohledávání
                N = 1000; %počet cyklů prohledávání
                sig = Inf; %počáteční kriteriální hodnota
                for i=1:N
                    s_rand = 0.0001+rand()*(sigma-0.0001); %náhodně vygenerované číslo v rozsahu <0.0001,sigma>
                    pom = mygauss([a u s_rand],x); %výpočet funkce
                    krit = sum((histogram'-pom).^2)/N; %výpočet kriteriální hodnoty
                    if krit<sig %pokud je kriteriální hodnota menši než doposud nejmenši
                        sig = krit; %přepsání
                        opt_s = s_rand; %uložení hodnoty sigma
                    end
                end
                curvefit = mygauss([a u opt_s],x); %výpočet optimalizované funkce
                threshold = u+12*opt_s; %stanovení prahu
                
            case 'lsq' %metoda nejmenších čverců
                lsq_options = optimset('lsqcurvefit'); %nastavení metody LSQ
                lsq_options.MaxFunEvals = 100; %max počet iterací
                lsq_options.TolFun = 1e-12; %max tolerance
                lb = [a-1e-6 u-1e-6 0.004]; %spodní hranice
                ub =[a+1e-6 u+1e-6 sigma+1e-6]; %horní hranice
                y = lsqcurvefit(mygauss,[a u sigma],curve,histogram',lb,ub,lsq_options); %řešení problému LSQ
                curvefit = mygauss(y,x); %výpočet optimalizované funkce
                threshold = u+12*y(3); %stanovení prahu
        end
    
    case 'pearson' %Pearsonovo rozdělění typu VII
        mypearson = @(x,xdata) x(1)*(1+(xdata-x(2)).^2/(x(3)*x(4).^2)).^-x(3); %Pearson VII function
        % x(1) = výška píku (c)
        % x(2) = centroid (u)
        % x(3) = koeficient špičatosti (m)
        % x(4) = šířka píku (a)
    
        [c u] = max(histogram);
        u = (u-1)/1023; % -1 index kvůli matlabovské indexaci od 1
        a = 0.015;
        m = 3;
        x = (0:1023)/1023;
        curve = mypearson([c u m a],x);
        
        switch optim_method %vybraná metoda optimalizace
            case 'exhaustive_search' %úplné prohledávání
                N = length(x); %délka osy x
                krit_best = Inf; %počáteční kriteriální hodnota
                a_vec = 0.0001:0.0001:a; %vektor prohledávaných hodnot a
                m_vec = m:0.1:7; %vektor prohledávaných hodnot m
                for i=1:length(a_vec)
                    for j=1:length(m_vec)
                        pom = mypearson([c u m_vec(j) a_vec(i)],x); %výpočet funkce
                        krit = sum((histogram'-pom).^2)/N; %výpočet kriteriální hodnoty
                        if krit<krit_best %pokud je kriteriální hodnota menší než doposud nejmenši
                            krit_best = krit; %přepsání
                            opt_a = a_vec(i); %uložení a
                            opt_m = m_vec(j); %uložení m
                        end
                    end
                end
                curvefit = mypearson([c u opt_m opt_a],x); %výpočet funkce
                threshold = u+10*opt_a; %stanovení prahu
                
            case 'random_search' %náhodné prohledávaní
                N = 1000; %počet cyklů prohledávání
                krit_best = Inf; %počáteční kriteriálni hodnota
                for i=1:N
                    a_rand = 0.0001+rand()*(a-0.0001); %náhodně generované a v rozsahu <0.0001,a>
                    m_rand = 0.1*randi([30 70]); %nahodně generované m v rozsahu <3,7>
                    pom = mypearson([c u m_rand a_rand],x); %výpočet funkce
                    krit = sum((histogram'-pom).^2)/N; %výpočet kriteriální hodnoty
                    if krit<krit_best %pokud je kriteriální hodnota menši než doposud nejmenši
                        krit_best = krit; %přepsání
                        opt_a = a_rand; %uložení a
                        opt_m = m_rand; %uložení m
                    end
                end
                curvefit = mypearson([c u opt_m opt_a],x); %výpočet funkce
                threshold = u+10*opt_a; %stanovení prahu
                
            case 'lsq' %metoda menjemších čtverců
                lsq_options = optimset('lsqcurvefit'); %LSQ nastavení
                lsq_options.MaxFunEvals = 100; %max počet iterací
                lsq_options.TolFun = 0.15; %max tolerance
                lb = [c-1e-6 u-1e-6 m-2 a-0.012]; %spodní hranice
                ub =[c+1e-6 u+1e-6 m+7 a+1e-6]; %horní hranice
                y = lsqcurvefit(mypearson,[c u m a],curve,histogram',lb,ub,lsq_options); %vyřešení LSQ problému
                curvefit = mypearson(y,x); %výpočet funkce
                threshold = u+10*y(4); %stanovení prahu
        end
end
end

