vol2mesh.m

function [node, elem, face, regions] = vol2mesh(img, ix, iy, iz, opt, maxvol, dofix, method, isovalues)
%
% [node,elem,face,regions]=vol2mesh(img,ix,iy,iz,opt,maxvol,dofix,method,isovalues)
%
% convert a binary (or multi-valued) volume to tetrahedral mesh
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    img: a volumetric binary image
%    ix,iy,iz: subvolume selection indices in x,y,z directions
%    opt: as defined in vol2surf.m
%    maxvol: target maximum tetrahedral elem volume
%                when method='cgalmesh', maxvol can specify the target
%                for each label (subregion index) by the following syntax
%                'label1=size1:label2=size2:...'
%    dofix: 1: perform mesh validation&repair, 0: skip repairing
%    method: 'cgalsurf' or omit: use CGAL surface mesher
%        'simplify': use binsurface and then simplify
%        'cgalmesh': use CGAL 3.5 3D mesher for direct mesh generation [new]
%
%        generally speaking, 'cgalmesh' is the most robust path
%        if you want to product meshes from binary or multi-region
%        volumes, however, its limitations include 1) only accept
%        uint8 volume, and 2) can not extract meshes from gray-scale
%        volumes. If ones goal is to process a gray-scale volume,
%        he/she should use the 'cgalsurf' option. 'simplify' approach
%        is not recommended unless other options has failed.
%    isovalues: a list of isovalues where the levelset is defined
%
% output:
%    node: output, node coordinates of the tetrahedral mesh
%    elem: output, element list of the tetrahedral mesh, the last
%          column is the region ID
%    face: output, mesh surface element list of the tetrahedral mesh
%          the last column denotes the boundary ID
%        region: optional output. if opt.autoregion is set to 1, region
%              saves the interior points for each closed surface component
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if (nargin >= 8)
    if (strcmp(method, 'cgalmesh'))
        vol = img(ix, iy, iz);
        if (length(unique(vol(:))) > 64 && dofix == 1)
            error(['it appears that you are processing a ' ...
                   'grayscale image. Currently cgalmesher ' ...
                   'does not support grayscale images. ' ...
                   'Please use "cgalsurf" method to mesh a grayscale ' ...
                   'volume. If you are certain to run cgalmesher ' ...
                   'on your data, please set dofix=0 and run this again.']);
        end
        [node, elem, face] = cgalv2m(vol, opt, maxvol);
        return
    end
end

% first, convert the binary volume into isosurfaces
if (nargin == 8)
    [no, el, regions, holes] = vol2surf(img, ix, iy, iz, opt, dofix, method);
elseif (nargin == 9)
    [no, el, regions, holes] = vol2surf(img, ix, iy, iz, opt, dofix, method, isovalues);
else
    [no, el, regions, holes] = vol2surf(img, ix, iy, iz, opt, dofix, 'cgalsurf');
end
% then, create volumetric mesh from the surface mesh
if (nargin >= 8)
    if (strcmp(method, 'cgalpoly'))
        [node, elem, face] = cgals2m(no(:, 1:3), el(:, 1:3), opt, maxvol);
        return
    end
end

[node, elem, face] = surf2mesh(no, el, [], [], 1, maxvol, regions, holes);