smoothbinvol

 vol=smoothbinvol(vol,layer)

 perform a memory-limited 3D image smoothing

 author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>

 input:
     vol: a 3D volumetric image to be smoothed
     layer: number of iterations for the smoothing

 output:
     vol: the volumetric image after smoothing

 -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)

0001 function vol=smoothbinvol(vol,layer)
0002 %
0003 % vol=smoothbinvol(vol,layer)
0004 %
0005 % perform a memory-limited 3D image smoothing
0006 %
0007 % author: Qianqian Fang <fangq at nmr.mgh.harvard.edu>
0008 %
0009 % input:
0010 %     vol: a 3D volumetric image to be smoothed
0011 %     layer: number of iterations for the smoothing
0012 %
0013 % output:
0014 %     vol: the volumetric image after smoothing
0015 %
0016 % -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
0017 %
0018 
0019 dim=size(vol);
0020 dxy=dim(1)*dim(2);
0021 fulllen=prod(dim);
0022 
0023 weight=1./6.;
0024 step=4000;
0025 
0026 % in case vol is a logical
0027 vol=double(vol);
0028 offs=[1,-1,dim(1), -dim(1),dxy, -dxy];
0029 
0030 for i=1:layer
0031   % find all non-zero values
0032   idx=find(vol);
0033   % get the neighbors of all the non-zero values
0034   % this may cause wrapping -- TODO
0035   val=vol(idx);
0036   for k=1:6
0037     nextidx=idx+offs(k);
0038     % find all 1-valued voxels that are located within the domain
0039       goodidx=find(nextidx>0 & nextidx<fulllen);
0040     % for all neighboring voxels, add a fraction from the non-0 voxels
0041     % problematic when running in parallel (racing)
0042     len=length(goodidx);
0043     % control granualarity with step
0044     if(len>step)
0045         for j=1:step:len-step
0046             vol(nextidx(goodidx(j:j+step-1)))=vol(nextidx(goodidx(j:j+step-1)))+weight*val(goodidx(j:j+step-1));
0047         end
0048         vol(nextidx(goodidx(j+step:end)))=vol(nextidx(goodidx(j+step:end)))+weight*val(goodidx(j+step:end));
0049     else
0050         vol(nextidx(goodidx))=vol(nextidx(goodidx))+weight*val(goodidx);
0051     end
0052     % the above line may change the values of the non-zero voxels, recover
0053     % them
0054   end
0055   vol(idx)=val;
0056 end