function [R2, grid] = manual_click_region_pts(R, grid, axis_handle)

if nargin==2
    figure(1); clf;
    axis_handle = gca;
end
cla(axis_handle);
imagesc(R, 'Parent', axis_handle); colormap copper;
hold(axis_handle, 'on');
pts = [];
% for i=2:length(grid.REGS)
%     region = grid.REGS{i};
%     plot(axis_handle, grid.control_pts(region.vertex,1), grid.control_pts(region.vertex,2), 'd', 'MarkerEdgeColor', 'g',  'MarkerFaceColor', 'g','HitTest', 'off');
%     plot(axis_handle, grid.shape(region.edge,1), grid.shape(region.edge,2), 'o','MarkerEdgeColor', 'w',  'MarkerFaceColor', 'w','HitTest', 'off');
%     pts = [pts; grid.control_pts(region.vertex,:); grid.shape(region.edge,:)];
% end
plot(axis_handle, grid.control_pts(:,1), grid.control_pts(:,2), 'd', 'MarkerEdgeColor', 'g',  'MarkerFaceColor', 'g','HitTest', 'off');
plot(axis_handle, grid.shape(:,1), grid.shape(:,2), 'o','MarkerEdgeColor', 'w',  'MarkerFaceColor', 'w','HitTest', 'off');
pts = [grid.shape; grid.control_pts];
pts = double(pts);
NEW_REGS ={};
R2(:,:,1) = R;
s = double(grid.shape);
for i=2:length(grid.REGS)
    region = grid.REGS{i};
    p1 = plot(axis_handle, grid.control_pts(region.vertex,1), grid.control_pts(region.vertex,2), 'd',  'MarkerEdgeColor', 'b', 'MarkerFaceColor', 'b','HitTest', 'off');
    button = 1;
    vh = [];
    regions_pts = [];
    fprintf('Click on the points of region %d\n', i-1);
    [x, y, button] = ginput(1);
    while button == 1
        d = sqrt(sum((pts - ones(size(pts,1),1)*[x y]).^2, 2));
        [val, ind] = min(d);
        vh(end+1) = plot(axis_handle, pts(ind,1), pts(ind,2), 'o', 'MarkerEdgeColor', 'r', 'MarkerFaceColor', 'r','HitTest', 'off'); 
        regions_pts = [regions_pts;  pts(ind,1), pts(ind,2)];

        [x, y, button] = ginput(1);
    end
    r2 = roipoly(zeros(size(R)), regions_pts(:,1), regions_pts(:,2));
    R2(:,:,end+1) = r2;
    mask = imdilate(r2, strel('disk', 10));
    [y, x] = find(bwperim(mask));
    [in, on] = find(inpolygon(round(s(:,1)), round(s(:,2)), x, y));
    grid.REGS{i}.edge = in;
    [in, on] = find(inpolygon(round(grid.control_pts(:,1)), round(grid.control_pts(:,2)), x, y));
    grid.REGS{i}.vertex = in;
end
return;

%     ng.edge = [];
%     eh = [];
%     fprintf('Now click on the edge points of region 1');
%     [x, y, button] = ginput(1);
%     while button == 1
%         d = sqrt(sum((double(grid.shape) - ones(size(grid.shape,1),1)*[x y]).^2, 2));
%         [val, ind] = min(d);
%         ng.edge = [ng.edge; ind];
%         vh(end+1) = plot(axis_handle, grid.shape(ind,1), grid.shape(ind,2), '*', 'MarkerEdgeColor', 'g', 'MarkerFaceColor', 'g','HitTest', 'off'); 
%         [x, y, button] = ginput(1);
%     end
%     NEW_REGS{i} = ng;
%     delete(eh);
%     delete(vh);  
% end
% grid.REGS = NEW_REGS;
% r2 = zeros(size(R));
% R2(:,:,1) = R;
% opts = [];
% for i=2:length(grid.REGS)
%     wg = grid.REGS{i};
%     pts = [grid.shape(wg.edge,1), grid.shape(wg.edge,2)];
%     pts = [pts; grid.control_pts(wg.vertex,1), grid.control_pts(wg.vertex,2)];
%     rpts = order_xy_pts(pts);
%     r2 = roipoly(zeros(size(R)), rpts(:,1), rpts(:,2));
%     R2(:,:,end+1) = r2;
% end
return;

