maskdist.m

function dist = maskdist(vol)
%
% dist=maskdist(vol)
%
% return the distance in each voxel towards the nearest label boundaries
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    vol: a 2D or 3D array
%
% output:
%    dist: an integer array, storing the distance, in voxel unit, towards
%          the nearest boundary between two distinct non-zero voxels, the
%          zero voxels in the domain and space outside of the array
%          are also treated as a unique non-zero label. If the goal is to
%          get the minimum distance measured from the center of the voxel,
%          one should use (dist-0.5).
%
% example:
%
%    a=ones(60,60,60);
%    a(:,:,1:10)=2;
%    a(:,:,11:20)=3;
%    im=maskdist(a);
%    imagesc(squeeze(im(:,30,:)))
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if (isempty(vol))
    error('input vol can not be empty');
end

vals = unique(vol(:));
if (length(vals) > 256)
    error('it appears that your input is a gray-scale image, you must convert it to binary or labels first');
end

newvol = ones(size(vol) + 2) * max(vals) + 1;
if (ndims(vol) == 2)
    newvol(2:end - 1, 2:end - 1) = vol;
elseif (ndims(vol) == 3)
    newvol(2:end - 1, 2:end - 1, 2:end - 1) = vol;
end

vals(end + 1) = newvol(1, 1, 1);
vals(vals == 0) = [];
newvol(newvol == 0) = newvol(1, 1, 1);

dist = ones(size(newvol)) * inf;

for i = 1:length(vals(:))
    vv = (newvol == vals(i));
    vdist = bwdist(vv);
    vdist(vdist == 0) = inf;
    dist = min(dist, vdist);
end

if (ndims(vol) == 2)
    dist = dist(2:end - 1, 2:end - 1);
elseif (ndims(vol) == 3)
    dist = dist(2:end - 1, 2:end - 1, 2:end - 1);
end