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
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
|
#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
#include <CGAL/Mesh_triangulation_3.h>
#include <CGAL/Mesh_complex_3_in_triangulation_3.h>
#include <CGAL/Mesh_criteria_3.h>
#include <CGAL/Labeled_mesh_domain_3.h>
#include <CGAL/Mesh_domain_with_polyline_features_3.h>
#include <CGAL/make_mesh_3.h>
#include <CGAL/Image_3.h>
// start PV //
#include <string>
#include <vector>
#include <CGAL/Mesh_constant_domain_field_3.h>
// end PV //
// Domain
struct K: public CGAL::Exact_predicates_inexact_constructions_kernel {};
typedef CGAL::Image_3 Image;
typedef CGAL::Labeled_mesh_domain_3<K> Image_domain;
typedef CGAL::Mesh_domain_with_polyline_features_3<Image_domain> Mesh_domain;
// Triangulation
typedef CGAL::Mesh_triangulation_3<Mesh_domain>::type Tr;
typedef CGAL::Mesh_complex_3_in_triangulation_3<Tr> C3t3;
// Mesh Criteria
typedef CGAL::Mesh_criteria_3<Tr> Mesh_criteria;
typedef Mesh_criteria::Facet_criteria Facet_criteria;
typedef Mesh_criteria::Cell_criteria Cell_criteria;
typedef Mesh_criteria::Edge_criteria Edge_criteria;
// start PV //
typedef CGAL::Mesh_constant_domain_field_3<Mesh_domain::R,
Mesh_domain::Index> Sizing_field_cell;
// end PV //
#include <CGAL/Mesh_3/polylines_to_protect.h> // undocumented header
bool add_1D_features(const CGAL::Image_3& image,
Mesh_domain& domain)
{
typedef K::Point_3 Point_3;
typedef unsigned char Word_type;
std::vector<std::vector<Point_3> > polylines_on_bbox;
CGAL::polylines_to_protect<Point_3, Word_type>(image, polylines_on_bbox);
domain.add_features(polylines_on_bbox.begin(), polylines_on_bbox.end());
return true;
}
void usage(char *exename){
printf("usage:\n\t%s input.inr output.mesh <angle|30> <surf-size|6> <approx|4> <rad-edge-ratio|3> \
<tetra-size|8> <randomseed|-1>\n\
example:\n\t%s input.inr output.mesh 30 6 4 3 8 123456789\n",exename,exename);
exit(1);
}
int main(int argc,char *argv[])
{
// Loads image
// start PV //
//float angle=30.f,ssize=6.f,approx=4.f,reratio=3.f,vsize=8.f;
float angle=30.f,ssize=6.f,approx=4.f,reratio=3.f,maxvol=0.f;
// end PV //
int labelid=0, lid;
printf("Volume/Surface Mesh Generation Utility (Based on CGAL %s)\n\
(modified for iso2mesh by Qianqian Fang and Peter Varga)\nhttp://iso2mesh.sf.net\n\n",CGAL_VERSION_STR);
if(argc!=3&&argc!=8&&argc!=9){
usage(argv[0]);
}
if(argc>=8){
angle=atof(argv[3]);
ssize=atof(argv[4]);
approx=atof(argv[5]);
reratio=atof(argv[6]);
// start PV //
//vsize=atof(argv[7]);
// end PV //
}
if(argc==9 && atoi(argv[8])>0){
printf("RNG seed=%d\n",atoi(argv[8]));
CGAL::Random rd(atoi(argv[8]));
CGAL::Random::State st;
rd.save_state(st);
CGAL::get_default_random().restore_state(st);
}
Image image;
image.read(argv[1]);
Mesh_domain domain = Mesh_domain::create_labeled_image_mesh_domain(image);
add_1D_features(image, domain);
// Mesh criteria
Facet_criteria facet_criteria(angle, ssize, approx); // angle, size, approximation
Edge_criteria edge_criteria(ssize);
// start PV //
//Cell_criteria cell_criteria(reratio, vsize); // radius-edge ratio, size
std::string vsize;
vsize=argv[7];
std::vector<double> vsize_vect;
std::vector<int> labels_vect;
std::stringstream stream_vsize(vsize);
std::string word_vsize;
while( std::getline(stream_vsize, word_vsize, ':') ){
int len=sscanf(word_vsize.c_str(), "%d=%f", &lid, &maxvol);
if(len==2){
labelid=lid;
if(maxvol<=0.f){
std::cerr << "cell volume must be positive" << std::endl;
exit(-2);
}
labels_vect.push_back(labelid);
vsize_vect.push_back(maxvol);
}else{
len=sscanf(word_vsize.c_str(), "%f", &maxvol);
if(len!=1 || maxvol<0 || word_vsize.find_first_of('=')!=std::string::npos){
std::cerr << "invalid sizing field label '" << word_vsize << "', please check your command" << std::endl;
exit(-1);
}
labels_vect.push_back(++labelid);
vsize_vect.push_back(maxvol);
}
}
float max_vsize = *max_element(vsize_vect.begin(),vsize_vect.end());
Sizing_field_cell vsize_cell(max_vsize);
int volume_dimension = 3;
for (int i=0; i<vsize_vect.size(); i++)
vsize_cell.set_size(vsize_vect[i], volume_dimension,
domain.index_from_subdomain_index(labels_vect[i]));
std::cout << "Mesh sizes are (label=size) ";
for (int i=0; i<vsize_vect.size(); i++){
std::cout << "(" << labels_vect[i] << "=" << vsize_vect[i] << ") ";
}
std::cout << std::endl;
Cell_criteria cell_criteria(reratio, vsize_cell);
// end PV //
Mesh_criteria criteria(edge_criteria, facet_criteria, cell_criteria);
// Meshing
C3t3 c3t3 = CGAL::make_mesh_3<C3t3>(domain, criteria);
// Output
std::ofstream medit_file;
if(argc>=8)
medit_file.open(argv[2]);
else
medit_file.open("output.mesh");
c3t3.output_to_medit(medit_file);
medit_file.close();
return 0;
}
|
