%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%  TP_FN_3DTree Function for getting an evaluation of the quality
% Input:  (1) name of the file generated by your method (see bellow for format)
%         (2) (BOOL default: false) the coordinate are from the tree top (used for the GPU result)
%         (3) (INT default: 0.5) the tolerance epsilon [0 1] of overlaps to detect a true positive
%         (4) (BOOL default: true) whether or not we save the result. 
% Output: (1) TruePositive FalsePositive FalseNegative average_overlap_(%) precision recall
%
%the data is formated as follow:
%x y z height_of_the_tree radius_of_the_tree with p(x,y,z) the tree bottom.
%
%zone1: epsilon 0.5
% 25.01.13 Yannick Verdie
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


function[] = TP_FN_3DTree(name_image,varargin)

%-----------------------------------parameters
narginchk(1, 4);

numvarargs = length(varargin);
if numvarargs > 4
    error('get_score:TooManyInputs','requires at most 3 optional inputs');
end

% set defaults for optional inputs
optargs = {false 0.5 true};
optargs(1:numvarargs) = varargin;
[z_up epsilon save_txt] = optargs{:};
epsilon = max(min(epsilon,1),0);
%-----------------------------------

cylinder_GT = dlmread('GT/GT_cylinder.txt');

[pathstr, name_file, ext] = fileparts(name_image);
cylinder = dlmread(name_image);

fprintf('nb Object GT:%d\n',length(cylinder_GT));
fprintf('nb Object Method:%d\n',length(cylinder));


%the mark is on top of the tree, so need to put it on bottom (need for GPU result)
if (z_up)
    cylinder = [cylinder(:,1) cylinder(:,2) cylinder(:,6) cylinder(:,3)-cylinder(:,6) cylinder(:,5)];
end

is_free = ones(length(cylinder_GT),1);
 TP = 0;
 FP = 0;
 
 
 sumAreaInter = 0;
 
for i=1:length(cylinder);%for each object detected
   r = cylinder(i,5);
   h = cylinder(i,4);
   z = cylinder(i,3);
   y = cylinder(i,2);
   x = cylinder(i,1);
   
   GT_r = cylinder_GT(:,5);
   GT_h = cylinder_GT(:,4);
   GT_z = cylinder_GT(:,3);
   GT_y = cylinder_GT(:,2);
   GT_x = cylinder_GT(:,1);
   ind = find((x-GT_x).^2+(y-GT_y).^2<=r+GT_r & is_free(:) == 1);
   if (~isempty(ind))
     %compute intersection for each and take the best
     AreaA = pi*r*r*h;
     bestAinter = 0;
     selectedId = 0;

     for j=1:length(ind)
          R = GT_r(ind(j));
         AreaB = pi*R*R*GT_h(ind(j));
         %area intersection
         d = sqrt( (x-GT_x(ind(j)))^2+(y-GT_y(ind(j))^2));
         areaInter = 0;
         if (d > abs(R-r)) %one is not included inside the other
             if d < r+R    %but they intersect
                 areaInter = r*r*acos((d*d+r*r-R*R)/(2*d*r))+R*R*acos((d*d+R*R-r*r)/(2*d*R))-0.5*sqrt((-d+r+R)*(d+r-R)*(d-r+R)*(d+r+R));
                 areaInter = areaInter*(min(z+h,GT_z(ind(j))+ GT_h(ind(j)))-max(z,GT_z(ind(j))));
            end 
         else              %one is included inside the other
            areaInter = pi*min(R*R,r*r)*(min(z+h,GT_z(ind(j))+ GT_h(ind(j)))-max(z,GT_z(ind(j))));
         end
         areaInter = areaInter/min(AreaA,AreaB);%between 0 and 1
         
          
         if ( areaInter > bestAinter)
             bestAinter = areaInter;
             selectedId = ind(j);
         end
     end
     
     if ( bestAinter > epsilon)

         is_free(selectedId,1) = 0;
          sumAreaInter = sumAreaInter + bestAinter; 
         TP = TP + 1;
     else
         FP = FP + 1; 
     end
     
  else
       FP = FP + 1;
  end
end

FN = length(find(is_free == 1));

precision = TP/(TP+FP);
recall = TP/(TP+FN);
 
data = [TP FP FN sumAreaInter/TP precision recall]

if (save_txt)
dlmwrite([name_file '_data.txt'],data,'delimiter',' ','precision','%10.6g');
end