cgalv2m.m

function [node, elem, face] = cgalv2m(vol, opt, maxvol)
%
% [node,elem,face]=cgalv2m(vol,opt,maxvol)
%
% wrapper for CGAL 3D mesher (CGAL 3.5 or up)
% convert a binary (or multi-valued) volume to tetrahedral mesh
%
% http://www.cgal.org/Manual/3.5/doc_html/cgal_manual/Mesh_3/Chapter_main.html
%
% author: Qianqian Fang (q.fang at neu.edu)
%
% input:
%    vol: a volumetric binary image
%    ix,iy,iz: subvolume selection indices in x,y,z directions
%    opt: parameters for CGAL mesher, if opt is a structure, then
%        opt.radbound: defines the maximum surface element size
%        opt.angbound: defines the miminum angle of a surface triangle
%        opt.distbound: defines the maximum distance between the
%        center of the surface bounding circle and center of the
%        element bounding sphere
%        opt.reratio:  maximum radius-edge ratio
%        if opt is a scalar, it only specifies radbound.
%    maxvol: target maximum tetrahedral elem volume
%
% 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
%         note: each triangle will appear twice in the face list with each
%           one attaches to each side of the interface. one can remove
%           the redundant triangles by unique(face(:,1:3),'rows')
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fprintf(1, 'creating surface and tetrahedral mesh from a multi-domain volume ...\n');

if (~(islogical(vol) || isa(vol, 'uint8')))
    error('cgalmesher can only handle uint8 volumes, you have to convert your image to unit8 first.');
end

if (~any(vol))
    error('no labeled regions found in the input volume.');
end

exesuff = getexeext;
exesuff = fallbackexeext(exesuff, 'cgalmesh');

ang = 30;
ssize = 6;
approx = 0.5;
reratio = 3;

if (~isstruct(opt))
    ssize = opt;
end

if (isstruct(opt) && length(opt) == 1)  % does not support settings for multiple labels
    ssize = jsonopt('radbound', ssize, opt);
    ang = jsonopt('angbound', ang, opt);
    approx = jsonopt('distbound', approx, opt);
    reratio = jsonopt('reratio', reratio, opt);
end

saveinr(vol, mwpath('pre_cgalmesh.inr'));
deletemeshfile(mwpath('post_cgalmesh.mesh'));

randseed = hex2dec('623F9A9E'); % "U+623F U+9A9E"

if (~isempty(getvarfrom('base', 'ISO2MESH_RANDSEED')))
    randseed = getvarfrom('base', 'ISO2MESH_RANDSEED');
end

if (ischar(maxvol))
    format_maxvol = '%s';
else
    format_maxvol = '%f';
end
cmd = sprintf(['"%s%s" "%s" "%s" %f %f %f %f ' format_maxvol ' %d'], mcpath('cgalmesh'), exesuff, ...
              mwpath('pre_cgalmesh.inr'), mwpath('post_cgalmesh.mesh'), ang, ssize, ...
              approx, reratio, maxvol, randseed);
system(cmd);
if (~exist(mwpath('post_cgalmesh.mesh'), 'file'))
    error(['output file was not found, failure was encountered when running command: \n', cmd]);
end
[node, elem, face] = readmedit(mwpath('post_cgalmesh.mesh'));

% if a transformation matrix/offset vector supplied, apply them
if (isstruct(opt) && length(opt) == 1)
    if (isfield(opt, 'A') && isfield(opt, 'B'))
        node(:, 1:3) = (opt.A * node(:, 1:3)' + repmat(opt.B(:), 1, size(node, 1)))';
    end
end

fprintf(1, 'node number:\t%d\ntriangles:\t%d\ntetrahedra:\t%d\nregions:\t%d\n', ...
        size(node, 1), size(face, 1), size(elem, 1), length(unique(elem(:, end))));
fprintf(1, 'surface and volume meshes complete\n');

if (size(node, 1) > 0)
    [node, elem, face] = sortmesh(node(1, :), node, elem, 1:4, face, 1:3);
end

node = node + 0.5;
elem(:, 1:4) = meshreorient(node(:, 1:3), elem(:, 1:4));