1
2 //this is an enriched Polyhedron with facet normals
3 #include "PolyhedralSurf.h"
4 #include "PolyhedralSurf_rings.h"
5 #include "compute_normals.h"
6 #include <CGAL/Ridges.h>
7
8 #include <CGAL/Umbilics.h>
9 #include <CGAL/Monge_via_jet_fitting.h>
10 #include <fstream>
11 #include <cassert>
12 #include <boost/program_options.hpp>
13
14 namespace po = boost::program_options;
15
16
17 typedef PolyhedralSurf::Traits Kernel;
18 typedef Kernel::FT FT;
19 typedef Kernel::Point_3 Point_3;
20 typedef Kernel::Vector_3 Vector_3;
21
22 typedef boost::graph_traits<PolyhedralSurf>::vertex_descriptor vertex_descriptor;
23 typedef boost::graph_traits<PolyhedralSurf>::vertex_iterator vertex_iterator;
24 typedef boost::graph_traits<PolyhedralSurf>::face_descriptor face_descriptor;
25
26 typedef T_PolyhedralSurf_rings<PolyhedralSurf> Poly_rings;
27 typedef CGAL::Monge_via_jet_fitting<Kernel> Monge_via_jet_fitting;
28 typedef Monge_via_jet_fitting::Monge_form Monge_form;
29
30
31 typedef std::map<vertex_descriptor, FT> VertexFT_map;
32 typedef boost::associative_property_map< VertexFT_map > VertexFT_property_map;
33
34 typedef std::map<vertex_descriptor, Vector_3> VertexVector_map;
35 typedef boost::associative_property_map< VertexVector_map > VertexVector_property_map;
36
37 typedef std::map<face_descriptor, Vector_3> Face2Vector_map;
38 typedef boost::associative_property_map< Face2Vector_map > Face2Vector_property_map;
39
40 //RIDGES
41 typedef CGAL::Ridge_line<PolyhedralSurf> Ridge_line;
42 typedef CGAL::Ridge_approximation < PolyhedralSurf,
43 VertexFT_property_map,
44 VertexVector_property_map > Ridge_approximation;
45 //UMBILICS
46 typedef CGAL::Umbilic<PolyhedralSurf> Umbilic;
47 typedef CGAL::Umbilic_approximation < PolyhedralSurf,
48 VertexFT_property_map,
49 VertexVector_property_map > Umbilic_approximation;
50
51 //create property maps
52 VertexFT_map vertex_k1_map, vertex_k2_map,
53 vertex_b0_map, vertex_b3_map,
54 vertex_P1_map, vertex_P2_map;
55 VertexVector_map vertex_d1_map, vertex_d2_map;
56 Face2Vector_map face2normal_map;
57
58 VertexFT_property_map vertex_k1_pm(vertex_k1_map), vertex_k2_pm(vertex_k2_map),
59 vertex_b0_pm(vertex_b0_map), vertex_b3_pm(vertex_b3_map),
60 vertex_P1_pm(vertex_P1_map), vertex_P2_pm(vertex_P2_map);
61 VertexVector_property_map vertex_d1_pm(vertex_d1_map), vertex_d2_pm(vertex_d2_map);
62 Face2Vector_property_map face2normal_pm(face2normal_map);
63
64 // default fct parameter values and global variables
65 unsigned int d_fitting = 3;
66 unsigned int d_monge = 3;
67 unsigned int nb_rings = 0;//seek min # of rings to get the required #pts
68 unsigned int nb_points_to_use = 0;//
69 CGAL::Ridge_order tag_order = CGAL::Ridge_order_3;
70 double umb_size = 2;
71 bool verbose = false;
72 unsigned int min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2;
73
74 /* gather points around the vertex v using rings on the
75 polyhedralsurf. the collection of points resorts to 3 alternatives:
76 1. the exact number of points to be used
77 2. the exact number of rings to be used
78 3. nothing is specified
79 */
80 template <typename VertexPointMap>
gather_fitting_points(vertex_descriptor v,std::vector<Point_3> & in_points,Poly_rings & poly_rings,VertexPointMap vpm)81 void gather_fitting_points(vertex_descriptor v,
82 std::vector<Point_3> &in_points,
83 Poly_rings& poly_rings,
84 VertexPointMap vpm)
85 {
86 //container to collect vertices of v on the PolyhedralSurf
87 std::vector<vertex_descriptor> gathered;
88 //initialize
89 in_points.clear();
90
91 //OPTION -p nb_points_to_use, with nb_points_to_use != 0. Collect
92 //enough rings and discard some points of the last collected ring to
93 //get the exact "nb_points_to_use"
94 if ( nb_points_to_use != 0 ) {
95 poly_rings.collect_enough_rings(v, nb_points_to_use, gathered);
96 if ( gathered.size() > nb_points_to_use ) gathered.resize(nb_points_to_use);
97 }
98 else { // nb_points_to_use=0, this is the default and the option -p is not considered;
99 // then option -a nb_rings is checked. If nb_rings=0, collect
100 // enough rings to get the min_nb_points required for the fitting
101 // else collect the nb_rings required
102 if ( nb_rings == 0 )
103 poly_rings.collect_enough_rings(v, min_nb_points, gathered);
104 else poly_rings.collect_i_rings(v, nb_rings, gathered);
105 }
106
107 //store the gathered points
108 std::vector<vertex_descriptor>::const_iterator
109 itb = gathered.begin(), ite = gathered.end();
110 CGAL_For_all(itb,ite) in_points.push_back(get(vpm,*itb));
111 }
112
113 /* Use the jet_fitting package and the class Poly_rings to compute
114 diff quantities.
115 */
compute_differential_quantities(PolyhedralSurf & P,Poly_rings & poly_rings)116 void compute_differential_quantities(PolyhedralSurf& P, Poly_rings& poly_rings)
117 {
118 //container for approximation points
119 std::vector<Point_3> in_points;
120
121 typedef boost::property_map<PolyhedralSurf,CGAL::vertex_point_t>::type VPM;
122 VPM vpm = get(CGAL::vertex_point,P);
123
124 //MAIN LOOP
125 vertex_iterator vitb = P.vertices_begin(), vite = P.vertices_end();
126 for (; vitb != vite; vitb++) {
127 //initialize
128 vertex_descriptor v = * vitb;
129 in_points.clear();
130 Monge_form monge_form;
131 Monge_via_jet_fitting monge_fit;
132
133 //gather points around the vertex using rings
134 gather_fitting_points(v, in_points, poly_rings, vpm);
135
136 //exit if the nb of points is too small
137 if ( in_points.size() < min_nb_points )
138 {std::cerr << "Too few points to perform the fitting" << std::endl; exit(1);}
139
140 //For Ridges we need at least 3rd order info
141 assert( d_monge >= 3);
142 // run the main fct : perform the fitting
143 monge_form = monge_fit(in_points.begin(), in_points.end(),
144 d_fitting, d_monge);
145
146 //switch min-max ppal curv/dir wrt the mesh orientation
147 const Vector_3 normal_mesh = computeFacetsAverageUnitNormal(P,v, face2normal_pm, Kernel());
148 monge_form.comply_wrt_given_normal(normal_mesh);
149
150 //Store monge data needed for ridge computations in property maps
151 vertex_d1_map[v] = monge_form.maximal_principal_direction();
152 vertex_d2_map[v] = monge_form.minimal_principal_direction();
153 vertex_k1_map[v] = monge_form.coefficients()[0];
154 vertex_k2_map[v] = monge_form.coefficients()[1];
155 vertex_b0_map[v] = monge_form.coefficients()[2];
156 vertex_b3_map[v] = monge_form.coefficients()[5];
157 if ( d_monge >= 4) {
158 //= 3*b1^2+(k1-k2)(c0-3k1^3)
159 vertex_P1_map[v] =
160 3*monge_form.coefficients()[3]*monge_form.coefficients()[3]
161 +(monge_form.coefficients()[0]-monge_form.coefficients()[1])
162 *(monge_form.coefficients()[6]
163 -3*monge_form.coefficients()[0]*monge_form.coefficients()[0]
164 *monge_form.coefficients()[0]);
165 //= 3*b2^2+(k2-k1)(c4-3k2^3)
166 vertex_P2_map[v] =
167 3*monge_form.coefficients()[4]*monge_form.coefficients()[4]
168 +(-monge_form.coefficients()[0]+monge_form.coefficients()[1])
169 *(monge_form.coefficients()[10]
170 -3*monge_form.coefficients()[1]*monge_form.coefficients()[1]
171 *monge_form.coefficients()[1]);
172 }
173 } //END FOR LOOP
174 }
175
176
177 ///////////////MAIN///////////////////////////////////////////////////////
178 #if defined(CGAL_USE_BOOST_PROGRAM_OPTIONS) && ! defined(DONT_USE_BOOST_PROGRAM_OPTIONS)
main(int argc,char * argv[])179 int main(int argc, char *argv[])
180 #else
181 int main()
182 #endif
183 {
184 std::string if_name, of_name;// of_name same as if_name with '/' -> '_'
185
186 try {
187 #if defined(CGAL_USE_BOOST_PROGRAM_OPTIONS) && ! defined(DONT_USE_BOOST_PROGRAM_OPTIONS)
188 unsigned int int_tag;
189 po::options_description desc("Allowed options");
190 desc.add_options()
191 ("help,h", "produce help message.")
192 ("input-file,f", po::value<std::string>(&if_name)->default_value("data/poly2x^2+y^2-0.062500.off"),
193 "name of the input off file")
194 ("degree-jet,d", po::value<unsigned int>(&d_fitting)->default_value(3),
195 "degree of the jet, 3 <= degre-jet <= 4")
196 ("degree-monge,m", po::value<unsigned int>(&d_monge)->default_value(3),
197 "degree of the Monge rep, 3<= degree-monge <= degree-jet")
198 ("nb-rings,a", po::value<unsigned int>(&nb_rings)->default_value(0),
199 "number of rings to collect neighbors. 0 means collect enough rings to make appro possible a>=1 fixes the nb of rings to be collected")
200 ("nb-points,p", po::value<unsigned int>(&nb_points_to_use)->default_value(0),
201 "number of neighbors to use. 0 means this option is not considered, this is the default p>=1 fixes the nb of points to be used")
202 ("ridge_order,t", po::value<unsigned int>(&int_tag)->default_value(3),
203 "Order of differential quantities used, must be 3 or 4")
204 ("umbilic-patch-size,u", po::value<double>(&umb_size)->default_value(2),
205 "size of umbilic patches (as multiple of 1ring size)")
206 ("verbose,v", po::value<bool>(&verbose)->default_value(false),
207 "verbose output on text file")
208 ;
209
210 po::variables_map vm;
211 po::store(po::parse_command_line(argc, argv, desc), vm);
212 po::notify(vm);
213
214 if (vm.count("help")) {
215 std::cerr << desc << "\n";
216 return 1;
217 }
218
219
220 if (vm.count("ridge_order")){
221 if ( int_tag == 3 ) tag_order = CGAL::Ridge_order_3;
222 if ( int_tag == 4 ) tag_order = CGAL::Ridge_order_4;
223 if ( int_tag != 3 && int_tag != 4 )
224 {std::cerr << "ridge_order must be CGAL::Ridge_order_3 or CGAL::Ridge_order_4";
225 return 1;}
226 }
227 #else
228 std::cerr << "Command-line options require Boost.ProgramOptions" << std::endl;
229 if_name = "data/poly2x^2+y^2-0.062500.off";
230 d_fitting = 3;
231 d_monge = 3;
232 nb_rings = 0;
233 nb_points_to_use = 0;
234 umb_size = 2;
235 verbose = false;
236 #endif
237 }
238
239 catch(std::exception& e) {
240 std::cerr << "error: " << e.what() << "\n";
241 return 1;
242 }
243 catch(...) {
244 std::cerr << "Exception of unknown type!\n";
245 }
246
247 //modify global variables
248 min_nb_points = (d_fitting + 1) * (d_fitting + 2) / 2;
249
250 //prepare output file names
251 assert(!if_name.empty());
252 of_name = if_name;
253 for(unsigned int i=0; i<of_name.size(); i++)
254 if (of_name[i] == '/') of_name[i]='_';
255 std::ostringstream str_4ogl;
256 str_4ogl << "data/"
257 << of_name << "RIDGES"
258 << "-d" << d_fitting
259 << "-m" << d_monge
260 << "-t" << tag_order
261 << "-a" << nb_rings
262 << "-p" << nb_points_to_use
263 << ".4ogl.txt";
264 std::cout << str_4ogl.str() << std::endl ;
265 std::ofstream out_4ogl(str_4ogl.str().c_str() , std::ios::out);
266
267 //if verbose only...
268 std::ostringstream str_verb;
269 str_verb << "data/"
270 << of_name << "RIDGES"
271 << "-d" << d_fitting
272 << "-m" << d_monge
273 << "-t" << tag_order
274 << "-a" << nb_rings
275 << "-p" << nb_points_to_use
276 << ".verb.txt";
277 std::cout << str_verb.str() << std::endl ;
278 std::ofstream out_verb(str_verb.str().c_str() , std::ios::out);
279
280 //load the model from <mesh.off>
281 PolyhedralSurf P;
282 std::ifstream stream(if_name.c_str());
283 stream >> P;
284 fprintf(stderr, "loadMesh %d Ves %d Facets\n",
285 (int)P.size_of_vertices(), (int)P.size_of_facets());
286 if(verbose)
287 out_verb << "Polysurf with " << P.size_of_vertices()
288 << " vertices and " << P.size_of_facets()
289 << " facets. " << std::endl;
290
291 //exit if not enough points in the model
292 if (min_nb_points > P.size_of_vertices())
293 {std::cerr << "not enough points in the model" << std::endl; exit(1);}
294
295 //initialize Polyhedral data : normal of facets
296 compute_facets_normals(P,face2normal_pm, Kernel());
297
298 //create a Poly_rings object
299 Poly_rings poly_rings(P);
300
301 std::cout << "Compute differential quantities via jet fitting..." << std::endl;
302 //initialize the diff quantities property maps
303 compute_differential_quantities(P, poly_rings);
304
305 //---------------------------------------------------------------------------
306 //Ridges
307 //--------------------------------------------------------------------------
308 std::cout << "Compute ridges..." << std::endl;
309 Ridge_approximation ridge_approximation(P,
310 vertex_k1_pm, vertex_k2_pm,
311 vertex_b0_pm, vertex_b3_pm,
312 vertex_d1_pm, vertex_d2_pm,
313 vertex_P1_pm, vertex_P2_pm );
314 std::vector<Ridge_line*> ridge_lines;
315 std::back_insert_iterator<std::vector<Ridge_line*> > ii(ridge_lines);
316
317 //Find MAX_RIDGE, MIN_RIDGE, CREST_RIDGES
318 // ridge_approximation.compute_max_ridges(ii, tag_order);
319 // ridge_approximation.compute_min_ridges(ii, tag_order);
320 ridge_approximation.compute_crest_ridges(ii, tag_order);
321
322 // or with the global function
323 CGAL::compute_max_ridges(P,
324 vertex_k1_pm, vertex_k2_pm,
325 vertex_b0_pm, vertex_b3_pm,
326 vertex_d1_pm, vertex_d2_pm,
327 vertex_P1_pm, vertex_P2_pm,
328 ii, tag_order);
329
330 std::vector<Ridge_line*>::iterator iter_lines = ridge_lines.begin(),
331 iter_end = ridge_lines.end();
332 //OpenGL output
333
334 typedef boost::property_map<PolyhedralSurf,CGAL::vertex_point_t>::type VPM;
335 VPM vpm = get(CGAL::vertex_point,P);
336
337 for (;iter_lines!=iter_end;iter_lines++) (*iter_lines)->dump_4ogl(out_4ogl, vpm);
338
339
340 for (iter_lines = ridge_lines.begin();iter_lines!=iter_end;iter_lines++){
341 //verbose txt output
342 if (verbose){
343 (*iter_lines)->dump_verbose(out_verb,vpm);
344 }
345 delete *iter_lines;
346 }
347
348 //---------------------------------------------------------------------------
349 // UMBILICS
350 //--------------------------------------------------------------------------
351 std::cout << "Compute umbilics..." << std::endl;
352 std::vector<Umbilic*> umbilics;
353 std::back_insert_iterator<std::vector<Umbilic*> > umb_it(umbilics);
354
355 //explicit construction of the class
356 // Umbilic_approximation umbilic_approximation(P,
357 // vertex_k1_pm, vertex_k2_pm,
358 // vertex_d1_pm, vertex_d2_pm);
359 // umbilic_approximation.compute(umb_it, umb_size);
360 //or global function call
361 CGAL::compute_umbilics(P,
362 vertex_k1_pm, vertex_k2_pm,
363 vertex_d1_pm, vertex_d2_pm,
364 umb_it, umb_size);
365
366 std::vector<Umbilic*>::iterator iter_umb = umbilics.begin(),
367 iter_umb_end = umbilics.end();
368 // output
369 std::cout << "nb of umbilics " << umbilics.size() << std::endl;
370 for (;iter_umb!=iter_umb_end;iter_umb++) std::cout << **iter_umb;
371
372 //verbose txt output
373 if (verbose) {
374 out_verb << "nb of umbilics " << umbilics.size() << std::endl;
375 }
376 for ( iter_umb = umbilics.begin();iter_umb!=iter_umb_end;iter_umb++){
377 if (verbose) {
378 out_verb << **iter_umb;
379 }
380 delete *iter_umb;
381 }
382 return 0;
383 }
384
385
386