function [index, stred] = spine_curve(gabor_spine_img, okno)
% funkce, s pomoci ktere lze ziskat levou a pravou hranici patere v
% sagitalnim MRI snimku

% Nasatvení proměnných a konstant
stred = round(mean(1:size(gabor_spine_img, 1)));
index = zeros(size(gabor_spine_img, 1), 2);
h = sum(gabor_spine_img(stred-okno:stred+okno, :), 1);
% krivka = cumsum(h); % samotnou krivku nepotrebuji pro vypocet

% Nalezení počátečních bodů
if mean(h) < 0.03 % heuristika, zajisti lepsi detekci okrajovych bodu
    min_peak_height = max(h)/25;
else
    min_peak_height = max(h)/6;
end

[~, temp] = findpeaks(h, 'MinPeakHeight', min_peak_height, 'MinPeakDistance', 40, 'NPeaks', 2);
index(stred, :) = temp;

sirka = temp(2) - temp(1);

% Výpočet křivky od středu k počátku
index = extend_curve(gabor_spine_img, 1, stred, okno, index, sirka);
% Výpočet křivky od středu ke konci
index = extend_curve(gabor_spine_img, -1, stred, okno, index, sirka);

end

function index = extend_curve(gabor_spine_img, up_down, stred, okno, index, sirka)
    % up_down is either 1 (going up the image) or -1 (going down the image)
    
    assert(up_down == 1 || up_down == -1)
    
    if up_down == 1
        indices = stred-1:-1:1;
    else
        indices = stred+1:size(gabor_spine_img, 1);
    end
    
    for i = indices
        start = i - okno;
        if start < 1
            start = 1;
        end
        last = i+ okno;
        if last > size(gabor_spine_img, 1)
            last = size(gabor_spine_img, 1);
        end
        
        h = sum(gabor_spine_img(start:last, :), 1);
        
%         krivka = cumsum(h);

        [~, temp] = findpeaks(h, 'MinPeakHeight', max(h)/50);
        [~, closest_left_border] = min(abs(temp-index(i+up_down, 1)));
        index(i, 1) = temp(closest_left_border);
        [~, closest_right_border] = min(abs(temp-index(i+up_down, 2)));
        index(i, 2) = temp(closest_right_border);
        
        if abs((index(i, 1) - index(i, 2))) < sirka/10 % pokud je kanal patere moc uzky
            if (index(i, 1) - index(i+up_down, 1)) < (index(i, 2) - index(i+up_down, 2)) % je levy nebo pravy okraj na spravnem miste?
                index(i, 2) = index(i, 2) + sirka; % pokud levy, posunu pravy okraj vice doprava
            else
                index(i, 1) = index(i, 1) - sirka; % pokud je to naopak, pak udelam opak
            end
        end
        
%         plot(h)
%         hold on;
%         plot(temp, h(temp), '*y')
%         plot(temp([closest_left_border, closest_right_border]), h(temp([closest_left_border, closest_right_border])), '*r');
%         hold off
%         pause
    end
    
end

