function [ dem_struct ] = dem2struct( data, R, aoi, params )

gsd = max(R.CellExtentInWorldX,R.CellExtentInWorldY);

R.XWorldLimits = R.XWorldLimits - R.CellExtentInWorldX * R.XIntrinsicLimits(1);
R.YWorldLimits = R.YWorldLimits - R.CellExtentInWorldY * R.YIntrinsicLimits(1);
bounds = [R.XWorldLimits;R.YWorldLimits];

if(strcmp(R.ColumnsStartFrom,'north'))
    data = flipud(data);
end

if(~isempty(aoi))
    inds = (1:numel(data))';
    [rs,cs] = ind2sub(size(data),inds);
    xs = (cs-1)*gsd+bounds(1);
    ys = (rs-1)*gsd+bounds(2);
    
    inside = isPointInPolygon([xs,ys],[aoi.X',aoi.Y']);
    data(inds(~inside)) = params.no_data;
end

%Maximum DEM
max_dem = double(data);
max_dem(data == params.no_data) = nan;
max_dem = demFill(max_dem,2);

%Minimum DEM
min_dem = max_dem;

%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)) = min(max_dem(:),[],'omitnan');
min_dem(isnan(min_dem)) = max(min_dem(:),[],'omitnan');

%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

