function [VM]= voidconst(I)
%% keep sphere shape as much as possible
%% smooth
%level=graythresh(I);
close all;
level=0.3;
bw=im2bw(I,level);
%figure; imshow(bw); title('bw, im2bw');
% smoothed (median)
sm_size = 2;       % size of the median filter
S = medfilt2(bw, [sm_size sm_size]);
%figure; imshow(S);title('S, median smoothing');


%% dilation and erosion
bh_radius = 1;
SE = strel('disk', bh_radius);
S1 = imerode(1-S, SE);
%figure; imshow(S1);title('SID, dilation');

bh_radius = 2;
SE = strel('disk', bh_radius);
S2 = imclose(S1,SE);
%figure; imshow(S2);title('SIDE, erotion after dilation');
S3 = imopen(S2,SE);
%figure; imshow(S3);title('SIDE, erotion after dilation');

bh_radius = 5;
SE = strel('disk', bh_radius);
S4 = imclose(S3,SE);
%figure; imshow(S4);title('SIDE, erotion after dilation');

VF=S4;
%% region prop & void matrix
cc = bwconncomp(VF, 8);
graindata = regionprops(cc, 'basic');
grain_areas = [graindata.Area];

%voidmin=100;
voidmin=100;
VM=false(size(I));
for n=1:1:cc.NumObjects;
    grain = false(size(I));
    if grain_areas(n) > voidmin;
        grain(cc.PixelIdxList{n}) = true;
        VM(grain(:,:)==1)=1;
    end
end
%figure, imshow(VM); title('VM, void map');
