%This is the code I used so far. For two slices of my images.
im1 = loop_new(:,:,1,1);
im2 = loop_new(:,:,2,1);
T1 = graythresh(im1);
T2 = graythresh(im2);
bw1 = imbinarize(im1,T1/2.8);
bw2 = imbinarize(im2,T2/2.8);
bw11 = medfilt2(bw1,[5 5]);
bw22 = medfilt2(bw2,[5 5]);
h = fspecial('gaussian',[5 5],2);
b11 = imfilter(double(bw11),h);
b22 = imfilter(double(bw22),h);
b11 = b11~=0;
b22 = b22~=0;
b111 = imfill(b11,'holes');
Mask1 = b111;
b222 = imfill(b22,'holes');
Mask2 = b222;
A = bwconncomp(b111);
B = bwconncomp(b222);
A1 = cell2mat(A.PixelIdxList(1));
B1 = cell2mat(B.PixelIdxList(1));
b111(A1) =0;
b222(B1) =0;
Mask(:,:,1) = Mask1 - b111;
Mask(:,:,2) = Mask2 - b222;
Mask = logical(Mask);
mask_brain =double(repmat(Mask(:,:,:),[1 1 1 120])); %128*128*2*120
output_brain = loop_new .* mask_brain; %128*128*2*120