Developer Documentation

adaptive_subdivider.cc
1 /* ========================================================================= *
2 * *
3 * OpenMesh *
4 * Copyright (c) 2001-2015, RWTH-Aachen University *
5 * Department of Computer Graphics and Multimedia *
6 * All rights reserved. *
7 * www.openmesh.org *
8 * *
9 *---------------------------------------------------------------------------*
10 * This file is part of OpenMesh. *
11 *---------------------------------------------------------------------------*
12 * *
13 * Redistribution and use in source and binary forms, with or without *
14 * modification, are permitted provided that the following conditions *
15 * are met: *
16 * *
17 * 1. Redistributions of source code must retain the above copyright notice, *
18 * this list of conditions and the following disclaimer. *
19 * *
20 * 2. Redistributions in binary form must reproduce the above copyright *
21 * notice, this list of conditions and the following disclaimer in the *
22 * documentation and/or other materials provided with the distribution. *
23 * *
24 * 3. Neither the name of the copyright holder nor the names of its *
25 * contributors may be used to endorse or promote products derived from *
26 * this software without specific prior written permission. *
27 * *
28 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
29 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
30 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A *
31 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
32 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
33 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
34 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
35 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
36 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
37 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
38 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
39 * *
40 * ========================================================================= */
41
42
43
44 // -------------------------------------------------------------- includes ----
45
46 // -------------------- OpenMesh
47 #include <OpenMesh/Core/IO/MeshIO.hh>
48 #include <OpenMesh/Core/Mesh/TriMesh_ArrayKernelT.hh>
49 #include <OpenMesh/Tools/Utils/Timer.hh>
50 #include <OpenMesh/Tools/Utils/getopt.h>
51 // -------------------- OpenMesh Adaptive Composite Subdivider
52 #include <OpenMesh/Tools/Subdivider/Adaptive/Composite/CompositeT.hh>
53 #include <OpenMesh/Tools/Subdivider/Adaptive/Composite/RulesT.hh>
54 // -------------------- STL
55 #include <iostream>
56 #include <fstream>
57 #include <sstream>
58 #include <limits>
59 #if defined(OM_CC_MIPS)
60 # include <math.h>
61 #else
62 # include <cmath>
63 using std::pow;
64 #endif
65
66
67 using OpenMesh::Subdivider::Adaptive::CompositeTraits;
68
69 // define mesh, rule interface, and subdivider types
70 typedef OpenMesh::TriMesh_ArrayKernelT<CompositeTraits> MyMesh;
71 typedef OpenMesh::Subdivider::Adaptive::RuleInterfaceT<MyMesh> Rule;
72 typedef OpenMesh::Subdivider::Adaptive::CompositeT<MyMesh> Subdivider;
73
74 // ----------------------------------------------------------------------------
75
76 using namespace OpenMesh::Subdivider;
77
78 // factory function to add a RULE to a subdivider
79 #define ADD_FN( RULE ) \
80 bool add_ ## RULE( Subdivider& _sub ) \
81 { return _sub.add< Adaptive:: RULE < MyMesh > >(); }
82
83 ADD_FN( Tvv3 );
84 ADD_FN( Tvv4 );
85 ADD_FN( VF );
86 ADD_FN( FF );
87 ADD_FN( FFc );
88 ADD_FN( FV );
89 ADD_FN( FVc );
90 ADD_FN( VV );
91 ADD_FN( VVc );
92 ADD_FN( VE );
93 ADD_FN( VdE );
94 ADD_FN( VdEc );
95 ADD_FN( EV );
96 ADD_FN( EVc );
97 ADD_FN( EF );
98 ADD_FN( FE );
99 ADD_FN( EdE );
100 ADD_FN( EdEc );
101
102 #undef ADD_FN
103
104 typedef bool (*add_rule_ft)( Subdivider& );
105
106 // map rule name to factory function
107 struct RuleMap : std::map< std::string, add_rule_ft >
108 {
109 RuleMap()
110 {
111 #define ADD( RULE ) \
112 (*this)[ #RULE ] = add_##RULE;
113
114 ADD( Tvv3 );
115 ADD( Tvv4 );
116 ADD( VF );
117 ADD( FF );
118 ADD( FFc );
119 ADD( FV );
120 ADD( FVc );
121 ADD( VV );
122 ADD( VVc );
123 ADD( VE );
124 ADD( VdE );
125 ADD( VdEc );
126 ADD( EV );
127 ADD( EVc );
128 ADD( EF );
129 ADD( FE );
130 ADD( EdE );
131 ADD( EdEc );
132
133 #undef ADD
134 }
135
136 } available_rules;
137
138
139 // ----------------------------------------------------------------------------
140
141 std::string basename( const std::string& _fname );
142 void usage_and_exit(const std::string& _fname, int xcode);
143
144 // ----------------------------------------------------------------------------
145
146
147 int main(int argc, char **argv)
148 {
149 size_t n_iter = 0; // n iteration
150 size_t max_nv = std::numeric_limits<size_t>::max(); // max. number of vertices in the end
151 std::string ifname; // input mesh
152 std::string ofname; // output mesh
153 std::string rule_sequence = "Tvv3 VF FF FVc"; // sqrt3 default
154 bool uniform = false;
155 int c;
156
157 // ---------------------------------------- evaluate command line
158 while ( (c=getopt(argc, argv, "hlm:n:r:sU"))!=-1 )
159 {
160 switch(c)
161 {
162 case 's': rule_sequence = "Tvv3 VF FF FVc"; break; // sqrt3
163 case 'l': rule_sequence = "Tvv4 VdE EVc VdE EVc"; break; // loop
164 case 'n': { std::stringstream s; s << optarg; s >> n_iter; } break;
165 case 'm': { std::stringstream s; s << optarg; s >> max_nv; } break;
166 case 'r': rule_sequence = optarg; break;
167 case 'U': uniform = true; break;
168 case 'h': usage_and_exit(argv[0],0); break;
169 case '?':
170 default: usage_and_exit(argv[0],1);
171 }
172 }
173
174 if ( optind == argc )
175 usage_and_exit(argv[0],2);
176
177 if ( optind < argc )
178 ifname = argv[optind++];
179
180 if ( optind < argc )
181 ofname = argv[optind++];
182
183 // if ( optind < argc ) // too many arguments
184
185 // ---------------------------------------- mesh and subdivider
186 MyMesh mesh;
187 Subdivider subdivider(mesh);
188
189
190 // -------------------- read mesh from file
191 std::cout << "Input mesh : " << ifname << std::endl;
192 if (!OpenMesh::IO::read_mesh(mesh, ifname))
193 {
194 std::cerr << " Error reading file!\n";
195 return 1;
196 }
197
198 // store orignal size of mesh
199 size_t n_vertices = mesh.n_vertices();
200 size_t n_edges = mesh.n_edges();
201 size_t n_faces = mesh.n_faces();
202
203 if ( n_iter > 0 )
204 std::cout << "Desired #iterations: " << n_iter << std::endl;
205
206 if ( max_nv < std::numeric_limits<size_t>::max() )
207 {
208 std::cout << "Desired max. #V : " << max_nv << std::endl;
209 if (!n_iter )
210 n_iter = std::numeric_limits<size_t>::max();
211 }
212
213
214 // -------------------- Setup rule sequence
215 {
216 std::stringstream s;
217 std::string token;
218
219 RuleMap::iterator it = available_rules.end();
220
221 for (s << rule_sequence; s >> token; )
222 {
223 if ( (it=available_rules.find( token )) != available_rules.end() )
224 {
225 it->second( subdivider );
226 }
227 else if ( token[0]=='(' && (subdivider.n_rules() > 0) )
228 {
229 std::string::size_type beg(1);
230 if (token.length()==1)
231 {
232 s >> token;
233 beg = 0;
234 }
235
236 std::string::size_type
237 end = token.find_last_of(')');
238 std::string::size_type
239 size = end==std::string::npos ? token.size()-beg : end-beg;
240
241 std::stringstream v;
242 MyMesh::Scalar coeff;
243 std::cout << " " << token << std::endl;
244 std::cout << " " << beg << " " << end << " " << size << std::endl;
245 v << token.substr(beg, size);
246 v >> coeff;
247 std::cout << " coeffecient " << coeff << std::endl;
248 subdivider.rule( subdivider.n_rules()-1 ).set_coeff(coeff);
249
250 if (end == std::string::npos)
251 {
252 s >> token;
253 if (token[0]!=')')
254 {
255 std::cerr << "Syntax error: Missing ')'\n";
256 return 1;
257 }
258 }
259 }
260 else
261 {
262 std::cerr << "Syntax error: " << token << "?\n";
263 return 1;
264 }
265 }
266 }
267
268 std::cout << "Rule sequence : "
269 << subdivider.rules_as_string() << std::endl;
270
271 // -------------------- Initialize subdivider
272 std::cout << "Initialize subdivider\n";
273 if (!subdivider.initialize())
274 {
275 std::cerr << " Error!\n";
276 return 1;
277 }
278
279 //
280 MyMesh::FaceFaceIter ff_it;
281 double quality(0.0);
282
283 // ---------------------------------------- subdivide
284 std::cout << "\nSubdividing...\n";
285
286 OpenMesh::Utils::Timer timer, timer2;
287 size_t i;
288
289 if ( uniform )
290 { // unifom
291 MyMesh::VertexHandle vh;
292 MyMesh::VertexIter v_it;
293 MyMesh::FaceHandle fh;
294 MyMesh::FaceIter f_it;
295
296
297 // raise all vertices to target state
298 timer.start();
299
300 size_t n = n_iter;
301 size_t n_rules = subdivider.n_rules();
302
303 i = 0;
304
305 // calculate target states for faces and vertices
306 size_t target1 = (n - 1) * n_rules + subdivider.subdiv_rule().number() + 1;
307 size_t target2 = n * n_rules;
308
309 for (f_it = mesh.faces_begin(); f_it != mesh.faces_end(); ++f_it) {
310
311 if (mesh.data(*f_it).state() < int(target1) ) {
312 ++i;
313 fh = *f_it;
314 timer2.start();
315 subdivider.refine(fh);
316 timer2.stop();
317 }
318 }
319
320 for (v_it = mesh.vertices_begin(); v_it != mesh.vertices_end(); ++v_it) {
321
322 if (mesh.data(*v_it).state() < int(target2) ) {
323 vh = *v_it;
324 timer2.cont();
325 subdivider.refine(vh);
326 timer2.stop();
327 }
328 }
329 timer.stop();
330 }
331 else
332 { // adaptive
333
334 MyMesh::FaceIter f_it;
335 MyMesh::FaceHandle fh;
336
337 std::vector<double> __acos;
338 size_t buckets(3000);
339 double range(2.0);
340 double range2bucket(buckets/range);
341
342 for (i = 0; i < buckets; ++i)
343 __acos.push_back( acos(-1.0 + i * range / buckets) );
344
345 timer.start(); // total time needed
346
347 // n iterations or until desired number of vertices reached approx.
348 for (i = 0; i < n_iter && mesh.n_vertices() < max_nv; ++i)
349 {
350 mesh.update_face_normals();
351
352 // calculate quality
353 quality = 0.0;
354
355 fh = *(mesh.faces_begin());
356
357 // check every face
358 for (f_it = mesh.faces_begin(); f_it != mesh.faces_end(); ++f_it) {
359
360 double face_quality = 0.0;
361 int valence = 0;
362
363 for (ff_it = mesh.ff_iter(*f_it); ff_it.is_valid(); ++ff_it) {
364
365 double temp_quality = OpenMesh::dot( mesh.normal(*f_it), mesh.normal(*ff_it) );
366
367 if (temp_quality >= 1.0)
368 temp_quality = .99;
369 else if (temp_quality <= -1.0)
370 temp_quality = -.99;
371 temp_quality = (1.0+temp_quality) * range2bucket;
372 face_quality += __acos[int(temp_quality+.5)];
373
374 ++valence;
375 }
376
377 face_quality /= valence;
378
379 // calaculate face area
380 MyMesh::Point p1, p2, p3;
381 MyMesh::Scalar area;
382
383 #define heh halfedge_handle
384 #define nheh next_halfedge_handle
385 #define tvh to_vertex_handle
386 #define fvh from_vertex_handle
387 p1 = mesh.point(mesh.tvh(mesh.heh(*f_it)));
388 p2 = mesh.point(mesh.fvh(mesh.heh(*f_it)));
389 p3 = mesh.point(mesh.tvh(mesh.nheh(mesh.heh(*f_it))));
390 #undef heh
391 #undef nheh
392 #undef tvh
393 #undef fvh
394
395 area = ((p2 - p1) % (p3 - p1)).norm();
396
397 // weight face_quality
398 face_quality *= pow(double(area), double(.1));
399 //face_quality *= area;
400
401 if (face_quality >= quality && !mesh.is_boundary(*f_it))
402 {
403 quality = face_quality;
404 fh = *f_it;
405 }
406 }
407
408 // Subdivide Face
409 timer2.cont();
410 subdivider.refine(fh);
411 timer2.stop();
412 }
413
414 // calculate time
415 timer.stop();
416
417 } // uniform/adaptive?
418
419 // calculate maximum refinement level
420 Adaptive::state_t max_level(0);
421
422 for (MyMesh::VertexIter v_it = mesh.vertices_begin();
423 v_it != mesh.vertices_end(); ++v_it)
424 {
425 if (mesh.data(*v_it).state() > max_level)
426 max_level = mesh.data(*v_it).state();
427 }
428
429
430 // output results
431 std::cout << "\nDid " << i << (uniform ? " uniform " : "" )
432 << " subdivision steps in "
433 << timer.as_string()
434 << ", " << i/timer.seconds() << " steps/s\n";
435 std::cout << " only refinement: " << timer2.as_string()
436 << ", " << i/timer2.seconds() << " steps/s\n\n";
437
438 std::cout << "Before: ";
439 std::cout << n_vertices << " Vertices, ";
440 std::cout << n_edges << " Edges, ";
441 std::cout << n_faces << " Faces. \n";
442
443 std::cout << "Now : ";
444 std::cout << mesh.n_vertices() << " Vertices, ";
445 std::cout << mesh.n_edges() << " Edges, ";
446 std::cout << mesh.n_faces() << " Faces. \n\n";
447
448 std::cout << "Maximum quality : " << quality << std::endl;
449 std::cout << "Maximum Subdivision Level: " << max_level/subdivider.n_rules()
450 << std::endl << std::endl;
451
452 // ---------------------------------------- write mesh to file
453 {
454 if ( ofname.empty() )
455 {
456 std::stringstream s;
457
458 s << "result." << subdivider.rules_as_string("_")
459 << "-" << i << "x.off";
460 s >> ofname;
461 }
462
463 std::cout << "Output file: '" << ofname << "'.\n";
464 if (!OpenMesh::IO::write_mesh(mesh, ofname, OpenMesh::IO::Options::Binary))
465 {
466 std::cerr << " Error writing file!\n";
467 return 1;
468 }
469 }
470 return 0;
471 }
472
473 // ----------------------------------------------------------------------------
474 // helper
475
476 void usage_and_exit(const std::string& _fname, int xcode)
477 {
478 using namespace std;
479
480 cout << endl
481 << "Usage: " << basename(_fname)
482 << " [Options] input-mesh [output-mesh]\n\n";
483 cout << "\tAdaptively refine an input-mesh. The refined mesh is stored in\n"
484 << "\ta file named \"result.XXX.off\" (binary .off), if not specified\n"
485 << "\texplicitely (optional 2nd parameter of command line).\n\n";
486 cout << "Options:\n\n";
487 cout << "-m <int>\n\tAdaptively refine up to approx. <int> vertices.\n\n"
488 << "-n <int>\n\tAdaptively refine <int> times.\n\n"
489 << "-r <rule sequence>\n\tDefine a custom rule sequence.\n\n"
490 << "-l\n\tUse rule sequence for adaptive Loop.\n\n"
491 << "-s\n\tUse rule sequence for adaptive sqrt(3).\n\n"
492 << "-U\n\tRefine mesh uniformly (simulates uniform subdivision).\n\n";
493
494 exit(xcode);
495 }
496
497 std::string basename(const std::string& _f)
498 {
499 std::string::size_type dot = _f.rfind("/");
500 if (dot == std::string::npos)
501 return _f;
502 return _f.substr(dot+1, _f.length()-(dot+1));
503 }
504
505 // ----------------------------------------------------------------------------
506 // end of file
507 // ============================================================================
OpenMesh::PolyMeshT::FaceFaceIter
Kernel::FaceFaceIter FaceFaceIter
Circulator.
Definition: PolyMeshT.hh:170
RulesT.hh
OpenMesh::IO::write_mesh
bool write_mesh(const Mesh &_mesh, const std::string &_filename, Options _opt=Options::Default, std::streamsize _precision=6)
Write a mesh to the file _filename.
Definition: MeshIO.hh:207
OpenMesh::PolyMeshT::Point
Kernel::Point Point
Coordinate type.
Definition: PolyMeshT.hh:112
OpenMesh::TriMesh_ArrayKernelT
Definition: PolyMesh_ArrayKernelT.hh:73
std
STL namespace.
CompositeT.hh
OpenMesh::Utils::Timer::as_string
std::string as_string(Format format=Automatic)
OpenMesh::Utils::Timer::seconds
double seconds(void) const
Returns measured time in seconds, if the timer is in state &#39;Stopped&#39;.
OpenMesh::Utils::Timer::stop
void stop(void)
Stop measurement.
OpenMesh::Utils::Timer::cont
void cont(void)
Continue measurement.
OpenMesh::IO::read_mesh
bool read_mesh(Mesh &_mesh, const std::string &_filename)
Read a mesh from file _filename.
Definition: MeshIO.hh:112
Timer.hh
OpenMesh::Subdivider::Adaptive::CompositeTraits
Definition: CompositeTraits.hh:81
OpenMesh::PolyMeshT::Scalar
Kernel::Scalar Scalar
Scalar type.
Definition: PolyMeshT.hh:110
OpenMesh::PolyMeshT::update_face_normals
void update_face_normals()
Update normal vectors for all faces.
Definition: PolyMeshT_impl.hh:335
OpenMesh::Subdivider
OpenMesh::Utils::Timer::start
void start(void)
Start measurement.
OpenMesh::Subdivider::Adaptive::state_t
CompositeTraits::state_t state_t
Definition: CompositeTraits.hh:250
OpenMesh::Utils::Timer
Definition: Timer.hh:82
RuleMap
Definition: adaptive_subdivider.cc:107
OpenMesh::Subdivider::Adaptive::RuleInterfaceT
Definition: CompositeT.hh:79
OpenMesh::IO::Options::Binary
Set binary mode for r/w.
Definition: Options.hh:100
OpenMesh::dot
osg::Vec3f::ValueType dot(const osg::Vec3f &_v1, const osg::Vec3f &_v2)
Adapter for osg vector member computing a scalar product.
Definition: VectorAdapter.hh:176
OpenMesh::Subdivider::Adaptive::CompositeT
Definition: CompositeT.hh:133
OpenMesh::PolyMeshT::VertexHandle
Kernel::VertexHandle VertexHandle
Handle for referencing the corresponding item.
Definition: PolyMeshT.hh:136