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 
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;
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 
543 
544  weights_t weights;
545 
546 };
547 
548 } // END_NS_UNIFORM
549 } // END_NS_SUBDIVIDER
550 } // END_NS_OPENMESH
551 #endif
552 
bool subdivide(MeshType &_m, size_t _n, const bool _update_points=true) override
Subdivide mesh _m _n times.
Definition: ModifiedButterFlyT.hh:177
Modified Butterfly subdivision algorithm.
Definition: ModifiedButterFlyT.hh:92
bool cleanup(mesh_t &_m) override
Cleanup mesh after usage, e.g. remove added properties.
Definition: ModifiedButterFlyT.hh:169
const char * name() const override
Return name of subdivision algorithm.
Definition: ModifiedButterFlyT.hh:120
Abstract base class for uniform subdivision algorithms.
Definition: SubdividerT.hh:88
bool prepare(mesh_t &_m) override
Prepare mesh, e.g. add properties.
Definition: ModifiedButterFlyT.hh:161
void init_weights(size_t _max_valence=30)
Pre-compute weights.
Definition: ModifiedButterFlyT.hh:124
Check integrity of mesh.
Definition: MeshCheckerT.hh:73
Contains all the mesh ingredients like the polygonal mesh, the triangle mesh, different mesh kernels ...
Definition: MeshItems.hh:59
Kernel::VertexHandle VertexHandle
Handle for referencing the corresponding item.
Definition: PolyMeshT.hh:136

Project OpenMesh, ©  Computer Graphics Group, RWTH Aachen. Documentation generated using doxygen .