% Copyright (c) 2015 The Johns Hopkins University/
% Applied Physics Laboratory. (JHU/APL)
% 
% This material was developed by JHU/APL under U.S. Government contract.
% The U.S. Government has Unlimited Rights in this material per
% DFARS 252.227-7013/7014, as appropriate.
% 
% For any other permission, please contact the Legal Office at JHU/APL.
%
function [metrics2D,  metrics3D,  metrics3Ds] = compute25DMetrics( testMask, testObj, refMask, refObj, refInfo, params)

% INPUTS
%   testMask   - binary mask of test objects 
%   testObj    - height image of test objects
%   refMask    - binary mask of reference objects 
%   refObj     - height image of reference objects
%   refInfo    - geotiffinfo struct of evalation space (used for GSD)
%   params     - See below
%
%
% OUTPUTS
%   metrics2D    - 2D metrics
%   metrics3D    - 3D metrics
%   metrics3Ds   - 3D metrics evaluated over 2D TP regions
%

DO_PLOTS        = false;
SAVE_FIGS       = false;
SAVE_PLOTS      = false;
SAVE_DIR        = '.';
MODEL_NAME      = 'Test Model';
RESULTS_FILE    = params.RESULTS_FILE;
if exist('params', 'var') &&  ~isempty(params)
    fn = fieldnames(params);
    for p = 1:numel(fn)
        if exist(fn{p}, 'var')  
            eval(sprintf('%s = params.(fn{p});',fn{p}) );
        end
    end
end


% Name used in plot titles
model_prefix_str = MODEL_NAME;

% String used in save filenames 
badChars = {'/', '\', ' '};
model_prefix = MODEL_NAME;
for p = 1:numel(badChars)
    model_prefix = strrep(model_prefix, badChars{p}, '');
end


%%
refOnlyMask  = double(refMask  .* abs(1-testMask));
testOnlyMask = double(testMask .* abs(1-refMask));
overlapMask  = double(refMask  .*  testMask);


unitArea = abs(refInfo.RefMatrix(1,2)*refInfo.RefMatrix(2,1));


% --- Hard Error ------------------------------------------------------
% -------- False Positive ---------------------------------------------
unitCountFP = nansum(nansum(testOnlyMask,1),2);
oobFP = nansum(nansum(testOnlyMask.*testObj,1),2).*unitArea;

% -------- False Negative ---------------------------------------------
unitCountFN = nansum(nansum(refOnlyMask,1),2);
oobFN = nansum(nansum(refOnlyMask.*refObj,1),2).*unitArea;


% --- Soft Error ------------------------------------------------------
overlap = overlapMask .* (testObj - refObj);
unitCountTP = sum(sum(~isnan(overlap),1),2);

% -------- False Positive -------------------------------------------------
FP =  nansum(nansum( (overlap > 0) .* overlap,1),2).*unitArea;

% -------- False Negative -------------------------------------------------
FN = -nansum(nansum( (overlap < 0) .* overlap,1),2).*unitArea;

% -------- True Positve ---------------------------------------------------
TP =  nansum(nansum( refObj .* overlapMask,1),2).*unitArea - FN;

tolFP = FP + oobFP;
tolFN = FN + oobFN;
tolTP = TP;

if false
    
    % Test model
    figure;
    imagesc(testObj); axis image off; colorbar; caxis([-2 22] );
    title('Performer Model')
    cm = parula(256); colormap( [0.4*[1 1 1]; cm])
    
    % Ref model
    imagesc(refObj); axis image off; colorbar; caxis([-2 22] );
    title('Truth Model')
    cm = parula(256); colormap( [0.4*[1 1 1]; cm])
    
    % FN height;
    tmp1 = abs(refOnlyMask.*refObj); tmp1(isnan(tmp1)) = 0;
    tmp2 = abs(overlap.*(overlap < 0)); tmp2(isnan(tmp2)) = 0;
    imagesc(tmp1+tmp2); axis image off; colorbar; caxis([0 5]); title('False Negative')
    cm = parula(256); colormap( [0.4*[1 1 1]; cm])

    % FP height;
    tmp1 = abs(testOnlyMask.*testObj); tmp1(isnan(tmp1)) = 0;
    tmp2 = abs(overlap.*(overlap > 0)); tmp2(isnan(tmp2)) = 0;
    imagesc(tmp1+tmp2); axis image off; colorbar; caxis([0 10]); title('False Positive')
    cm = parula(256); colormap( [0.4*[1 1 1]; cm])

    % TP height;
    tmp1 = abs(refObj.*(overlap >= 0)); tmp1(isnan(tmp1)) = 0;
    tmp2 = abs(testObj.*(overlap < 0)); tmp2(isnan(tmp2)) = 0;
    imagesc(tmp1+tmp2); axis image off; colorbar; caxis([0 20]); title('True Positive')
    cm = parula(256); colormap( [0.4*[1 1 1]; cm])
    
    %%
    
    % Overlap
    imagesc(overlap)
    caxis([-1 1]*2); colorbar; 
    cm = [parula(265); flipud(parula(265))]; colormap( [0.4*[1 1 1]; cm])
    title({'Correctly Labeled Buildings' 'Height Error'})
    axis image off
    
    imagesc(overlap .* (overlap < 0))
    caxis([-1 0]*2); colorbar; 
    cm = [flipud(parula(265))]; colormap([ 0.4*[1 1 1]; cm ; 0.4*[1 1 1]])
    title({'Correctly Labeled Buildings' 'Under Estimated Height'})
    axis image off
    
    imagesc(overlap .* (overlap > 0))
    caxis([0 1]*2); colorbar; 
    cm = [parula(265)]; colormap([ 0.4*[1 1 1]; cm ; 0.4*[1 1 1]])
    title({'Correctly Labeled Buildings' 'Over Estimated Height'})
    axis image off
    
    
end

metrics2D     = calcMops(unitCountTP, unitCountFN, unitCountFP);
metrics3D     = calcMops(tolTP, tolFN, tolFP);
metrics3Ds    = calcMops(TP, FN, FP);

if ~isempty(RESULTS_FILE)
    fid = fopen(RESULTS_FILE,'w');
    if fid == -1
        warning('Unable to open results files: %s\n', RESULTS_FILE)
    end
    fprintf(fid,'metric,2D,3D\n');
    fn = fieldnames(metrics2D);
    for p = 1:numel(fn)
        fprintf(fid,'%s, %0.5f, %.5f\n', fn{p}, metrics2D.(fn{p}), metrics3D.(fn{p}));   
    end
    fclose(fid)
end
  

typeMap = zeros(size(refObj));
typeMap(refOnlyMask == 1)  = 1;   % False Negetive
typeMap(testOnlyMask == 1) = 2;   % False Positive
typeMap(overlapMask == 1)  = 3;   % True Positive


fprintf('\n====================================\n');
fprintf('2D Metrics\n');

fprintf('\tFalse Positive =\t% 9.1f (m^2)\n', unitCountFP.*unitArea);
fprintf('\tFalse Negetive =\t% 9.1f (m^2)\n', unitCountFN.*unitArea);
fprintf('\tTrue Positive  =\t% 9.1f (m^2)\n', unitCountTP.*unitArea);
fprintf('\tCompleteness =\t%0.5f\n', metrics2D.completeness);
fprintf('\tCorrectness =\t%0.5f\n',  metrics2D.correctness);

fprintf('\n3D Metrics\n');
fprintf('\tFalse Positive =\t% 9.1f (m^3)\n', tolFP);
fprintf('\tFalse Negetive =\t% 9.1f (m^3)\n', tolFN);
fprintf('\tTrue Positive  =\t% 9.1f (m^3)\n', tolTP);

fprintf('\tCompleteness =\t%0.5f\n', metrics3D.completeness);
fprintf('\tCorrectness =\t%0.5f\n',  metrics3D.correctness);


fprintf('\n====================================\n');


% ---------- Make some plots ----------------------------------------------
if DO_PLOTS

    fprintf('\tMaking plots of test data...'); btic = tic;

    figure(DO_PLOTS); clf; cla; ax = [];
    ax(1) = subplot(1,3,1); imagesc(refOnlyMask); axis image ij; title('False Negative'); axis off
    ax(2) = subplot(1,3,2); imagesc(testOnlyMask); axis image ij; title('False Positive'); axis off
    ax(3) = subplot(1,3,3); imagesc(overlapMask); axis image ij; title('True Positive'); axis off
    linkaxes(ax);
    colormap summer
    if SAVE_PLOTS
        str = sprintf('%s_result', model_prefix);
        saveas(gcf, fullfile(SAVE_DIR, sprintf('%s.%s', str, 'png')))
        if SAVE_FIGS
            saveas(gcf, fullfile(SAVE_DIR, sprintf('%s.%s', str, 'fig')))
        end
    end

    figure(DO_PLOTS+1); clf; cla
    imagesc(typeMap); caxis([0 3]); colorbar; colormap([ [0.1 0.3 0.4].*1.2; 1 0.3 0.1; 1 1 0.3; 0.2 0.6 0.3])
    colorbar('Ticks',linspace(0.4,2.6,4),...
        'TickLabels',{'TN','FN','FP','TP'})
    title('Error - Categories')
    axis image ij off
    ax207 = gca;


    if SAVE_PLOTS
        str = sprintf('%s_result_label', model_prefix);
        saveas(gcf, fullfile(SAVE_DIR, sprintf('%s.%s', str, 'png')))
        if SAVE_FIGS
            saveas(gcf, fullfile(SAVE_DIR, sprintf('%s.%s', str, 'fig')))
        end
    end

    errorMap = nan(size(refObj));
    errorMap(testOnlyMask == 1) = testObj(testOnlyMask == 1);  % False Positive
    errorMap(refOnlyMask == 1)  = -refObj(refOnlyMask == 1);   % False Negetive
    errorMap(overlapMask == 1)  = overlap(overlapMask == 1);

    figure(DO_PLOTS+2); cla
    imagesc(errorMap); colorbar
    cm = colormap;
    axis image ij
    title('Error - Height (meters)')
    ax208 = gca;


    if SAVE_PLOTS
        str = sprintf('%s_result_error', model_prefix);
        saveas(gcf, fullfile(SAVE_DIR, sprintf('%s.%s', str, 'png')))
        if SAVE_FIGS
            saveas(gcf, fullfile(SAVE_DIR, sprintf('%s.%s', str, 'fig')))
        end
    end



    upperLim =  5;
    lowerLim = -5;
    errorMap(errorMap > upperLim) = upperLim;
    errorMap(errorMap < lowerLim) = lowerLim;

    cm = jet(128);
    errorMap( isnan(errorMap) ) = lowerLim - (upperLim - lowerLim) / (size(colormap,1) - 1 );

    figure(DO_PLOTS+3)
    imagesc(errorMap); colormap(cm); colorbar

    axis image ij
    title('Error - Height (meters)')
    ax208 = gca;

    cm(1,:) = [1 1 1]*0.2;
    colormap(cm);

    if SAVE_PLOTS
        str = sprintf('%s_result_error_clipped', model_prefix);
        saveas(gcf, fullfile(SAVE_DIR, sprintf('%s.%s', str, 'png')))
        if SAVE_FIGS
            saveas(gcf, fullfile(SAVE_DIR, sprintf('%s.%s', str, 'fig')))
        end
    end
    fprintf('Done (%0.3f)\n', toc(btic));
end




