function plot_prcurve(res_folder, gt_folder)

% folders must end with /

close all

% parameters

%d = 32.573;         % a detection within d pixels of the ground truth is considered correct

d = 20.358;

%d = 10;

threshold = 0.1:0.01:1;

precision = zeros(1,size(threshold,2));

recall = zeros(1,size(threshold,2));

f1 = zeros(1,size(threshold,2));

res_listing = dir([res_folder '*.csv']);

g = waitbar(0);

for h=1:size(threshold,2)

    true_pos = 0;

    false_neg = 0;

    total_pos = 0;

    for i=1:size(res_listing)

        results = csvread([res_folder res_listing(i).name]);

        % cut results off at the threshold

        for j=1:size(results,1)

            if results(j,3) < threshold(h)

                results = results(1:j-1,:);

                break

            end

        end

        total_pos = total_pos + size(results,1);

        [slide, remain] = strtok(res_listing(i).name, '_');

        subfolder = [gt_folder slide '_v2/'];

        ground_truth = csvread([subfolder res_listing(i).name]);

        % find the centroids

        len_gt = size(ground_truth,1);

        tmp_gt = zeros(len_gt,2);

        for j=1:len_gt;

            tmp = ground_truth(j,:);

            tmp = tmp(tmp~=0);

            len_tmp = size(tmp,2);

            sum_X = 0;

            sum_Y = 0;

            for k=1:len_tmp/2

                sum_X = sum_X + tmp(2*k-1);

                sum_Y = sum_Y + tmp(2*k);

            end

            mean_X = 2*sum_X/len_tmp;

            mean_Y = 2*sum_Y/len_tmp;

            tmp_gt(j,:) = [mean_Y mean_X];

        end

        ground_truth = round(tmp_gt);

        len_res = size(results,1);

        res_tally = zeros(len_res,1);

        for j=1:len_gt

            false_neg = false_neg + 1;

            for k=1:len_res

                if (ground_truth(j,1) - results(k,1))^2 + ...

                        (ground_truth(j,2) - results(k,2))^2 < d^2 ...

                        && res_tally(k) == 0

                    true_pos = true_pos + 1;

                    false_neg = false_neg - 1;

                    res_tally(k) = 1;

                    break

                end

            end

        end

    end

    precision(h) = true_pos/total_pos;

    recall(h) = true_pos/(true_pos + false_neg);

    f1(h) = 2*precision(h)*recall(h)/(precision(h) + recall(h));

    waitbar(h/size(threshold,2));

end

close(g)

[M,I] = max(f1);

disp(['The F1 score reaches a maximum of ' num2str(M) ' when the threshold is ' num2str(threshold(I))])

figure

plot(precision, recall)

title('precision-recall curve')

xlabel('precision')

ylabel('recall')

xlim([0 1])

ylim([0 1])

figure

plot(threshold, f1)

title('f1 score')

xlabel('threshold')

ylabel('f1 score')

xlim([0 1])

ylim([0 1])

figure

plot(threshold, precision)

title('precision')

xlabel('threshold')

ylabel('precision')

xlim([0 1])

ylim([0 1])

figure

plot(threshold, recall)

title('recall')

xlabel('threshold')

ylabel('recall')

xlim([0 1])

ylim([0 1])