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