| 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 explicit 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 override { 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 ) override
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 ) override
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) override
178 {
179
181
182 // Do _n subdivisions
183 for (size_t i=0; i < _n; ++i)
184 {
185
186 // This is an interpolating scheme, old vertices remain the same.
187 typename mesh_t::VertexIter initialVerticesEnd = _m.vertices_end();
188 for ( auto vh : _m.vertices())
189 _m.property( vp_pos_, vh ) = _m.point(vh);
190
191 // Compute position for new vertices and store them in the edge property
192 for (auto eh : _m.edges())
193 compute_midpoint( _m, eh);
194
195
196 // Split each edge at midpoint and store precomputed positions (stored in
197 // edge property ep_pos_) in the vertex property vp_pos_;
198
199 // Attention! Creating new edges, hence make sure the loop ends correctly.
200 for (auto eh : _m.edges())
201 split_edge(_m, eh );
202
203
204 // Commit changes in topology and reconsitute consistency
205
206 // Attention! Creating new faces, hence make sure the loop ends correctly.
207 for (auto fh : _m.faces())
208 split_face(_m, fh );
209
210
211 // Commit changes in geometry
212 for ( auto vh : _m.vertices())
213 _m.set_point(vh, _m.property( vp_pos_, vh ) );
214
215 #if defined(_DEBUG) || defined(DEBUG)
216 // Now we have an consistent mesh!
217 assert( OpenMesh::Utils::MeshCheckerT<mesh_t>(_m).check() );
218 #endif
219 }
220
221 return true;
222 }
223
224 private: // topological modifiers
225
226 void split_face(mesh_t& _m, const typename mesh_t::FaceHandle& _fh)
227 {
228 typename mesh_t::HalfedgeHandle
229 heh1(_m.halfedge_handle(_fh)),
230 heh2(_m.next_halfedge_handle(_m.next_halfedge_handle(heh1))),
231 heh3(_m.next_halfedge_handle(_m.next_halfedge_handle(heh2)));
232
233 // Cutting off every corner of the 6_gon
234 corner_cutting( _m, heh1 );
235 corner_cutting( _m, heh2 );
236 corner_cutting( _m, heh3 );
237 }
238
239
240 void corner_cutting(mesh_t& _m, const typename mesh_t::HalfedgeHandle& _he)
241 {
242 // Define Halfedge Handles
243 typename mesh_t::HalfedgeHandle
244 heh1(_he),
245 heh5(heh1),
246 heh6(_m.next_halfedge_handle(heh1));
247
248 // Cycle around the polygon to find correct Halfedge
249 for (; _m.next_halfedge_handle(_m.next_halfedge_handle(heh5)) != heh1;
250 heh5 = _m.next_halfedge_handle(heh5))
251 {}
252
253 typename mesh_t::VertexHandle
254 vh1 = _m.to_vertex_handle(heh1),
255 vh2 = _m.to_vertex_handle(heh5);
256
257 typename mesh_t::HalfedgeHandle
258 heh2(_m.next_halfedge_handle(heh5)),
259 heh3(_m.new_edge( vh1, vh2)),
260 heh4(_m.opposite_halfedge_handle(heh3));
261
262 /* Intermediate result
263 *
264 * *
265 * 5 /|\
266 * /_ \
267 * vh2> * *
268 * /|\3 |\
269 * /_ \|4 \
270 * *----\*----\*
271 * 1 ^ 6
272 * vh1 (adjust_outgoing halfedge!)
273 */
274
275 // Old and new Face
276 typename mesh_t::FaceHandle fh_old(_m.face_handle(heh6));
277 typename mesh_t::FaceHandle fh_new(_m.new_face());
278
279
280 // Re-Set Handles around old Face
281 _m.set_next_halfedge_handle(heh4, heh6);
282 _m.set_next_halfedge_handle(heh5, heh4);
283
284 _m.set_face_handle(heh4, fh_old);
285 _m.set_face_handle(heh5, fh_old);
286 _m.set_face_handle(heh6, fh_old);
287 _m.set_halfedge_handle(fh_old, heh4);
288
289 // Re-Set Handles around new Face
290 _m.set_next_halfedge_handle(heh1, heh3);
291 _m.set_next_halfedge_handle(heh3, heh2);
292
293 _m.set_face_handle(heh1, fh_new);
294 _m.set_face_handle(heh2, fh_new);
295 _m.set_face_handle(heh3, fh_new);
296
297 _m.set_halfedge_handle(fh_new, heh1);
298 }
299
300
301 void split_edge(mesh_t& _m, const typename mesh_t::EdgeHandle& _eh)
302 {
303 typename mesh_t::HalfedgeHandle
304 heh = _m.halfedge_handle(_eh, 0),
305 opp_heh = _m.halfedge_handle(_eh, 1);
306
307 typename mesh_t::HalfedgeHandle new_heh, opp_new_heh, t_heh;
308 typename mesh_t::VertexHandle vh;
309 typename mesh_t::VertexHandle vh1(_m.to_vertex_handle(heh));
310 typename mesh_t::Point zero(0,0,0);
311
312 // new vertex
313 vh = _m.new_vertex( zero );
314
315 // memorize position, will be set later
316 _m.property( vp_pos_, vh ) = _m.property( ep_pos_, _eh );
317
318
319 // Re-link mesh entities
320 if (_m.is_boundary(_eh))
321 {
322 for (t_heh = heh;
323 _m.next_halfedge_handle(t_heh) != opp_heh;
324 t_heh = _m.opposite_halfedge_handle(_m.next_halfedge_handle(t_heh)))
325 {}
326 }
327 else
328 {
329 for (t_heh = _m.next_halfedge_handle(opp_heh);
330 _m.next_halfedge_handle(t_heh) != opp_heh;
331 t_heh = _m.next_halfedge_handle(t_heh) )
332 {}
333 }
334
335 new_heh = _m.new_edge(vh, vh1);
336 opp_new_heh = _m.opposite_halfedge_handle(new_heh);
337 _m.set_vertex_handle( heh, vh );
338
339 _m.set_next_halfedge_handle(t_heh, opp_new_heh);
340 _m.set_next_halfedge_handle(new_heh, _m.next_halfedge_handle(heh));
341 _m.set_next_halfedge_handle(heh, new_heh);
342 _m.set_next_halfedge_handle(opp_new_heh, opp_heh);
343
344 if (_m.face_handle(opp_heh).is_valid())
345 {
346 _m.set_face_handle(opp_new_heh, _m.face_handle(opp_heh));
347 _m.set_halfedge_handle(_m.face_handle(opp_new_heh), opp_new_heh);
348 }
349
350 _m.set_face_handle( new_heh, _m.face_handle(heh) );
351 _m.set_halfedge_handle( vh, new_heh);
352
353 // We cant reconnect a non existing face, so we skip this here if necessary
354 if ( !_m.is_boundary(heh) )
355 _m.set_halfedge_handle( _m.face_handle(heh), heh );
356
357 _m.set_halfedge_handle( vh1, opp_new_heh );
358
359 // Never forget this, when playing with the topology
360 _m.adjust_outgoing_halfedge( vh );
361 _m.adjust_outgoing_halfedge( vh1 );
362 }
363
364 private: // geometry helper
365
366 void compute_midpoint(mesh_t& _m, const typename mesh_t::EdgeHandle& _eh)
367 {
368 typename mesh_t::HalfedgeHandle heh, opp_heh;
369
370 heh = _m.halfedge_handle( _eh, 0);
371 opp_heh = _m.halfedge_handle( _eh, 1);
372
373 typename mesh_t::Point pos(0,0,0);
374
375 typename mesh_t::VertexHandle a_0(_m.to_vertex_handle(heh));
376 typename mesh_t::VertexHandle a_1(_m.to_vertex_handle(opp_heh));
377
378 // boundary edge: 4-point scheme
379 if (_m.is_boundary(_eh) )
380 {
381 pos = _m.point(a_0);
382 pos += _m.point(a_1);
383 pos *= static_cast<RealType>(9.0/16.0);
384 typename mesh_t::Point tpos;
385 if(_m.is_boundary(heh))
386 {
387 tpos = _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
388 tpos += _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh))));
389 }
390 else
391 {
392 assert(_m.is_boundary(opp_heh));
393 tpos = _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(opp_heh)));
394 tpos += _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(opp_heh))));
395 }
396 tpos *= static_cast<RealType>(-1.0/16.0);
397 pos += tpos;
398 }
399 else
400 {
401 int valence_a_0 = _m.valence(a_0);
402 int valence_a_1 = _m.valence(a_1);
403 assert(valence_a_0>2);
404 assert(valence_a_1>2);
405
406 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
407 {
408 real_t alpha = real_t(1.0/2);
409 real_t beta = real_t(1.0/8);
410 real_t gamma = real_t(-1.0/16);
411
412 //get points
413 typename mesh_t::VertexHandle b_0, b_1, c_0, c_1, c_2, c_3;
414 typename mesh_t::HalfedgeHandle t_he;
415
416 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(heh));
417 b_0 = _m.to_vertex_handle(t_he);
418 if(!_m.is_boundary(_m.opposite_halfedge_handle(t_he)))
419 {
420 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
421 c_0 = _m.to_vertex_handle(t_he);
422 }
423
424 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh));
425 b_1 = _m.to_vertex_handle(t_he);
426 if(!_m.is_boundary(t_he))
427 {
428 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
429 c_1 = _m.to_vertex_handle(t_he);
430 }
431
432 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(opp_heh));
433 assert(b_1.idx()==_m.to_vertex_handle(t_he).idx());
434 if(!_m.is_boundary(_m.opposite_halfedge_handle(t_he)))
435 {
436 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
437 c_2 = _m.to_vertex_handle(t_he);
438 }
439
440 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(opp_heh));
441 assert(b_0==_m.to_vertex_handle(t_he));
442 if(!_m.is_boundary(t_he))
443 {
444 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
445 c_3 = _m.to_vertex_handle(t_he);
446 }
447
448 //compute position.
449 //a0,a1,b0,b1 must exist.
450 assert(a_0.is_valid());
451 assert(a_1.is_valid());
452 assert(b_0.is_valid());
453 assert(b_1.is_valid());
454 //The other vertices may be created from symmetry is they are on the other side of the boundary.
455
456 pos = _m.point(a_0);
457 pos += _m.point(a_1);
458 pos *= alpha;
459
460 typename mesh_t::Point tpos ( _m.point(b_0) );
461 tpos += _m.point(b_1);
462 tpos *= beta;
463 pos += tpos;
464
465 typename mesh_t::Point pc_0, pc_1, pc_2, pc_3;
466 if(c_0.is_valid())
467 pc_0 = _m.point(c_0);
468 else //create the point by symmetry
469 {
470 pc_0 = _m.point(a_1) + _m.point(b_0) - _m.point(a_0);
471 }
472 if(c_1.is_valid())
473 pc_1 = _m.point(c_1);
474 else //create the point by symmetry
475 {
476 pc_1 = _m.point(a_1) + _m.point(b_1) - _m.point(a_0);
477 }
478 if(c_2.is_valid())
479 pc_2 = _m.point(c_2);
480 else //create the point by symmetry
481 {
482 pc_2 = _m.point(a_0) + _m.point(b_1) - _m.point(a_1);
483 }
484 if(c_3.is_valid())
485 pc_3 = _m.point(c_3);
486 else //create the point by symmetry
487 {
488 pc_3 = _m.point(a_0) + _m.point(b_0) - _m.point(a_1);
489 }
490 tpos = pc_0;
491 tpos += pc_1;
492 tpos += pc_2;
493 tpos += pc_3;
494 tpos *= gamma;
495 pos += tpos;
496 }
497 else //at least one endpoint is [irregular and not in boundary]
498 {
499 RealType normFactor = static_cast<RealType>(0.0);
500
501 if(valence_a_0!=6 && !_m.is_boundary(a_0))
502 {
503 assert((int)weights[valence_a_0].size()==valence_a_0+1);
504 typename mesh_t::HalfedgeHandle t_he = opp_heh;
505 for(int i = 0; i < valence_a_0 ; t_he=_m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he)), ++i)
506 {
507 pos += weights[valence_a_0][i] * _m.point(_m.to_vertex_handle(t_he));
508 }
509 assert(t_he==opp_heh);
510
511 //add irregular vertex:
512 pos += weights[valence_a_0][valence_a_0] * _m.point(a_0);
513 ++normFactor;
514 }
515
516 if(valence_a_1!=6 && !_m.is_boundary(a_1))
517 {
518 assert((int)weights[valence_a_1].size()==valence_a_1+1);
519 typename mesh_t::HalfedgeHandle t_he = heh;
520 for(int i = 0; i < valence_a_1 ; t_he=_m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he)), ++i)
521 {
522 pos += weights[valence_a_1][i] * _m.point(_m.to_vertex_handle(t_he));
523 }
524 assert(t_he==heh);
525 //add irregular vertex:
526 pos += weights[valence_a_1][valence_a_1] * _m.point(a_1);
527 ++normFactor;
528 }
529
530 assert(normFactor>0.1); //normFactor should be 1 or 2
531
532 //if both vertices are irregular, average positions:
533 pos /= normFactor;
534 }
535 }
536 _m.property( ep_pos_, _eh ) = pos;
537 }
538
539 private: // data
540
541 OpenMesh::VPropHandleT< typename mesh_t::Point > vp_pos_;
542 OpenMesh::EPropHandleT< typename mesh_t::Point > ep_pos_;
543
544 weights_t weights;
545
546 };
547
548 } // END_NS_UNIFORM
549 } // END_NS_SUBDIVIDER
550 } // END_NS_OPENMESH
551 #endif
552
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::subdivide
bool subdivide(MeshType &_m, size_t _n, const bool _update_points=true) override
Subdivide mesh _m _n times.
Definition: ModifiedButterFlyT.hh:177
OpenMesh::VPropHandleT< typename mesh_t::Point >
OpenMesh::Subdivider::Uniform::ModifiedButterflyT
Modified Butterfly subdivision algorithm.
Definition: ModifiedButterFlyT.hh:92
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::cleanup
bool cleanup(mesh_t &_m) override
Cleanup mesh after usage, e.g. remove added properties.
Definition: ModifiedButterFlyT.hh:169
OpenMesh::EPropHandleT< typename mesh_t::Point >
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::name
const char * name() const override
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::Subdivider::Uniform::ModifiedButterflyT::prepare
bool prepare(mesh_t &_m) override
Prepare mesh, e.g. add properties.
Definition: ModifiedButterFlyT.hh:161
OpenMesh::Subdivider::Uniform::ModifiedButterflyT::init_weights
void init_weights(size_t _max_valence=30)
Pre-compute weights.
Definition: ModifiedButterFlyT.hh:124
OpenMesh::Utils::MeshCheckerT
Check integrity of mesh.
Definition: MeshCheckerT.hh:73
SubdividerT.hh
OpenMesh
Contains all the mesh ingredients like the polygonal mesh, the triangle mesh, different mesh kernels ...
Definition: MeshItems.hh:59
OpenMesh::PolyMeshT::VertexHandle
Kernel::VertexHandle VertexHandle
Handle for referencing the corresponding item.
Definition: PolyMeshT.hh:136