function [control_ind, control_pts] = calc_region_corners(R1, reg_x, reg_y)
control_ind = {};
xr=range(find(any(R1)));
yr=range(find(any(R1')));
dx=ceil(diff(xr)/reg_x);
dy=ceil(diff(yr)/reg_y);
control_pts = [];

    for j=1:reg_y
for i=1:reg_x

        control_pts = [control_pts; [xr(1)+dx*(i-1)+.5, yr(1)+dy*(j-1)+.5; xr(1)+dx*(i-1)+.5, yr(1)+dy*(j-1)+dy+.5; xr(1)+dx*(i-1)+dx+.5, yr(1)+dy*(j-1)+dy+.5; xr(1)+dx*(i-1)+dx+.5, yr(1)+dy*(j-1)+.5;]];
        control_ind{end+1} = [xr(1)+dx*(i-1)+.5, yr(1)+dy*(j-1)+.5; xr(1)+dx*(i-1)+.5, yr(1)+dy*(j-1)+dy+.5; xr(1)+dx*(i-1)+dx+.5, yr(1)+dy*(j-1)+dy+.5; xr(1)+dx*(i-1)+dx+.5, yr(1)+dy*(j-1)+.5;];
    end
end
control_pts = unique(control_pts, 'rows');
for i=1:length(control_ind)
    c = control_ind{i};
    ind = [];
    for ii=1:size(c,1)
      p = c(ii,:);
      d = sqrt(sum((control_pts - ones(size(control_pts,1),1)*p).^2,2));
      [val, iii] = min(d);
      ind = [ind; iii];
   end
   control_ind{i} = ind;
end

return;