% 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 [refObj, refMask, refDTM, refInfo] = loadReferenceData(refDSMFilename, refDTMFilename, refMaskFilename, params)


% INPUTS
%   refDSMFilename      - filename of reference surface model geotiff 
%   refDTMFilename      - filename of reference terrain model geotiff 
%   refMaskFilename     - filename of reference object mask geotiff 
%   params              - See below
%
%   NOTE: All geotiffs need to be georeferenced in UTM to the same UTM zone
% 
% OUTPUTS
%   refObj      refenence object height image (non-object pixels are set to NaN)
%   refMask     mask of reference objects
%   refDTM      terrain model image
%   refInfo     geotiffinfo struct of evalaution image space

% Parameters:
QUANTIZE        = false;
REF_IMAGE       = 3; % Image space to perform evaluation (1:DSM, 2:DTM, 3:MASK)
DSM_FILL_VALUE  = -9999;
DTM_FILL_VALUE  = -9999;
MASK_TRUE       = 1;
DO_PLOTS        = 0;
SAVE_FIGS       = false;
SAVE_PLOTS      = false;
SAVE_DIR        = '.';
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
fprintf('\n')

% ----- Load Reference Data -----------------------------------------------
fprintf('\tLoading reference data...'); btic = tic;

refDSMImg   = geotiffread(refDSMFilename);
refDSMInfo  = geotiffinfo(refDSMFilename);

refDTMImg   = geotiffread(refDTMFilename);
refDTMInfo  = geotiffinfo(refDTMFilename);

refMaskImg  = geotiffread(refMaskFilename);
refMaskInfo = geotiffinfo(refMaskFilename);

fprintf('     (%0.3f)\n', toc(btic));



% ---------- Align Images to DSM ------------------------------------------
fprintf('\tAligning reference data...'); btic = tic;

% Align all reference data to selected image space
switch REF_IMAGE
    case 1
        refInfo = refDSMInfo;
    case 2
        refInfo = refDTMInfo;
    case 3
        refInfo = refMaskInfo;
end

refDTM    = imwarpGeoTiff(refDTMImg, refDTMInfo, refInfo, DTM_FILL_VALUE);
refDSM    = imwarpGeoTiff(refDSMImg, refDSMInfo, refInfo, DSM_FILL_VALUE);

refMask  = imwarpGeoTiff(double(refMaskImg==MASK_TRUE) , refMaskInfo, refInfo, double(false));
refMask = double(~isnan(refMask));

refInvalid  = isnan(refDTM) | isnan(refDSM) | isnan(refMask);

refDTM(refInvalid)  = nan;
refDSM(refInvalid)  = nan;
refMask(refInvalid) = nan;

fprintf('     (%0.3f)\n', toc(btic));

% ---------- Process Reference Data ---------------------------------------
fprintf('\tProcessing reference data...'); btic = tic;


if QUANTIZE
   res = mean([abs(refInfo.RefMatrix(1,2))  abs(refInfo.RefMatrix(2,1))]);
   refDSM   = round(refDSM./res).*res;
   refDTM   = round(refDTM./res).*res;
end

% Calculate relative heights
refHgt   = refDSM - refDTM;

% Create height map of only objective pixels
tmp = double(refMask);
tmp(tmp == 0) = nan;
refObj = refHgt.*tmp;
clear tmp

fprintf('     (%0.3f)\n', toc(btic));



% ---------- Make some plots ----------------------------------------------
if DO_PLOTS
    fprintf('\tMaking plots of reference data...'); btic = tic;
    
    % refDSM
    figure(DO_PLOTS);
    imagesc(refDSM); colorbar; axis image ij
    title('Reference - DSM (refDSM)'); %caxis(caxDTM)
    
    if SAVE_PLOTS
        str = 'ref_dsm';
        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
    
    
    % refDTM
    figure(DO_PLOTS+1);
    imagesc(refDTM); colorbar; axis image ij
    title('Reference - DTM (refDTM)'); %caxis(caxDTM)
    
    if SAVE_PLOTS
        str = 'ref_dtm';
        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
    
    
    % refHgt
    figure(DO_PLOTS+2);
    imagesc(refHgt); colorbar; axis image ij
    title('Reference - Height (refHgt)'); %caxis(caxObj)
    
    if SAVE_PLOTS
        str = 'ref_hgt';
        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
    
    
    % refMask
    figure(DO_PLOTS+3)
    imagesc(refMask); axis image ij; colorbar
    title('Reference - Mask (refMask)');
    
    if SAVE_PLOTS
        str = 'ref_mask';
        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+4)
    imagesc(refObj); axis image ij
    title('Reference - Objects (refObj)');
    colorbar; %caxis(caxObj);
    
    if SAVE_PLOTS
        str = 'ref_obj';
        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('     (%0.3f)\n', toc(btic));
end

end

