surf2volz.m

function img = surf2volz(node, face, xi, yi, zi)
%
% img=surf2volz(node,face,xi,yi,zi)
%
% convert a triangular surface to a shell of voxels in a 3D image
% along the z-axis
%
% 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
%    xi,yi,zi: x/y/z grid for the resulting volume
%
% output:
%    img: a volumetric binary image at position of ndgrid(xi,yi,zi)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

ne = size(face, 1);
img = zeros(length(xi), length(yi), length(zi), 'uint8');

dx0 = min(abs(diff(xi)));
dx = dx0 / 2;
dy0 = min(abs(diff(yi)));
dy = dy0 / 2;
dz0 = min(abs(diff(zi)));

dl = sqrt(dx * dx + dy * dy);

minz = min(node(:, 3));
maxz = max(node(:, 3));
iz = hist([minz, maxz], zi);
hz = find(iz);
iz = hz(1):min(length(zi), hz(end) + 1);

for i = 1:length(iz)
    plane = [0 100 zi(iz(i)); 100 0 zi(iz(i)); 0 0 zi(iz(i))];
    [bcutpos, bcutvalue, bcutedges] = qmeshcut(face(:, 1:3), node, node(:, 1), plane);
    if (isempty(bcutpos))
        continue
    end
    enum = size(bcutedges, 1);
    for j = 1:enum
        e0 = bcutpos(bcutedges(j, 1), 1:2);
        e1 = bcutpos(bcutedges(j, 2), 1:2);
        len = ceil(sum(abs(e1 - e0)) / (abs(dx) + abs(dy))) + 1;
        dd = (e1 - e0) / len;

        posx = floor(((e0(1) + (0:len) * dd(1) - xi(1))) / dx0)';
        posy = floor(((e0(2) + (0:len) * dd(2) - yi(1))) / dy0)';
        pos = [posx, posy];
        pos(find(posx > length(xi) | posy > length(yi) | posx <= 0 | posy <= 0), :) = [];

        if (length(pos) > 0)
            zz = floor(((zi(iz(i)) - zi(1))) / dz0);
            for k = 1:size(pos, 1)
                img(pos(k, 1), pos(k, 2), zz) = 1;
            end
            % img(sub2ind(size(img),pos(:,1),pos(:,2),i*ones(size(pos,1),1)))=1;
        end
    end
end