function [R, vertex, edge, edge_petal] = autoRemoveSmallRegions(R, vertex, edge, edge_petal, region_size_threshold)
if isempty(R)
    return;
end
for i=1:size(R,3)
    a(i) = sum(sum(R(:,:,i)));
end
thr = region_size_threshold*median(a(find(a>0)));
ind = find((a<thr) &(a>0));
%ind = find(a<thr);

if isempty(ind)
    return;
end
wbar_handle = waitbar(0, 'Removing small regions, please wait...');
for i=ind
    r = R(:,:,i);
    rd = imdilate(r, ones(5));
    [x, y] = find(rd == 1);
    rd = rd(min(x):max(x),min(y):max(y));

    area = 0;
    m = 0;
    ind2 = setdiff(2:size(R,3), ind);
    for ii=ind2
        r2 = R(min(x):max(x),min(y):max(y),ii);
        aa = sum(sum(r2&rd));
        if aa>area
            area = aa;
            m = ii;
        end
    end
    R(:,:,m) = R(:,:,m) | r;
    edge{m} = [edge{m}; edge{i}];
    waitbar(i/length(ind), wbar_handle);
end
close(wbar_handle);
R(:,:,ind) = [];
edge(ind) = [];
a = [];
for i=1:size(R,3)
    a(i) = sum(sum(R(:,:,i)));
end
thr = region_size_threshold*median(a(find(a>0)));
ind = find(a==0);
R(:,:,ind) = [];
return;