function [ dem_struct ] = demGenerate( pc, aoi, params )
disp('Create DEM');
gsd = params.gsd;

if(~isempty(aoi))
    inside = isPointInPolygon([pc.x,pc.y],[aoi.X',aoi.Y']);
    pc.x(~inside) = [];
    pc.y(~inside) = [];
    pc.z(~inside) = [];
    pc.header.min_x = min(pc.x);
    pc.header.max_x = max(pc.x);
    pc.header.min_y = min(pc.y);
    pc.header.max_y = max(pc.y);
    pc.header.min_z = min(pc.z);
    pc.header.max_z = max(pc.z);
end
bounds = [pc.header.min_x,pc.header.max_x;pc.header.min_y,pc.header.max_y];

ref_xi = (pc.x-bounds(1))/gsd+1;
ref_yi = (pc.y-bounds(2))/gsd+1;
indices = round([ref_yi,ref_xi]);

%Maximum DEM
max_dem = accumarray(indices,pc.z,[],@max,nan);
max_dem = demFill(max_dem,2);

%Minimum DEM
min_dem = accumarray(indices,pc.z,[],@min,nan);
min_dem = demFill(min_dem,2);

%DEM boundary
boundary = bwboundaries(~isnan(max_dem));
boundary = fliplr(boundary{1});
boundary = simplifyPolygon(boundary,1);
boundary = (boundary - 1)*gsd;
boundary(:,1) = boundary(:,1)+bounds(1);
boundary(:,2) = boundary(:,2)+bounds(2);
[t1,t2] = poly2cw(boundary(:,1),boundary(:,2));
boundary = [t1,t2];

%Remove NaNs
max_dem(isnan(max_dem)) = pc.header.min_z;
min_dem(isnan(min_dem)) = pc.header.max_z;

%Results
dem_struct = [];
dem_struct.gsd = gsd;
dem_struct.bounds = bounds;
dem_struct.boundary = boundary;
dem_struct.max_dem = max_dem;
dem_struct.min_dem = min_dem;
end

