BG1 = mean(ag,3);
BGstd = std(ag,0,3);

for i=1:size(ag,3)
    FG = BGthres(ag(:,:,i), BG1, BGstd);
    v(i) = size(find(FG),1);
    FGag(:,:,i) = FG;
end

vm = mean(v);
vlow = find(v < vm);
BG2 = mean(ag(:,:,vlow),3);
BG2std = std(ag(:,:,vlow),0,3);
BGdiff = (abs(BG2-BG1));
BGhd = histeq(BGdiff);

for i=1:size(ag,3)
    FG = BGthres(ag(:,:,i), BG2, BG2std);
%    FG2ag(:,:,i) = FG;
    se=strel('square',5);
    FG3ag(:,:,i) = imclose(imclose(imopen(imopen(FG,se),se),se),se);
%    FG3ag(:,:,i) = imclose(imopen(FG,se),se);
    v2(i) = size(find(FG3ag(:,:,i)),1);
    ED(:,:,i) = imdilate(FG3ag(:,:,i),se);
end
FG4 = ED - FG3ag;
FG4 = (FG4 == 1);
