% 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 [testObj, testMask, testRaster] = loadTestData(testModelFilename, testMaskFilename, refInfo, refDTM, offsets, params)


% INPUTS
%   testModelFilename   - filename of point cloud (LAS or LAZ)  
%   testMaskFilename    - filename of geotiff mask of indicating regions to test.  
%   refInfo             - geotiffinfo struct representing image space to compute metrics
%   refDTM              - terrain image. Must be in space defined by refInfo
%   offsets             - registration offsets [East, North, Height] in meters 
%   params              - See below
%
%   NOTE: All data need to be georeferenced in UTM to the same UTM zone
%
%
% OUTPUTS
%   testObj      raster of test objects relative to DTM  
%                (Excluded pixels are set to NaN)
%   testMask     binary mask of test object footprints
%   testRaster   raster image of test model point cloud
%                (No Data pixels are set to NaN)
%

% Parameters:
QUANTIZE        = false;
DSM_FILL_VALUE  = -9999;
DSM_METHOD      = 1;
MASK_TRUE       = 1;
MASK_APPLY_OFFSETS  = true;
MODEL_APPLY_OFFSETS = true;
DO_PLOTS        = false;
SAVE_FIGS       = false;
SAVE_PLOTS      = false;
SAVE_DIR        = '.';
MODEL_NAME      = 'Test Model';
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


if ~exist('offsets', 'var') ||  isempty(offsets)  
    offsets = [0 0 0];
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

% ---------------------------------------------------------------------
fprintf('\tLoading test data...');    btic = tic;
   
% Object mask
testMaskImg  = geotiffread(testMaskFilename);
testMaskInfo = geotiffinfo(testMaskFilename);
    

% Transform object mask to reference space and apply registration offsets
if MASK_APPLY_OFFSETS
    testMaskInfo.RefMatrix(3,1) = testMaskInfo.RefMatrix(3,1) + offsets(1);  
    testMaskInfo.RefMatrix(3,2) = testMaskInfo.RefMatrix(3,2) + offsets(2);
end
testMask  = imwarpGeoTiff(double(testMaskImg == MASK_TRUE), testMaskInfo, refInfo, double(false));


% Load the model  (LAS or GEOTIFF)
[~,~,ext] = fileparts(testModelFilename);
IS_LAS = any(strcmpi({'.las', '.laz'}, ext));
if IS_LAS
    testModel = lasdata(testModelFilename); 
    fprintf('     (%0.3f)\n', toc(btic));

    fprintf('\tProcessing test data... ');    btic = tic;
    
    % Apply registation offsets
    if MODEL_APPLY_OFFSETS
        testModel.x = testModel.x+offsets(1);
        testModel.y = testModel.y+offsets(2);
        testModel.z = testModel.z+offsets(3);
    end

    % Align point cloud with reference image space
    [row, col]   = map2pix(refInfo.RefMatrix, testModel.x, testModel.y);

    % Bin the test cloud to create a DSM
    outSize = [refInfo.Height refInfo.Width];
    
    % If test mask is supplied, pixels with no samples are inpained
    [testRaster, ~, numPtsPerBin] = rasterData([col row], testModel.z, outSize, ~isnan(testMask), DSM_METHOD);

else
    testModelImg   = geotiffread(testModelFilename);
    testModelInfo  = geotiffinfo(testModelFilename);
    fprintf('     (%0.3f)\n', toc(btic));
    
    fprintf('\tProcessing test data... ');    btic = tic;
    
    % Apply registation offsets
    if MODEL_APPLY_OFFSETS
        testModelInfo.RefMatrix(3,1) = testModelInfo.RefMatrix(3,1) + offsets(1);  
        testModelInfo.RefMatrix(3,2) = testModelInfo.RefMatrix(3,2) + offsets(2);
    end
    
    % Project to reference image space
    testRaster     = imwarpGeoTiff(testModelImg, testModelInfo, refInfo, DSM_FILL_VALUE);
    testRaster     = testRaster + offsets(3);
end

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


% Make object only DSM
testMask = ~isnan(testRaster) & (testMask==1);
testObj = testRaster - refDTM;
testObj(testMask == 0) = nan;


fprintf('%0.3f\n', toc(btic));  % END PROCESSING  TEST DATA
      


% ---------- Make some plots ----------------------------------------------
if DO_PLOTS
    fprintf('\tMaking plots of test data...'); btic = tic;

    figure(DO_PLOTS)
    imagesc(testRaster); axis image ij; colorbar; %caxis(caxDTM);
    title(sprintf('%s - Mask (testRaster)', model_prefix_str));

    if SAVE_PLOTS
        str = sprintf('%s_mask', 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

    if IS_LAS
        figure(DO_PLOTS+1)
        imagesc(numPtsPerBin); axis image ij;  colorbar;
        title(sprintf('%s - Samples Per Pixel (N)', model_prefix_str));

        if SAVE_PLOTS
            str = sprintf('%s_counts', 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
    end

    figure(DO_PLOTS+2)
    imagesc(testMask); axis image ij;  colorbar
    title(sprintf('%s - Mask (testMask)', model_prefix_str));

    if SAVE_PLOTS
        str = sprintf('%s_mask', 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+3)
    imagesc(testObj); axis image ij; colorbar; %caxis(caxObj);
    title(sprintf('%s - Objects (testObj)', model_prefix_str));

    if SAVE_PLOTS
        str = sprintf('%s_mask', 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('%0.3f\n', toc(btic));
end



