1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
|
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));
|
