function [rec, prec, ap, sc] = rgbdaprecision(bboxtrue, bbox, categorytag)
% compute precision-recall curve and averaged precision

lapthreshold = 0.5;
% get predicted bounding boxes
it = 0;
for i = 1:length(bbox)
    if length(bbox{i})
       for j = 1:length(bbox{i})
           it = it + 1;
           BB(it,1) = bbox{i}(j).left;
           BB(it,2) = bbox{i}(j).top;
           BB(it,3) = bbox{i}(j).right;
           BB(it,4) = bbox{i}(j).bottom;
           confidence(it,1) = bbox{i}(j).score;
           ids(it,1) = i;
           allcategory{it} = bbox{i}(j).category;
           allinstancenum(it,1) = bbox{i}(j).instance;
           allinstance{it} = [allcategory{it} '_' num2str(bbox{i}(j).instance)];
       end
    end
end
categorycell = unique(allcategory);
instancecell = unique(allinstance);

npos = 0;
for i = 1:length(bboxtrue)
    if length(bboxtrue{i})
       for j  = 1:length(bboxtrue{i})
           if categorytag
              for s = 1:length(categorycell)
                  if strcmp(bboxtrue{i}(j).category, categorycell{s})
                     npos = npos + 1;
                  end
              end
           else
              for s = 1:length(instancecell)
                  ind = find(instancecell{s} == '_');
                  category = instancecell{s}(1:ind(end)-1);
                  instancenum = str2num(instancecell{s}(ind(end)+1:end));
                  if strcmp(bboxtrue{i}(j).category, category) & bboxtrue{i}(j).instance == instancenum
                     npos = npos + 1;
                  end
              end
           end
        end
    end
end

% sort detections by decreasing confidence
[sc,si] = sort(-confidence);
sc = -sc;
ids = ids(si);
BB = BB(si,:);
allcategory = allcategory(si);
allinstancenum = allinstancenum(si);

% assign detections to ground truth objects
nd = length(confidence);
tp = zeros(nd,1);
fp = zeros(nd,1);
tag = zeros(length(bboxtrue),30);
for d = 1:nd
    % assign detection to ground truth object if any
    bb = BB(d,:);
    i = ids(d);
    ovmax = -inf;
    if length(bboxtrue{i})
       for j = 1:length(bboxtrue{i})
           if categorytag
              if strcmp(allcategory{d}, bboxtrue{i}(j).category)
                 bbgt(1,1) = bboxtrue{i}(j).left;
                 bbgt(1,2) = bboxtrue{i}(j).top;
                 bbgt(1,3) = bboxtrue{i}(j).right;
                 bbgt(1,4) = bboxtrue{i}(j).bottom;

                 bi = [max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
                 iw = bi(3)-bi(1)+1;
                 ih = bi(4)-bi(2)+1;
                 if iw > 0 & ih > 0                
                 % compute overlap as area of intersection / area of union
                 ua = (bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+(bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-iw*ih;
                 ov = iw*ih/ua;
                    if ov > ovmax
                       ovmax = ov;
                       jmax = j;
                    end
                 end
              end
           else
              if strcmp(allcategory{d}, bboxtrue{i}(j).category) & allinstancenum(d) == bboxtrue{i}(j).instance
                 bbgt(1,1) = bboxtrue{i}(j).left;
                 bbgt(1,2) = bboxtrue{i}(j).top;
                 bbgt(1,3) = bboxtrue{i}(j).right;
                 bbgt(1,4) = bboxtrue{i}(j).bottom;

                 bi = [max(bb(1),bbgt(1)) ; max(bb(2),bbgt(2)) ; min(bb(3),bbgt(3)) ; min(bb(4),bbgt(4))];
                 iw = bi(3)-bi(1)+1;
                 ih = bi(4)-bi(2)+1;
                 if iw > 0 & ih > 0
                 % compute overlap as area of intersection / area of union
                 ua = (bb(3)-bb(1)+1)*(bb(4)-bb(2)+1)+(bbgt(3)-bbgt(1)+1)*(bbgt(4)-bbgt(2)+1)-iw*ih;
                 ov = iw*ih/ua;
                    if ov > ovmax
                       ovmax = ov;
                       jmax = j;
                    end
                 end
              end
           end
       end
    else
       ovmax = 0.0;
       jmax = j;
    end
    % assign detection as true positive/don't care/false positive
    if ovmax >= lapthreshold
       if ~tag(i,jmax)
           tp(d) = 1;            % true positive
           tag(i,jmax) = 1;
       else
           fp(d) = 1;            % false positive (multiple detection)
       end
    else
        fp(d) = 1;                    % false positive
    end
end

% compute precision/recall
fp = cumsum(fp);
tp = cumsum(tp);
rec = tp/npos;
prec = tp./(fp+tp);

ap = vocap(rec,prec);

function ap = vocap(rec, prec)

mrec = [0 ; rec ; 1];
mpre = [0 ; prec ; 0];
for i = numel(mpre)-1:-1:1
    mpre(i) = max(mpre(i),mpre(i+1));
end
i = find(mrec(2:end)~=mrec(1:end-1))+1;
ap = sum((mrec(i)-mrec(i-1)).*mpre(i));


