% 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.

% Reads the lidar las file, removes non-unique entries and sorts data using
% X coordinates, triangulates to create mesh file, and reads AOI mask
% according 
%
% Version History
%   2017/03/28  Marc Bosch       Initial version

function meshreference=get_referencedata(referencepcpath,gtAOIpath,pointspacing,Surfacemetrics)

referencepc = lasdata(referencepcpath);
%referencepc = lasdataread(referencepc,referencepcpath);
tic;

pcMat=[referencepc.x referencepc.y];
[pcMatu,ppo]=unique(pcMat,'rows','stable');
[referencepc.y,posorted]=sort(pcMatu(:,2));
referencepc.x=pcMatu(posorted,1);
temp=referencepc.z;
referencepc.z=temp(ppo);
temp=referencepc.z;
referencepc.z=temp(posorted);
a=referencepc.x(1:end);
b=[referencepc.x(2:end);min(referencepc.x(1:end)); ];
estimated_pointspacing = median(abs(a-b));
DS_downs=floor(pointspacing/estimated_pointspacing);
if DS_downs>1
    temp=referencepc.x(1:DS_downs:end);
    referencepc.x=temp;
    temp=referencepc.y(1:DS_downs:end);
    referencepc.y=temp;
    temp=referencepc.z(1:DS_downs:end);
    referencepc.z=temp;
end

%read AOI mask and extract region of Lidar gt
if ~isempty(gtAOIpath)==1
    AOImask=alignAOImask(referencepc,gtAOIpath,1);
else
    Classbuilding=6;% class label for buildings in AOI
    AOImask=Classbuilding*ones(length(referencepc.x),1);
end
Ntiles=1;Npointtile=length(referencepc.x);
if length(referencepc.x)>1.5e9
    %tile processing
    Ntiles = floor(length(referencepc.x)/1.5e9);
    Npointtile = floor(length(referencepc.x)/Ntiles);
end

%triangulation for point cloud to mesh conversion
for tid=1:Ntiles
    tic;
    start_point=(tid-1)*Npointtile+1;
    end_point=tid*Npointtile;
    if tid==Ntiles
        end_point = length(referencepc.x);
    end
    if Surfacemetrics==1
        trireference = delaunay(referencepc.x(start_point:end_point),referencepc.y(start_point:end_point));
    else
        trireference=zeros(3);
    end
    %ground truth mesh
    meshtilereference=mesh2struct(trireference', ...
        [referencepc.x(start_point:end_point)';referencepc.y(start_point:end_point)';referencepc.z(start_point:end_point)'], ...
        pointspacing,1);
    meshreference(tid) = meshtilereference;
    %to be done once mesh data is available externally as part of the input solution
    %trireference = populatefromFile(referencepc.x(1:end),referencepc.y(1:end));
    %ground truth mesh
    toc;
end
%rest of data
for tid=1:Ntiles
    start_point=(tid-1)*Npointtile+1;
    end_point=tid*Npointtile;
    if tid==Ntiles
        end_point = length(referencepc.x);
    end
    meshreference(tid).AOImask = AOImask(start_point:end_point)';
    % write_mesh2ply(meshreference.node_xyz',meshreference.triangle_node','C:\Users\boschm1\src\Algorithms\Core3D\metrics\Matlab\reference.ply');
    meshreference(tid).referencepcpath=referencepcpath;
    meshreference(tid).gtAOIpath=gtAOIpath;
    meshreference(tid).pointspacing=pointspacing;
    meshreference(tid).estimated_pointspacing=estimated_pointspacing;
    
end