vol2restrictedtri.m

function [node, elem] = vol2restrictedtri(vol, thres, cent, brad, ang, radbound, distbound, maxnode)
%
% [node,elem]=vol2restrictedtri(vol,thres,cent,brad,ang,radbound,distbound,maxnode)
%
% surface mesh extraction using CGAL mesher
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2009/01/06
%
% input:
%       vol: a 3D volumetric image
%       thres: a scalar as the threshold of of the extraction
%       cent: a 3d position (x,y,z) which locates inside the resulting
%             mesh, this is automatically computed from vol2surf
%       brad: maximum bounding sphere squared of the resulting mesh
%       ang: minimum angular constrains of the resulting tranglar elements
%            (in degrees)
%       radbound: maximum triangle delaunay circle radius
%       distbound: maximum delaunay sphere distances
%       maxnode: maximum number of surface nodes (even radbound is not reached)
% output:
%       node: the list of 3d nodes in the resulting surface (x,y,z)
%       elem: the element list of the resulting mesh (3 columns of integers)
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

if (radbound < 1)
    warning(['You are meshing the surface with sub-pixel size. If this ' ...
             'is not your your intent, please check if you set ' ...
             '"opt.radbound" correctly for the default meshing method.']);
end

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

saveinr(vol, mwpath('pre_extract.inr'));
deletemeshfile(mwpath('post_extract.off'));

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

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

initnum = 50;
if (~isempty(getvarfrom({'caller', 'base'}, 'ISO2MESH_INITSIZE')))
    initnum = getvarfrom({'caller', 'base'}, 'ISO2MESH_INITSIZE');
end

system([' "' mcpath('cgalsurf') exesuff '" "' mwpath('pre_extract.inr') ...
        '" ' sprintf('%.16f %.16f %.16f %.16f %.16f %.16f %.16f %.16f %d ', thres, cent, brad, ang, radbound, distbound, maxnode) ...
        ' "' mwpath('post_extract.off') '" ' sprintf('%.0f %d', randseed, initnum)]);
[node, elem] = readoff(mwpath('post_extract.off'));

[node, elem] = meshcheckrepair(node, elem);

% assuming the origin [0 0 0] is located at the lower-bottom corner of the image
node = node + 0.5;