This is well done!
I have a minor suggestion for the visualization of the result (nargout == 0). For large No. of slices in a 3D data set, the output of montage() is very hard to grasp. My suggestion is to present only image slices that are part of the segmented volume:
Below line 89, add
% reduce montage size by selecting the interesting slices, only
slices = squeeze( sum(sum( dMask(:,:,1,:) ,1),2) > 0 );
And change line 90 to:
actually there is nothing really new in this algorithm. I only decided to trade a little bit of accuracy for a lot of speed. This is done by using the differen queues which can be regarded as intensity difference bins. I saw a similar approach in a paper by Johan Berglund (2010 in "Magn Reson Imaging") in a slightly differen context, however, I do not know whether he came up with it himself.