function mask_thresh = fmri_mask_thresh(fmri_file);

%FMRI_MASK_THRESH
%
% Tries to find a good threshold that defines the brain from fmri data.
%
% MASK_THRESH = FMRI_MASK_THRESH( FMRI_FILE );
%
% MASK_THRESH = the first 'dip' in the histogram of the first frame of
% FMRI_FILE by searching up to the largest local max; 
% uses max(FMRI_FILE)/4 if unsuccessful.

nbin=100;
dbin=10;
dm=fmris_read_image(fmri_file,0,0);
numslices=dm.dim(3)
dm=fmris_read_image(fmri_file,1:numslices,1);
[freq, mask]=hist(dm.data(:),nbin);
fm=freq(1:nbin-2*dbin);
f0=freq(1+dbin:nbin-dbin);
fp=freq(1+2*dbin:nbin);
h=(abs(f0-fm)+abs(f0-fp)).*(f0>fp).*(f0>fm);
if any(h)
    mh=min(find(h==max(h))+dbin);
    mask_thresh=max(mask(find(freq==min(freq(1:mh)))))
else
    mask_thresh=max(mask)/4
end

return