% Copyright 2017 The Johns Hopkins University Applied Physics Laboratory.
% Licensed under the MIT License. See LICENSE.txt in the project root for full license information.
%
%
% Reference:M. Bosch, A. Leichtman, D. Chilcott, H. Goldberg, M. Brown, 
% METRIC EVALUATION PIPELINE FOR 3D MODELING OF URBAN SCENES",
% International Archives of the Photogrammetry, Remote Sensing and 
% Spatial Information Sciences, Volume XLII-1/W1, 2017.

%Systems tested on:
%Windows7 16GB RAM (minimum required) Visual Studio 10, 12, 13, Matlab2015b
%or newer
%Matlab toolboxes needed: image_toolbox, map_toolbox
%Linux RedHat 128GB RAM (registration code only tested on Windows)

warning('off','all')
%clear all
close all
%add dependencies
%Set-up paths for registration and triangular subdivision
addpath(genpath('pipelinePCS'))
addpath(genpath('pipelineVectors'))
addpath(genpath('pipeline25D'))
addpath(genpath('3rdparty'))
addpath(genpath('Registration'))

%Input data
pathinputpc='./data_example/JHUAPL_MasterProvisional2.las';
referencepcpath='./data_example/LIDAR_MasterProvisional2.las';
gtAOIpath='./data_example/SanFernando_AOI2_GT_refined.tif';

%Configuration
pointspacing = 0.3;%point spacing of ground truth model used
visualflag=false;
metrics25flag=1;% turn to 0 to disabled them
metricsvectorsflag=1;% turn to 0 to disabled them
metricspcsflag=1;% turn to 0 to disabled them


%read input data, subdivide and register
fprintf('Reading solution data and creating mesh...\n');
pathinputpc=get_fullpath(pathinputpc);
writelas=1;writeply=0;
meshinput=get_inputdata(pathinputpc,gtAOIpath,pointspacing,writelas,writeply);
%perform registration
fprintf('Registering dense mesh vertices to ground truth Lidar...\n');
node_xyz=meshinput.node_xyz;%input vertices
inputfilename=meshinput.filenameLAS;
referencepcpath=get_fullpath(referencepcpath);
if ~isempty(gtAOIpath)
    gtAOIpath=get_fullpath(gtAOIpath);
end
[filenameLAS_registered,XYZoffset,registeredpc] = XYZregistration(inputfilename,referencepcpath,node_xyz,pointspacing);
meshinput.filenameLAS_registered=filenameLAS_registered;
%replace mesh vertices coordinates with registered ones, faces and edges remain the
%same
meshinput.node_xyz(1,:) = registeredpc.x;
meshinput.node_xyz(2,:) = registeredpc.y;
meshinput.node_xyz(3,:) = registeredpc.z;  

    
    
if metrics25flag~=0
	% %2.5D metrics. It calculates volumetric overlap with ground truth data
	files25D.refDSMFilename      = './data_example/SanFernando_AOI2_DSM.tif';
	files25D.refDTMFilename      = './data_example/SanFernando_AOI2_DTM.tif';
	files25D.refMaskFilename     = './data_example/SanFernando_AOI2_GT_refined.tif';
	files25D.testModelFilename   = meshinput.filenameLAS_registered;
	files25D.testMaskFilename    = './data_example/JHUAPL_MasterProvisional2_mask.tif';
	files25D.testOffsetFilename  = './data_example/JHUAPL_MasterProvisional2_offsets.txt';
	params25D.DSM_FILL_VALUE_REF   = -10000; 
	params25D.DTM_FILL_VALUE_REF   = -10000; 
	params25D.BUILDING_VALUE_REF   = 6; 
	params25D.BUILDING_VALUE_TEST  = 240; 
	params25D.DO_PLOTS        = true;
	params25D.SAVE_PLOTS      = false;
	params25D.SAVE_FIGS       = false;
	params25D.SAVE_DIR        = '.';
	params25D.MODEL_NAME      = 'JHUAPL';
	params25D.REF_IMAGE       = 3;
	params25D.RESULTS_FILE    = 'pipeline25D_results.txt';
	[metrics2D,  metrics3D] = pipeline25D_run(files25D, params25D);
end
%vector metrics. It calculates vertical and horizontal accuracy errors
%(reference free)
if metricsvectorsflag~=0
 aoi_filename = '.\data_example\AOI_MasterProvisional2.kml';
 data_filename = meshinput.filenameLAS_registered;%'.\data_example\LIDAR_MasterProvisional2.las';
 bldg_filename = '.\data_example\Buildings_MasterProvisional2.kml';
 road_filename = '.\data_example\Roads_MasterProvisional2.kml';%[];
 params = get_default_params();
 params.dem.gsd = pointspacing;
 params.vrms.rms_thr = 2; %For lidar point cloud , use 0.5
 params.hrms.display = visualflag;
 params.vrms.display = visualflag;
 [hrms,hrms_uncertainty, vrms,vrms_uncertainty] = vector_metrics_run(data_filename,aoi_filename,bldg_filename,road_filename,params);
end


%PCS metrics. It calculates positional accuracy, surface and curvature
%accuracies (perceptual metrics)
if metricspcsflag~=0
	registration_done=1;
	% Summary report file
	ss=fix(clock);
	timedata=[num2str(ss(1)) num2str(ss(2)) num2str(ss(3)) num2str(ss(4)) num2str(ss(5)) num2str(ss(6))];
	fid=fopen(['PCSmetric_results_' timedata '.txt'],'a');

	[outmesh] = PCSmetric_run(meshinput,referencepcpath,gtAOIpath,pointspacing,registration_done,visualflag,fid);

	fclose(fid);
end