binsurface.m

function [node, elem] = binsurface(img, nface)
%
% [node,elem]=binsurface(img,nface)
%
% fast isosurface extraction from 3D binary images
%
% author: Qianqian Fang, <q.fang at neu.edu>
%
% input:
%   img:  a 3D binary image
%   nface: nface=3 or ignored - for triangular faces,
%          nface=4 - square faces
%          nface=0 - return a boundary mask image via node
%
% output:
%   elem: integer array with dimensions of NE x nface, each row represents
%         a surface mesh face element
%   node: node coordinates, 3 columns for x, y and z respectively
%
% the outputs of this subroutine can be easily plotted using plotmesh()
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

dim = size(img);
if (length(dim) < 3)
    dim(3) = 1;
end
newdim = dim + [1 1 1];

% find the jumps (0->1 or 1->0) for all directions
d1 = diff(img, 1, 1);
d2 = diff(img, 1, 2);
d3 = diff(img, 1, 3);
[ix, iy] = find(d1 == 1 | d1 == -1);
[jx, jy] = find(d2 == 1 | d2 == -1);
[kx, ky] = find(d3 == 1 | d3 == -1);

% compensate the dim. reduction due to diff, and
% wrap the indices in a bigger array (newdim)
ix = ix + 1;
[iy, iz] = ind2sub(dim(2:3), iy);
iy = sub2ind(newdim(2:3), iy, iz);

[jy, jz] = ind2sub([dim(2) - 1, dim(3)], jy);
jy = jy + 1;
jy = sub2ind(newdim(2:3), jy, jz);

[ky, kz] = ind2sub([dim(2), dim(3) - 1], ky);
kz = kz + 1;
ky = sub2ind(newdim(2:3), ky, kz);

id1 = sub2ind(newdim, ix, iy);
id2 = sub2ind(newdim, jx, jy);
id3 = sub2ind(newdim, kx, ky);

if (nargin == 2 && nface == 0)
    elem = [id1; id2; id3];
    node = zeros(newdim);
    node(elem) = 1;
    node = node(2:end - 1, 2:end - 1, 2:end - 1) - 1;
    return
end

% populate all the triangles
xy = newdim(1) * newdim(2);

if (nargin == 1 || nface == 3)  % create triangular elements
    elem = [id1 id1 + newdim(1) id1 + newdim(1) + xy; id1 id1 + newdim(1) + xy id1 + xy];
    elem = [elem; id2 id2 + 1 id2 + 1 + xy; id2 id2 + 1 + xy id2 + xy];
    elem = [elem; id3 id3 + 1 id3 + 1 + newdim(1); id3 id3 + 1 + newdim(1) id3 + newdim(1)];
else                       % create box elements
    elem = [id1 id1 + newdim(1) id1 + newdim(1) + xy id1 + xy];
    elem = [elem; id2 id2 + 1 id2 + 1 + xy id2 + xy];
    elem = [elem; id3 id3 + 1 id3 + 1 + newdim(1) id3 + newdim(1)];
end
% compress the node indices
nodemap = zeros(max(elem(:)), 1);
nodemap(elem(:)) = 1;
id = find(nodemap);

nodemap = 0;
nodemap(id) = 1:length(id);
elem = nodemap(elem);

% create the coordiniates
[xi, yi, zi] = ind2sub(newdim, id);

% assuming the origin [0 0 0] is located at the lower-bottom corner of the image
node = [xi(:), yi(:), zi(:)] - 1;

if (nargin == 1 || nface == 3)
    [node, elem] = meshcheckrepair(node, elem);
end