Developer Documentation

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 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
360 // We cant reconnect a non existing face, so we skip this here if necessary
361 if ( !_m.is_boundary(heh) )
362 _m.set_halfedge_handle( _m.face_handle(heh), heh );
363
364 _m.set_halfedge_handle( vh1, opp_new_heh );
365
366 // Never forget this, when playing with the topology
367 _m.adjust_outgoing_halfedge( vh );
368 _m.adjust_outgoing_halfedge( vh1 );
369 }
370
371 private: // geometry helper
372
373 void compute_midpoint(mesh_t& _m, const typename mesh_t::EdgeHandle& _eh)
374 {
375 typename mesh_t::HalfedgeHandle heh, opp_heh;
376
377 heh = _m.halfedge_handle( _eh, 0);
378 opp_heh = _m.halfedge_handle( _eh, 1);
379
380 typename mesh_t::Point pos(0,0,0);
381
382 typename mesh_t::VertexHandle a_0(_m.to_vertex_handle(heh));
383 typename mesh_t::VertexHandle a_1(_m.to_vertex_handle(opp_heh));
384
385 // boundary edge: 4-point scheme
386 if (_m.is_boundary(_eh) )
387 {
388 pos = _m.point(a_0);
389 pos += _m.point(a_1);
390 pos *= static_cast<RealType>(9.0/16.0);
391 typename mesh_t::Point tpos;
392 if(_m.is_boundary(heh))
393 {
394 tpos = _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
395 tpos += _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh))));
396 }
397 else
398 {
399 assert(_m.is_boundary(opp_heh));
400 tpos = _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(opp_heh)));
401 tpos += _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(opp_heh))));
402 }
403 tpos *= static_cast<RealType>(-1.0/16.0);
404 pos += tpos;
405 }
406 else
407 {
408 int valence_a_0 = _m.valence(a_0);
409 int valence_a_1 = _m.valence(a_1);
410 assert(valence_a_0>2);
411 assert(valence_a_1>2);
412
413 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
414 {
415 real_t alpha = real_t(1.0/2);
416 real_t beta = real_t(1.0/8);
417 real_t gamma = real_t(-1.0/16);
418
419 //get points
420 typename mesh_t::VertexHandle b_0, b_1, c_0, c_1, c_2, c_3;
421 typename mesh_t::HalfedgeHandle t_he;
422
423 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(heh));
424 b_0 = _m.to_vertex_handle(t_he);
425 if(!_m.is_boundary(_m.opposite_halfedge_handle(t_he)))
426 {
427 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
428 c_0 = _m.to_vertex_handle(t_he);
429 }
430
431 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh));
432 b_1 = _m.to_vertex_handle(t_he);
433 if(!_m.is_boundary(t_he))
434 {
435 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
436 c_1 = _m.to_vertex_handle(t_he);
437 }
438
439 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(opp_heh));
440 assert(b_1.idx()==_m.to_vertex_handle(t_he).idx());
441 if(!_m.is_boundary(_m.opposite_halfedge_handle(t_he)))
442 {
443 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
444 c_2 = _m.to_vertex_handle(t_he);
445 }
446
447 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(opp_heh));
448 assert(b_0==_m.to_vertex_handle(t_he));
449 if(!_m.is_boundary(t_he))
450 {
451 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
452 c_3 = _m.to_vertex_handle(t_he);
453 }
454
455 //compute position.
456 //a0,a1,b0,b1 must exist.
457 assert(a_0.is_valid());
458 assert(a_1.is_valid());
459 assert(b_0.is_valid());
460 assert(b_1.is_valid());
461 //The other vertices may be created from symmetry is they are on the other side of the boundary.
462
463 pos = _m.point(a_0);
464 pos += _m.point(a_1);
465 pos *= alpha;
466
467 typename mesh_t::Point tpos ( _m.point(b_0) );
468 tpos += _m.point(b_1);
469 tpos *= beta;
470 pos += tpos;
471
472 typename mesh_t::Point pc_0, pc_1, pc_2, pc_3;
473 if(c_0.is_valid())
474 pc_0 = _m.point(c_0);
475 else //create the point by symmetry
476 {
477 pc_0 = _m.point(a_1) + _m.point(b_0) - _m.point(a_0);
478 }
479 if(c_1.is_valid())
480 pc_1 = _m.point(c_1);
481 else //create the point by symmetry
482 {
483 pc_1 = _m.point(a_1) + _m.point(b_1) - _m.point(a_0);
484 }
485 if(c_2.is_valid())
486 pc_2 = _m.point(c_2);
487 else //create the point by symmetry
488 {
489 pc_2 = _m.point(a_0) + _m.point(b_1) - _m.point(a_1);
490 }
491 if(c_3.is_valid())
492 pc_3 = _m.point(c_3);
493 else //create the point by symmetry
494 {
495 pc_3 = _m.point(a_0) + _m.point(b_0) - _m.point(a_1);
496 }
497 tpos = pc_0;
498 tpos += pc_1;
499 tpos += pc_2;
500 tpos += pc_3;
501 tpos *= gamma;
502 pos += tpos;
503 }
504 else //at least one endpoint is [irregular and not in boundary]
505 {
506 RealType normFactor = static_cast<RealType>(0.0);
507
508 if(valence_a_0!=6 && !_m.is_boundary(a_0))
509 {
510 assert((int)weights[valence_a_0].size()==valence_a_0+1);
511 typename mesh_t::HalfedgeHandle t_he = opp_heh;
512 for(int i = 0; i < valence_a_0 ; t_he=_m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he)), ++i)
513 {
514 pos += weights[valence_a_0][i] * _m.point(_m.to_vertex_handle(t_he));
515 }
516 assert(t_he==opp_heh);
517
518 //add irregular vertex:
519 pos += weights[valence_a_0][valence_a_0] * _m.point(a_0);
520 ++normFactor;
521 }
522
523 if(valence_a_1!=6 && !_m.is_boundary(a_1))
524 {
525 assert((int)weights[valence_a_1].size()==valence_a_1+1);
526 typename mesh_t::HalfedgeHandle t_he = heh;
527 for(int i = 0; i < valence_a_1 ; t_he=_m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he)), ++i)
528 {
529 pos += weights[valence_a_1][i] * _m.point(_m.to_vertex_handle(t_he));
530 }
531 assert(t_he==heh);
532 //add irregular vertex:
533 pos += weights[valence_a_1][valence_a_1] * _m.point(a_1);
534 ++normFactor;
535 }
536
537 assert(normFactor>0.1); //normFactor should be 1 or 2
538
539 //if both vertices are irregular, average positions:
540 pos /= normFactor;
541 }
542 }
543 _m.property( ep_pos_, _eh ) = pos;
544 }
545
546 private: // data
547
548 OpenMesh::VPropHandleT< typename mesh_t::Point > vp_pos_;
549 OpenMesh::EPropHandleT< typename mesh_t::Point > ep_pos_;
550
551 weights_t weights;
552
553 };
554
555 } // END_NS_UNIFORM
556 } // END_NS_SUBDIVIDER
557 } // END_NS_OPENMESH
558 #endif
559
OpenMesh::VPropHandleT< typename mesh_t::Point >
SubdividerT.hh
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::EPropHandleT< typename mesh_t::Point >
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::name
const char * name() const override
Return name of subdivision algorithm.
Definition: ModifiedButterFlyT.hh:120
OpenMesh::Utils::MeshCheckerT
Definition: MeshCheckerT.hh:73
OpenMesh
Definition: MeshItems.hh:59
OpenMesh::Subdivider::Uniform::SubdividerT
Definition: SubdividerT.hh:88
OpenMesh::Subdivider::Uniform::ModifiedButterflyT
Definition: ModifiedButterFlyT.hh:92
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::cleanup
bool cleanup(mesh_t &_m) override
Cleanup mesh after usage, e.g. remove added properties.
Definition: ModifiedButterFlyT.hh:169