surf2mesh.m

function [node, elem, face] = surf2mesh(v, f, p0, p1, keepratio, maxvol, regions, holes, forcebox, method, cmdopt)
%
% [node,elem,face]=surf2mesh(v,f,p0,p1,keepratio,maxvol,regions,holes,forcebox)
%
% create quality volumetric mesh from isosurface patches
%
% author: Qianqian Fang, <q.fang at neu.edu>
% date: 2007/11/24
%
% input parameters:
%      v: input, isosurface node list, dimension (nn,3)
%         if v has 4 columns, the last column specifies mesh density near each node
%      f: input, isosurface face element list, dimension (be,3)
%      p0: input, coordinates of one corner of the bounding box, p0=[x0 y0 z0]
%      p1: input, coordinates of the other corner of the bounding box, p1=[x1 y1 z1]
%      keepratio: input, percentage of elements being kept after the simplification
%      maxvol: input, maximum tetrahedra element volume
%      regions: list of regions, specifying by an internal point for each region
%      holes: list of holes, similar to regions
%      forcebox: 1: add bounding box, 0: automatic
%
% outputs:
%      node: output, node coordinates of the tetrahedral mesh
%      elem: output, element list of the tetrahedral mesh
%      face: output, mesh surface element list of the tetrahedral mesh
%             the last column denotes the boundary ID
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

fprintf(1, 'generating tetrahedral mesh from closed surfaces ...\n');

if (nargin < 10)
    method = 'tetgen';
end

exesuff = getexeext;
exesuff = fallbackexeext(exesuff, method);

if (keepratio > 1 || keepratio < 0)
    warn(['The "keepratio" parameter is required to be between 0 and 1. '...
          'Your input is out of this range. surf2mesh will not perform '...
          'simplification. Please double check to correct this.']);
end

% first, resample the surface mesh with cgal
if (keepratio < 1 - 1e-9 && ~iscell(f))
    fprintf(1, 'resampling surface mesh ...\n');
    [no, el] = meshresample(v(:, 1:3), f(:, 1:3), keepratio);
    el = unique(sort(el, 2), 'rows');

    % then smooth the resampled surface mesh (Laplace smoothing)

    %% edges=surfedge(el);  % disable on 12/05/08, very slow on octave
    %% mask=zeros(size(no,1),1);
    %% mask(unique(edges(:)))=1;  % =1 for edge nodes, =0 otherwise
    % [conn,connnum,count]=meshconn(el,length(no));
    % no=smoothsurf(no,mask,conn,2);

    % remove end elements (all nodes are edge nodes)
    % el=delendelem(el,mask);
else
    no = v;
    el = f;
end
if (nargin == 6)
    regions = [];
    holes = [];
elseif (nargin == 7)
    holes = [];
end

if (size(regions, 2) >= 4 && ~isempty(maxvol))
    warning('you specified both maxvol and the region based volume constraint,the maxvol setting will be ignored');
    maxvol = [];
end

dobbx = 0;
if (nargin >= 9)
    dobbx = forcebox;
end

% dump surface mesh to .poly file format
if (~iscell(el) && ~isempty(no) && ~isempty(el))
    saveoff(no(:, 1:3), el(:, 1:3), mwpath('post_vmesh.off'));
end
deletemeshfile(mwpath('post_vmesh.mtr'));
savesurfpoly(no, el, holes, regions, p0, p1, mwpath('post_vmesh.poly'), dobbx);

moreopt = '';
if (size(no, 2) == 4)
    moreopt = [moreopt ' -m '];
end
% call tetgen to create volumetric mesh
deletemeshfile(mwpath('post_vmesh.1.*'));
fprintf(1, 'creating volumetric mesh from a surface mesh ...\n');

if (nargin < 11)
    try
        cmdopt = evalin('caller', 'ISO2MESH_TETGENOPT');
    catch
        try
            cmdopt = evalin('base', 'ISO2MESH_TETGENOPT');
        catch
            cmdopt = '';
        end
    end
end

if (isempty(cmdopt))
    [status, cmdout] = system([' "' mcpath(method, exesuff) '" -A -q1.414a' num2str(maxvol) ' ' moreopt ' "' mwpath('post_vmesh.poly') '"']);
else
    [status, cmdout] = system([' "' mcpath(method, exesuff) '" ' cmdopt ' "' mwpath('post_vmesh.poly') '"']);
end

if (status ~= 0)
    error('Tetgen command failed:\n%s\n', cmdout);
end

% read in the generated mesh
[node, elem, face] = readtetgen(mwpath('post_vmesh.1'));

fprintf(1, 'volume mesh generation is complete\n');