mesh2vol.m

function [mask weight] = mesh2vol(node, elem, xi, yi, zi)
%
% [mask weight]=mesh2vol(node,face,Nxyz)
% [mask weight]=mesh2vol(node,face,[Nx,Ny,Nz])
% [mask weight]=mesh2vol(node,face,xi,yi,zi,hf)
%   or
% newval=mesh2vol(node_val,face,...)
%
% fast rasterization of a 3D mesh to a volume with tetrahedron index labels
%
% author: Qianqian Fang <q.fang at neu.edu>
% date for initial version: Feb 10,2014
%
% input:
%      node: node coordinates, dimension N by 2 or N by 3 array
%      nodeval: a 4-column array, the first 3 columns are the node coordinates,
%            the last column denotes the values associated with each node
%      face: a triangle surface, N by 3 or N by 4 array
%      Nx,Ny,Nxy: output image in x/y dimensions, or both
%      xi,yi: linear vectors for the output pixel center positions in x/y
%      hf: the handle of a pre-created figure window for faster rendering
%
% output:
%      mask: a 3D image, the value of each pixel is the index of the
%            enclosing triangle, if the pixel is outside of the mesh, NaN
%      weight: (optional) a 3 by Nx by Ny array, where Nx/Ny are the dimensions
%            for the mask
%      newval: when the node has 4 columns, the last column represents the
%            values (color) at each node, the output newval is the rasterized
%            mesh value map over the specified grid.
%
% note: This function only works for matlab
%
% example:
%
%   [no,el]=meshgrid6(0:5,0:5,0:3);
%   mask=mesh2vol(no,el(:,1:4),0:0.1:5,0:0.1:5,0:0.1:4);
%   imagesc(mask(:,:,3))
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

nodeval = [];
if (size(node, 2) == 4)
    nodeval = node(:, 4);
    node(:, 4) = [];
end

if (nargin == 3 && length(xi) == 1 && xi > 0)
    mn = min(node);
    mx = max(node);
    df = (mx(1:min(3, size(node, 2))) - mn(1:min(3, size(node, 2)))) / xi;
elseif (nargin == 3 && length(xi) == 3 && all(xi > 0))
    mn = min(node);
    mx = max(node);
    df = (mx(1:min(3, size(node, 2))) - mn(1:min(3, size(node, 2)))) ./ (xi(:)');
elseif (nargin == 5)
    mx = [max(xi) max(yi) max(zi)];
    mn = [min(xi) min(yi) min(zi)];
    df = [min(diff(xi(:))) min(diff(yi(:))) min(diff(zi(:)))];
else
    error('you must give at least xi input');
end

xi = mn(1):df(1):mx(1);
yi = mn(2):df(2):mx(2);
zi = mn(3):df(3):mx(3);

if (size(node, 2) ~= 3 || size(elem, 2) <= 3)
    error('node must have 3 columns; face can not have less than 3 columns');
end

nz = length(zi);
mask = zeros(length(xi) - 1, length(yi) - 1, length(zi) - 1);
if (nargout > 1 || ~isempty(nodeval))
    weight = zeros(4, length(xi) - 1, length(yi) - 1, length(zi) - 1);
end

hf = figure('visible', 'on');
for i = 1:nz
    if (~isempty(nodeval))
        [cutpos, cutvalue, facedata, elemid] = qmeshcut(elem, node, nodeval, ['z=' num2str(zi(i))]);
    else
        [cutpos, cutvalue, facedata, elemid] = qmeshcut(elem, node, node(:, 1), ['z=' num2str(zi(i))]);
    end
    if (isempty(cutpos))
        continue
    end
    if (nargout > 1 || ~isempty(nodeval))
        [maskz, weightz] = mesh2mask(cutpos, facedata, xi, yi, hf);
        weight(:, :, :, i) = weightz;
    else
        maskz = mesh2mask(cutpos, facedata, xi, yi, hf);
    end
    idx = find(~isnan(maskz));
    if (~isempty(nodeval))
        eid = facedata(maskz(idx), :);
        maskz(idx) = cutvalue(eid(:, 1)) .* weightz(1, idx)' + cutvalue(eid(:, 2)) .* weightz(2, idx)' + ...
            cutvalue(eid(:, 3)) .* weightz(3, idx)' + cutvalue(eid(:, 4)) .* weightz(4, idx)';
    else
        maskz(idx) = elemid(maskz(idx));
    end
    mask(:, :, i) = maskz;
end
close(hf);