surf2vol.m

function [img, v2smap] = surf2vol(node, face, xi, yi, zi, varargin)
%
% [img, v2smap]=surf2vol(node,face,xi,yi,zi,'options',values,...)
%
% convert a triangular surface to a shell of voxels in a 3D image
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    node: node list of the triangular surface, 3 columns for x/y/z
%    face: triangle node indices, each row is a triangle
%              if face contains the 4th column, it indicates the label of
%              the face triangles (each face componment must be closed); if
%              face contains 5 columns, it stores a tetrahedral mesh with
%              labels, where the first 4 columns are the element list and
%              the last column is the element label;
%    xi,yi,zi: x/y/z grid for the resulting volume
%        options: 'fill', if set to 1, the enclosed voxels are labeled by 1
%                 'label', if set to 1, the enclosed voxels are labeled by
%                          the corresponding label of the face or element;
%                          setting 'label' to 1 also implies 'fill'.
%
% output:
%    img: a volumetric binary image at position of ndgrid(xi,yi,zi)
%        v2smap (optional): a 4x4 matrix denoting the Affine transformation to map
%             the voxel coordinates back to the mesh space.
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fprintf(1, 'converting a closed surface to a volumetric binary image ...\n');

opt = varargin2struct(varargin{:});
label = jsonopt('label', 0, opt);

elabel = 1;
if (size(face, 2) >= 4)
    elabel = unique(face(:, end));
    if (size(face, 2) == 5)
        label = 1;
        el = face;
        face = [];
        for i = 1:length(elabel)
            fc = volface(el(el(:, 5) == elabel(i), 1:4));
            fc(:, 4) = elabel(i);
            face = [face; fc];
        end
    end
else
    fc = face;
end

img = zeros(length(xi), length(yi), length(zi), class(elabel));

for i = 1:length(elabel)
    if (size(face, 2) == 4)
        fc = face(face(:, 4) == elabel(i), 1:3);
    end
    im = surf2volz(node(:, 1:3), fc(:, 1:3), xi, yi, zi);
    im = im | shiftdim(surf2volz(node(:, [3 1 2]), fc(:, 1:3), zi, xi, yi), 1);
    im = im | shiftdim(surf2volz(node(:, [2 3 1]), fc(:, 1:3), yi, zi, xi), 2);

    v2smap = [];

    % here we assume the grid is uniform; surf2vol can handle non-uniform grid,
    % but the affine output is not correct in this case

    if (jsonopt('fill', 0, opt) || label)
        im = fillholes3d(im);
        if (label)
            im = cast(im, class(elabel)) * elabel(i);
        end
    end
    img = max(cast(im, class(img)), img);
end

if (nargout > 1)
    dlen = abs([xi(2) - xi(1) yi(2) - yi(1) zi(2) - zi(1)]);
    p0 = min(node);
    offset = p0;
    v2smap = diag(abs(dlen));
    v2smap(4, 4) = 1;
    v2smap(1:3, 4) = offset';
end