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

% 2D subsampling of input mesh

% Version History
%   2017/03/28  Marc Bosch       Initial version
function [inputpcout]=decimationgrid(inputpc1,gtAOIpath,DSlevel)
inputpcout=inputpc1;
if DSlevel>1
    %convert2dgrid
    [aoigt,aoiRef]=geotiffread(gtAOIpath);
    aoiinfo=geotiffinfo(gtAOIpath);
    [xy]=round(map2pix(aoiinfo.RefMatrix,inputpc1.x,inputpc1.y));  
    ofy=min(0,min(xy(:,2)));
    ofx=min(0,min(xy(:,1)));
    grid2dx=zeros(max(xy(:,2))-ofy+1,max(xy(:,1))-ofx+1);
    grid2dy=zeros(max(xy(:,2))-ofy+1,max(xy(:,1))-ofx+1);
    grid2dz=zeros(max(xy(:,2))-ofy+1,max(xy(:,1))-ofx+1);
    for ij=1:length(inputpc1.x)
        grid2dz(xy(ij,2)-ofy+1,xy(ij,1)-ofx+1)=double(inputpc1.z(ij));
        grid2dy(xy(ij,2)-ofy+1,xy(ij,1)-ofx+1)=double(inputpc1.y(ij));
        grid2dx(xy(ij,2)-ofy+1,xy(ij,1)-ofx+1)=double(inputpc1.x(ij));
    end
%     grid2dx=medfilt2(grid2dx,[4 4]);
%     grid2dy=medfilt2(grid2dy,[4 4]);
%     grid2dz=medfilt2(grid2dz,[4 4]);
    temp=grid2dx(1:DSlevel:size(grid2dx,1),1:DSlevel:size(grid2dx,2));
    p=find(temp(:)~=0);
    inputpc1.x=temp(p);
    temp=grid2dy(1:DSlevel:size(grid2dx,1),1:DSlevel:size(grid2dx,2));
    inputpc1.y=temp(p);
    temp=grid2dz(1:DSlevel:size(grid2dx,1),1:DSlevel:size(grid2dx,2));
    inputpc1.z=temp(p);
    [inputpcout.x,posorted]=sort(inputpc1.x);    
    inputpcout.y=inputpc1.y(posorted);
    inputpcout.z=inputpc1.z(posorted);
end