function stage_grid = validateGridForStage(stage_grid, xydata,min_number_sectors)

% Extract the triangulation information
T = stage_grid.tri_info.T;
pts = stage_grid.tri_info.pts;
regions_ind = stage_grid.region_tri_ind;



[density, def_regions, non_empty_reg] = calcDeficientRegions(T, pts, regions_ind, xydata, min_number_sectors);

% for each non-empty region
while ~isempty( def_regions)
    i = def_regions(1);

    neigh_tri = findTouchingTriangles(T,  regions_ind{i});
    neighbours = unique(findRegionContainingTri(regions_ind, neigh_tri));
    neighbours(neighbours == i) = [];

    % estimate the density of sectors in those regions
    density_of_neighbours = density(neighbours);

    % only keep the neigbours who are empty
    empty_neighbours = neighbours(density_of_neighbours<min_number_sectors);

    if isempty(empty_neighbours)
        [val, ind] = sort(density_of_neighbours, 'ascend');
        empty_neighbours =neighbours(ind(1));
    end

    % merge those empty neighbours with the non-empty region
    regions_ind{i} = cat(2, regions_ind{i}, [regions_ind{empty_neighbours}]);

    % update the density of the triangles
    regions_ind(empty_neighbours) = [];
    
    % for each region that is not empty, merge surrounding regions if they are empty.
    [density, def_regions, non_empty_reg] = calcDeficientRegions(T, pts, regions_ind, xydata, min_number_sectors);

end
stage_grid.region_tri_ind = regions_ind;
stage_grid.density = density;
return