| OpenMesh

ModifiedButterFlyT.hh
Go to the documentation of this file.
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
51 //=============================================================================
52 //
53 // CLASS ModifiedButterflyT
54 //
55 //=============================================================================
56
57
58 #ifndef SP_MODIFIED_BUTTERFLY_H
59 #define SP_MODIFIED_BUTTERFLY_H
60
61 #include <OpenMesh/Tools/Subdivider/Uniform/SubdividerT.hh>
62 #include <OpenMesh/Core/Utils/vector_cast.hh>
63 #include <OpenMesh/Core/Utils/Property.hh>
64 // -------------------- STL
65 #include <vector>
66 #if defined(OM_CC_MIPS)
67 # include <math.h>
68 #else
69 # include <cmath>
70 #endif
71
72
73 //== NAMESPACE ================================================================
74
75 namespace OpenMesh { // BEGIN_NS_OPENMESH
76 namespace Subdivider { // BEGIN_NS_DECIMATER
77 namespace Uniform { // BEGIN_NS_UNIFORM
78
79
80 //== CLASS DEFINITION =========================================================
81
82
91 template <typename MeshType, typename RealType = double>
92 class ModifiedButterflyT : public SubdividerT<MeshType, RealType>
93 {
94 public:
95
96 typedef RealType real_t;
97 typedef MeshType mesh_t;
98 typedef SubdividerT< mesh_t, real_t > parent_t;
99
100 typedef std::vector< std::vector<real_t> > weights_t;
101 typedef std::vector<real_t> weight_t;
102
103 public:
104
105
106 ModifiedButterflyT() : parent_t()
107 { init_weights(); }
108
109
110 ModifiedButterflyT( mesh_t& _m) : parent_t(_m)
111 { init_weights(); }
112
113
114 ~ModifiedButterflyT() {}
115
116
117 public:
118
119
120 const char *name() const { return "Uniform Spectral"; }
121
122
124 void init_weights(size_t _max_valence=30)
125 {
126 weights.resize(_max_valence);
127
128 //special case: K==3, K==4
129 weights[3].resize(4);
130 weights[3][0] = real_t(5.0)/12;
131 weights[3][1] = real_t(-1.0)/12;
132 weights[3][2] = real_t(-1.0)/12;
133 weights[3][3] = real_t(3.0)/4;
134
135 weights[4].resize(5);
136 weights[4][0] = real_t(3.0)/8;
137 weights[4][1] = 0;
138 weights[4][2] = real_t(-1.0)/8;
139 weights[4][3] = 0;
140 weights[4][4] = real_t(3.0)/4;
141
142 for(unsigned int K = 5; K<_max_valence; ++K)
143 {
144 weights[K].resize(K+1);
145 // s(j) = ( 1/4 + cos(2*pi*j/K) + 1/2 * cos(4*pi*j/K) )/K
146 double invK = 1.0/static_cast<double>(K);
147 real_t sum = 0;
148 for(unsigned int j=0; j<K; ++j)
149 {
150 weights[K][j] = static_cast<real_t>((0.25 + cos(2.0*M_PI*static_cast<double>(j)*invK) + 0.5*cos(4.0*M_PI*static_cast<double>(j)*invK))*invK);
151 sum += weights[K][j];
152 }
153 weights[K][K] = static_cast<real_t>(1.0) - sum;
154 }
155 }
156
157
158 protected:
159
160
161 bool prepare( mesh_t& _m )
162 {
163 _m.add_property( vp_pos_ );
164 _m.add_property( ep_pos_ );
165 return true;
166 }
167
168
169 bool cleanup( mesh_t& _m )
170 {
171 _m.remove_property( vp_pos_ );
172 _m.remove_property( ep_pos_ );
173 return true;
174 }
175
176
177 bool subdivide( MeshType& _m, size_t _n , const bool _update_points = true)
178 {
179
181
182 typename mesh_t::FaceIter fit, f_end;
183 typename mesh_t::EdgeIter eit, e_end;
184 typename mesh_t::VertexIter vit;
185
186 // Do _n subdivisions
187 for (size_t i=0; i < _n; ++i)
188 {
189
190 // This is an interpolating scheme, old vertices remain the same.
191 typename mesh_t::VertexIter initialVerticesEnd = _m.vertices_end();
192 for ( vit = _m.vertices_begin(); vit != initialVerticesEnd; ++vit)
193 _m.property( vp_pos_, *vit ) = _m.point(*vit);
194
195 // Compute position for new vertices and store them in the edge property
196 for (eit=_m.edges_begin(); eit != _m.edges_end(); ++eit)
197 compute_midpoint( _m, *eit );
198
199
200 // Split each edge at midpoint and store precomputed positions (stored in
201 // edge property ep_pos_) in the vertex property vp_pos_;
202
203 // Attention! Creating new edges, hence make sure the loop ends correctly.
204 e_end = _m.edges_end();
205 for (eit=_m.edges_begin(); eit != e_end; ++eit)
206 split_edge(_m, *eit );
207
208
209 // Commit changes in topology and reconsitute consistency
210
211 // Attention! Creating new faces, hence make sure the loop ends correctly.
212 f_end = _m.faces_end();
213 for (fit = _m.faces_begin(); fit != f_end; ++fit)
214 split_face(_m, *fit );
215
216
217 // Commit changes in geometry
218 for ( vit = /*initialVerticesEnd;*/_m.vertices_begin();
219 vit != _m.vertices_end(); ++vit)
220 _m.set_point(*vit, _m.property( vp_pos_, *vit ) );
221
222 #if defined(_DEBUG) || defined(DEBUG)
223 // Now we have an consistent mesh!
224 assert( OpenMesh::Utils::MeshCheckerT<mesh_t>(_m).check() );
225 #endif
226 }
227
228 return true;
229 }
230
231 private: // topological modifiers
232
233 void split_face(mesh_t& _m, const typename mesh_t::FaceHandle& _fh)
234 {
235 typename mesh_t::HalfedgeHandle
236 heh1(_m.halfedge_handle(_fh)),
237 heh2(_m.next_halfedge_handle(_m.next_halfedge_handle(heh1))),
238 heh3(_m.next_halfedge_handle(_m.next_halfedge_handle(heh2)));
239
240 // Cutting off every corner of the 6_gon
241 corner_cutting( _m, heh1 );
242 corner_cutting( _m, heh2 );
243 corner_cutting( _m, heh3 );
244 }
245
246
247 void corner_cutting(mesh_t& _m, const typename mesh_t::HalfedgeHandle& _he)
248 {
249 // Define Halfedge Handles
250 typename mesh_t::HalfedgeHandle
251 heh1(_he),
252 heh5(heh1),
253 heh6(_m.next_halfedge_handle(heh1));
254
255 // Cycle around the polygon to find correct Halfedge
256 for (; _m.next_halfedge_handle(_m.next_halfedge_handle(heh5)) != heh1;
257 heh5 = _m.next_halfedge_handle(heh5))
258 {}
259
260 typename mesh_t::VertexHandle
261 vh1 = _m.to_vertex_handle(heh1),
262 vh2 = _m.to_vertex_handle(heh5);
263
264 typename mesh_t::HalfedgeHandle
265 heh2(_m.next_halfedge_handle(heh5)),
266 heh3(_m.new_edge( vh1, vh2)),
267 heh4(_m.opposite_halfedge_handle(heh3));
268
269 /* Intermediate result
270 *
271 * *
272 * 5 /|\
273 * /_ \
274 * vh2> * *
275 * /|\3 |\
276 * /_ \|4 \
277 * *----\*----\*
278 * 1 ^ 6
279 * vh1 (adjust_outgoing halfedge!)
280 */
281
282 // Old and new Face
283 typename mesh_t::FaceHandle fh_old(_m.face_handle(heh6));
284 typename mesh_t::FaceHandle fh_new(_m.new_face());
285
286
287 // Re-Set Handles around old Face
288 _m.set_next_halfedge_handle(heh4, heh6);
289 _m.set_next_halfedge_handle(heh5, heh4);
290
291 _m.set_face_handle(heh4, fh_old);
292 _m.set_face_handle(heh5, fh_old);
293 _m.set_face_handle(heh6, fh_old);
294 _m.set_halfedge_handle(fh_old, heh4);
295
296 // Re-Set Handles around new Face
297 _m.set_next_halfedge_handle(heh1, heh3);
298 _m.set_next_halfedge_handle(heh3, heh2);
299
300 _m.set_face_handle(heh1, fh_new);
301 _m.set_face_handle(heh2, fh_new);
302 _m.set_face_handle(heh3, fh_new);
303
304 _m.set_halfedge_handle(fh_new, heh1);
305 }
306
307
308 void split_edge(mesh_t& _m, const typename mesh_t::EdgeHandle& _eh)
309 {
310 typename mesh_t::HalfedgeHandle
311 heh = _m.halfedge_handle(_eh, 0),
312 opp_heh = _m.halfedge_handle(_eh, 1);
313
314 typename mesh_t::HalfedgeHandle new_heh, opp_new_heh, t_heh;
315 typename mesh_t::VertexHandle vh;
316 typename mesh_t::VertexHandle vh1(_m.to_vertex_handle(heh));
317 typename mesh_t::Point zero(0,0,0);
318
319 // new vertex
320 vh = _m.new_vertex( zero );
321
322 // memorize position, will be set later
323 _m.property( vp_pos_, vh ) = _m.property( ep_pos_, _eh );
324
325
326 // Re-link mesh entities
327 if (_m.is_boundary(_eh))
328 {
329 for (t_heh = heh;
330 _m.next_halfedge_handle(t_heh) != opp_heh;
331 t_heh = _m.opposite_halfedge_handle(_m.next_halfedge_handle(t_heh)))
332 {}
333 }
334 else
335 {
336 for (t_heh = _m.next_halfedge_handle(opp_heh);
337 _m.next_halfedge_handle(t_heh) != opp_heh;
338 t_heh = _m.next_halfedge_handle(t_heh) )
339 {}
340 }
341
342 new_heh = _m.new_edge(vh, vh1);
343 opp_new_heh = _m.opposite_halfedge_handle(new_heh);
344 _m.set_vertex_handle( heh, vh );
345
346 _m.set_next_halfedge_handle(t_heh, opp_new_heh);
347 _m.set_next_halfedge_handle(new_heh, _m.next_halfedge_handle(heh));
348 _m.set_next_halfedge_handle(heh, new_heh);
349 _m.set_next_halfedge_handle(opp_new_heh, opp_heh);
350
351 if (_m.face_handle(opp_heh).is_valid())
352 {
353 _m.set_face_handle(opp_new_heh, _m.face_handle(opp_heh));
354 _m.set_halfedge_handle(_m.face_handle(opp_new_heh), opp_new_heh);
355 }
356
357 _m.set_face_handle( new_heh, _m.face_handle(heh) );
358 _m.set_halfedge_handle( vh, new_heh);
359 _m.set_halfedge_handle( _m.face_handle(heh), heh );
360 _m.set_halfedge_handle( vh1, opp_new_heh );
361
362 // Never forget this, when playing with the topology
363 _m.adjust_outgoing_halfedge( vh );
364 _m.adjust_outgoing_halfedge( vh1 );
365 }
366
367 private: // geometry helper
368
369 void compute_midpoint(mesh_t& _m, const typename mesh_t::EdgeHandle& _eh)
370 {
371 typename mesh_t::HalfedgeHandle heh, opp_heh;
372
373 heh = _m.halfedge_handle( _eh, 0);
374 opp_heh = _m.halfedge_handle( _eh, 1);
375
376 typename mesh_t::Point pos(0,0,0);
377
378 typename mesh_t::VertexHandle a_0(_m.to_vertex_handle(heh));
379 typename mesh_t::VertexHandle a_1(_m.to_vertex_handle(opp_heh));
380
381 // boundary edge: 4-point scheme
382 if (_m.is_boundary(_eh) )
383 {
384 pos = _m.point(a_0);
385 pos += _m.point(a_1);
386 pos *= static_cast<typename mesh_t::Point::value_type>(9.0/16.0);
387 typename mesh_t::Point tpos;
388 if(_m.is_boundary(heh))
389 {
390 tpos = _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
391 tpos += _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh))));
392 }
393 else
394 {
395 assert(_m.is_boundary(opp_heh));
396 tpos = _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(opp_heh)));
397 tpos += _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(opp_heh))));
398 }
399 tpos *= static_cast<typename mesh_t::Point::value_type>(-1.0/16.0);
400 pos += tpos;
401 }
402 else
403 {
404 int valence_a_0 = _m.valence(a_0);
405 int valence_a_1 = _m.valence(a_1);
406 assert(valence_a_0>2);
407 assert(valence_a_1>2);
408
409 if( (valence_a_0==6 && valence_a_1==6) || (_m.is_boundary(a_0) && valence_a_1==6) || (_m.is_boundary(a_1) && valence_a_0==6) || (_m.is_boundary(a_0) && _m.is_boundary(a_1)) )// use 8-point scheme
410 {
411 real_t alpha = real_t(1.0/2);
412 real_t beta = real_t(1.0/8);
413 real_t gamma = real_t(-1.0/16);
414
415 //get points
416 typename mesh_t::VertexHandle b_0, b_1, c_0, c_1, c_2, c_3;
417 typename mesh_t::HalfedgeHandle t_he;
418
419 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(heh));
420 b_0 = _m.to_vertex_handle(t_he);
421 if(!_m.is_boundary(_m.opposite_halfedge_handle(t_he)))
422 {
423 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
424 c_0 = _m.to_vertex_handle(t_he);
425 }
426
427 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh));
428 b_1 = _m.to_vertex_handle(t_he);
429 if(!_m.is_boundary(t_he))
430 {
431 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
432 c_1 = _m.to_vertex_handle(t_he);
433 }
434
435 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(opp_heh));
436 assert(b_1.idx()==_m.to_vertex_handle(t_he).idx());
437 if(!_m.is_boundary(_m.opposite_halfedge_handle(t_he)))
438 {
439 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
440 c_2 = _m.to_vertex_handle(t_he);
441 }
442
443 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(opp_heh));
444 assert(b_0==_m.to_vertex_handle(t_he));
445 if(!_m.is_boundary(t_he))
446 {
447 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
448 c_3 = _m.to_vertex_handle(t_he);
449 }
450
451 //compute position.
452 //a0,a1,b0,b1 must exist.
453 assert(a_0.is_valid());
454 assert(a_1.is_valid());
455 assert(b_0.is_valid());
456 assert(b_1.is_valid());
457 //The other vertices may be created from symmetry is they are on the other side of the boundary.
458
459 pos = _m.point(a_0);
460 pos += _m.point(a_1);
461 pos *= alpha;
462
463 typename mesh_t::Point tpos ( _m.point(b_0) );
464 tpos += _m.point(b_1);
465 tpos *= beta;
466 pos += tpos;
467
468 typename mesh_t::Point pc_0, pc_1, pc_2, pc_3;
469 if(c_0.is_valid())
470 pc_0 = _m.point(c_0);
471 else //create the point by symmetry
472 {
473 pc_0 = _m.point(a_1) + _m.point(b_0) - _m.point(a_0);
474 }
475 if(c_1.is_valid())
476 pc_1 = _m.point(c_1);
477 else //create the point by symmetry
478 {
479 pc_1 = _m.point(a_1) + _m.point(b_1) - _m.point(a_0);
480 }
481 if(c_2.is_valid())
482 pc_2 = _m.point(c_2);
483 else //create the point by symmetry
484 {
485 pc_2 = _m.point(a_0) + _m.point(b_1) - _m.point(a_1);
486 }
487 if(c_3.is_valid())
488 pc_3 = _m.point(c_3);
489 else //create the point by symmetry
490 {
491 pc_3 = _m.point(a_0) + _m.point(b_0) - _m.point(a_1);
492 }
493 tpos = pc_0;
494 tpos += pc_1;
495 tpos += pc_2;
496 tpos += pc_3;
497 tpos *= gamma;
498 pos += tpos;
499 }
500 else //at least one endpoint is [irregular and not in boundary]
501 {
502 typename mesh_t::Point::value_type normFactor = static_cast<typename mesh_t::Point::value_type>(0.0);
503
504 if(valence_a_0!=6 && !_m.is_boundary(a_0))
505 {
506 assert((int)weights[valence_a_0].size()==valence_a_0+1);
507 typename mesh_t::HalfedgeHandle t_he = opp_heh;
508 for(int i = 0; i < valence_a_0 ; t_he=_m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he)), ++i)
509 {
510 pos += weights[valence_a_0][i] * _m.point(_m.to_vertex_handle(t_he));
511 }
512 assert(t_he==opp_heh);
513
514 //add irregular vertex:
515 pos += weights[valence_a_0][valence_a_0] * _m.point(a_0);
516 ++normFactor;
517 }
518
519 if(valence_a_1!=6 && !_m.is_boundary(a_1))
520 {
521 assert((int)weights[valence_a_1].size()==valence_a_1+1);
522 typename mesh_t::HalfedgeHandle t_he = heh;
523 for(int i = 0; i < valence_a_1 ; t_he=_m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he)), ++i)
524 {
525 pos += weights[valence_a_1][i] * _m.point(_m.to_vertex_handle(t_he));
526 }
527 assert(t_he==heh);
528 //add irregular vertex:
529 pos += weights[valence_a_1][valence_a_1] * _m.point(a_1);
530 ++normFactor;
531 }
532
533 assert(normFactor>0.1); //normFactor should be 1 or 2
534
535 //if both vertices are irregular, average positions:
536 pos /= normFactor;
537 }
538 }
539 _m.property( ep_pos_, _eh ) = pos;
540 }
541
542 private: // data
543
544 OpenMesh::VPropHandleT< typename mesh_t::Point > vp_pos_;
545 OpenMesh::EPropHandleT< typename mesh_t::Point > ep_pos_;
546
547 weights_t weights;
548
549 };
550
551 } // END_NS_UNIFORM
552 } // END_NS_SUBDIVIDER
553 } // END_NS_OPENMESH
554 #endif
555
OpenMesh::VPropHandleT< typename mesh_t::Point >
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::cleanup
bool cleanup(mesh_t &_m)
Cleanup mesh after usage, e.g. remove added properties.
Definition: ModifiedButterFlyT.hh:169
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::subdivide
bool subdivide(MeshType &_m, size_t _n, const bool _update_points=true)
Subdivide mesh _m _n times.
Definition: ModifiedButterFlyT.hh:177
SubdividerT.hh
OpenMesh::PolyMeshT::VertexHandle
Kernel::VertexHandle VertexHandle
Handle for referencing the corresponding item.
Definition: PolyMeshT.hh:136
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::prepare
bool prepare(mesh_t &_m)
Prepare mesh, e.g. add properties.
Definition: ModifiedButterFlyT.hh:161
OpenMesh
Contains all the mesh ingredients like the polygonal mesh, the triangle mesh, different mesh kernels ...
Definition: MeshItems.hh:59
OpenMesh::EPropHandleT< typename mesh_t::Point >
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::name
const char * name() const
Return name of subdivision algorithm.
Definition: ModifiedButterFlyT.hh:120
OpenMesh::Subdivider::Uniform::SubdividerT
Abstract base class for uniform subdivision algorithms.
Definition: SubdividerT.hh:88
OpenMesh::Utils::MeshCheckerT
Check integrity of mesh.
Definition: MeshCheckerT.hh:73
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::init_weights
void init_weights(size_t _max_valence=30)
Pre-compute weights.
Definition: ModifiedButterFlyT.hh:124
OpenMesh::Subdivider::Uniform::ModifiedButterflyT
Modified Butterfly subdivision algorithm.
Definition: ModifiedButterFlyT.hh:92