58 #ifndef SP_MODIFIED_BUTTERFLY_H 59 #define SP_MODIFIED_BUTTERFLY_H 62 #include <OpenMesh/Core/Utils/vector_cast.hh> 63 #include <OpenMesh/Core/Utils/Property.hh> 66 #if defined(OM_CC_MIPS) 76 namespace Subdivider {
91 template <
typename MeshType,
typename RealType =
double>
96 typedef RealType real_t;
97 typedef MeshType mesh_t;
100 typedef std::vector< std::vector<real_t> > weights_t;
101 typedef std::vector<real_t> weight_t;
120 const char *
name()
const override {
return "Uniform Spectral"; }
126 weights.resize(_max_valence);
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;
135 weights[4].resize(5);
136 weights[4][0] = real_t(3.0)/8;
138 weights[4][2] = real_t(-1.0)/8;
140 weights[4][4] = real_t(3.0)/4;
142 for(
unsigned int K = 5; K<_max_valence; ++K)
144 weights[K].resize(K+1);
146 double invK = 1.0/
static_cast<double>(K);
148 for(
unsigned int j=0; j<K; ++j)
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];
153 weights[K][K] =
static_cast<real_t
>(1.0) - sum;
163 _m.add_property( vp_pos_ );
164 _m.add_property( ep_pos_ );
171 _m.remove_property( vp_pos_ );
172 _m.remove_property( ep_pos_ );
177 bool subdivide( MeshType& _m,
size_t _n ,
const bool _update_points =
true)
override 182 typename mesh_t::FaceIter fit, f_end;
183 typename mesh_t::EdgeIter eit, e_end;
184 typename mesh_t::VertexIter vit;
187 for (
size_t i=0; i < _n; ++i)
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);
196 for (eit=_m.edges_begin(); eit != _m.edges_end(); ++eit)
197 compute_midpoint( _m, *eit );
204 e_end = _m.edges_end();
205 for (eit=_m.edges_begin(); eit != e_end; ++eit)
206 split_edge(_m, *eit );
212 f_end = _m.faces_end();
213 for (fit = _m.faces_begin(); fit != f_end; ++fit)
214 split_face(_m, *fit );
218 for ( vit = _m.vertices_begin();
219 vit != _m.vertices_end(); ++vit)
220 _m.set_point(*vit, _m.property( vp_pos_, *vit ) );
222 #if defined(_DEBUG) || defined(DEBUG) 233 void split_face(mesh_t& _m,
const typename mesh_t::FaceHandle& _fh)
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)));
241 corner_cutting( _m, heh1 );
242 corner_cutting( _m, heh2 );
243 corner_cutting( _m, heh3 );
247 void corner_cutting(mesh_t& _m,
const typename mesh_t::HalfedgeHandle& _he)
250 typename mesh_t::HalfedgeHandle
253 heh6(_m.next_halfedge_handle(heh1));
256 for (; _m.next_halfedge_handle(_m.next_halfedge_handle(heh5)) != heh1;
257 heh5 = _m.next_halfedge_handle(heh5))
260 typename mesh_t::VertexHandle
261 vh1 = _m.to_vertex_handle(heh1),
262 vh2 = _m.to_vertex_handle(heh5);
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));
283 typename mesh_t::FaceHandle fh_old(_m.face_handle(heh6));
284 typename mesh_t::FaceHandle fh_new(_m.new_face());
288 _m.set_next_halfedge_handle(heh4, heh6);
289 _m.set_next_halfedge_handle(heh5, heh4);
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);
297 _m.set_next_halfedge_handle(heh1, heh3);
298 _m.set_next_halfedge_handle(heh3, heh2);
300 _m.set_face_handle(heh1, fh_new);
301 _m.set_face_handle(heh2, fh_new);
302 _m.set_face_handle(heh3, fh_new);
304 _m.set_halfedge_handle(fh_new, heh1);
308 void split_edge(mesh_t& _m,
const typename mesh_t::EdgeHandle& _eh)
310 typename mesh_t::HalfedgeHandle
311 heh = _m.halfedge_handle(_eh, 0),
312 opp_heh = _m.halfedge_handle(_eh, 1);
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);
320 vh = _m.new_vertex( zero );
323 _m.property( vp_pos_, vh ) = _m.property( ep_pos_, _eh );
327 if (_m.is_boundary(_eh))
330 _m.next_halfedge_handle(t_heh) != opp_heh;
331 t_heh = _m.opposite_halfedge_handle(_m.next_halfedge_handle(t_heh)))
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) )
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 );
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);
351 if (_m.face_handle(opp_heh).is_valid())
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);
357 _m.set_face_handle( new_heh, _m.face_handle(heh) );
358 _m.set_halfedge_handle( vh, new_heh);
361 if ( !_m.is_boundary(heh) )
362 _m.set_halfedge_handle( _m.face_handle(heh), heh );
364 _m.set_halfedge_handle( vh1, opp_new_heh );
367 _m.adjust_outgoing_halfedge( vh );
368 _m.adjust_outgoing_halfedge( vh1 );
373 void compute_midpoint(mesh_t& _m,
const typename mesh_t::EdgeHandle& _eh)
375 typename mesh_t::HalfedgeHandle heh, opp_heh;
377 heh = _m.halfedge_handle( _eh, 0);
378 opp_heh = _m.halfedge_handle( _eh, 1);
380 typename mesh_t::Point pos(0,0,0);
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));
386 if (_m.is_boundary(_eh) )
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))
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))));
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))));
403 tpos *=
static_cast<RealType
>(-1.0/16.0);
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);
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)) )
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);
420 typename mesh_t::VertexHandle b_0, b_1, c_0, c_1, c_2, c_3;
421 typename mesh_t::HalfedgeHandle t_he;
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)))
427 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
428 c_0 = _m.to_vertex_handle(t_he);
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))
435 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
436 c_1 = _m.to_vertex_handle(t_he);
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)))
443 t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
444 c_2 = _m.to_vertex_handle(t_he);
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))
451 t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
452 c_3 = _m.to_vertex_handle(t_he);
457 assert(a_0.is_valid());
458 assert(a_1.is_valid());
459 assert(b_0.is_valid());
460 assert(b_1.is_valid());
464 pos += _m.point(a_1);
467 typename mesh_t::Point tpos ( _m.point(b_0) );
468 tpos += _m.point(b_1);
472 typename mesh_t::Point pc_0, pc_1, pc_2, pc_3;
474 pc_0 = _m.point(c_0);
477 pc_0 = _m.point(a_1) + _m.point(b_0) - _m.point(a_0);
480 pc_1 = _m.point(c_1);
483 pc_1 = _m.point(a_1) + _m.point(b_1) - _m.point(a_0);
486 pc_2 = _m.point(c_2);
489 pc_2 = _m.point(a_0) + _m.point(b_1) - _m.point(a_1);
492 pc_3 = _m.point(c_3);
495 pc_3 = _m.point(a_0) + _m.point(b_0) - _m.point(a_1);
506 RealType normFactor =
static_cast<RealType
>(0.0);
508 if(valence_a_0!=6 && !_m.is_boundary(a_0))
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)
514 pos += weights[valence_a_0][i] * _m.point(_m.to_vertex_handle(t_he));
516 assert(t_he==opp_heh);
519 pos += weights[valence_a_0][valence_a_0] * _m.point(a_0);
523 if(valence_a_1!=6 && !_m.is_boundary(a_1))
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)
529 pos += weights[valence_a_1][i] * _m.point(_m.to_vertex_handle(t_he));
533 pos += weights[valence_a_1][valence_a_1] * _m.point(a_1);
537 assert(normFactor>0.1);
543 _m.property( ep_pos_, _eh ) = pos;