Pixyco-XYZ/openmesh
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added Sqrt3InterpolatingSubdividerLabsikGreinerT and ModifiedButterFlyT (Thanks to Clément Courbet for providing the code)

Browse Source 
git-svn-id: http://www.openmesh.org/svnrepo/OpenMesh/trunk@341 fdac6126-5c0c-442c-9429-916003d36597
This commit is contained in:
Jan Möbius
2010-11-15 08:39:00 +00:00
parent eba601f190
commit 0a9f2815ae
2 changed files with 1146 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

543
src/OpenMesh/Tools/Subdivider/Uniform/ModifiedButterFlyT.hh Normal file
View File

@@ -0,0 +1,543 @@
/*===========================================================================*\
* *
* OpenMesh *
* Copyright (C) 2001-2010 by Computer Graphics Group, RWTH Aachen *
* www.openmesh.org *
* *
*---------------------------------------------------------------------------*
* This file is part of OpenMesh. *
* *
* OpenMesh is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as *
* published by the Free Software Foundation, either version 3 of *
* the License, or (at your option) any later version with the *
* following exceptions: *
* *
* If other files instantiate templates or use macros *
* or inline functions from this file, or you compile this file and *
* link it with other files to produce an executable, this file does *
* not by itself cause the resulting executable to be covered by the *
* GNU Lesser General Public License. This exception does not however *
* invalidate any other reasons why the executable file might be *
* covered by the GNU Lesser General Public License. *
* *
* OpenMesh is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU LesserGeneral Public *
* License along with OpenMesh. If not, *
* see <http://www.gnu.org/licenses/>. *
* *
\*==========================================================================*/
/*==========================================================================*\
* *
* $Revision: 410 $ *
* $Date: 2010-06-17 12:45:58 +0200 (Do, 17. Jun 2010) $ *
* *
\*==========================================================================*/
/** \file ModifiedButterflyT.hh
The modified butterfly scheme of Denis Zorin, Peter Schröder and Wim Sweldens,
``Interpolating subdivision for meshes with arbitrary topology,'' in Proceedings
of SIGGRAPH 1996, ACM SIGGRAPH, 1996, pp. 189-192.
Clement Courbet - clement.courbet@ecp.fr
*/
//=============================================================================
//
// CLASS ModifiedButterflyT
//
//=============================================================================
#ifndef SP_MODIFIED_BUTTERFLY_H
#define SP_MODIFIED_BUTTERFLY_H
#include <OpenMesh/Tools/Subdivider/Uniform/SubdividerT.hh>
#include <OpenMesh/Core/Utils/vector_cast.hh>
#include <OpenMesh/Core/Utils/Property.hh>
// -------------------- STL
#include <vector>
#if defined(OM_CC_MIPS)
# include <math.h>
#else
# include <cmath>
#endif
//== NAMESPACE ================================================================
namespace OpenMesh { // BEGIN_NS_OPENMESH
namespace Subdivider { // BEGIN_NS_DECIMATER
namespace Uniform { // BEGIN_NS_UNIFORM
//== CLASS DEFINITION =========================================================
template <typename MeshType, typename RealType = float>
class ModifiedButterflyT : public SubdividerT<MeshType, RealType>
{
public:
typedef RealType real_t;
typedef MeshType mesh_t;
typedef SubdividerT< mesh_t, real_t > parent_t;
typedef std::vector< std::vector<real_t> > weights_t;
typedef std::vector<real_t> weight_t;
public:
ModifiedButterflyT() : parent_t()
{ init_weights(); }
ModifiedButterflyT( mesh_t& _m) : parent_t(_m)
{ init_weights(); }
~ModifiedButterflyT() {}
public:
const char *name() const { return "Uniform Spectral"; }
/// Pre-compute weights
void init_weights(size_t _max_valence=20)
{
weights.resize(_max_valence);
//special case: K==3, K==4
weights[3].resize(4);
weights[3][0] = real_t(5.0)/12;
weights[3][1] = real_t(-1.0)/12;
weights[3][2] = real_t(-1.0)/12;
weights[3][3] = real_t(3.0)/4;
weights[4].resize(5);
weights[4][0] = real_t(3.0)/8;
weights[4][1] = 0;
weights[4][2] = real_t(-1.0)/8;
weights[4][3] = 0;
weights[4][4] = real_t(3.0)/4;
for(unsigned int K = 5; K<_max_valence; ++K)
{
weights[K].resize(K+1);
// s(j) = ( 1/4 + cos(2*pi*j/K) + 1/2 * cos(4*pi*j/K) )/K
real_t invK = 1.0/real_t(K);
real_t sum = 0;
for(unsigned int j=0; j<K; ++j)
{
weights[K][j] = (0.25 + cos(2.0*M_PI*j*invK) + 0.5*cos(4.0*M_PI*j*invK))*invK;
sum += weights[K][j];
}
weights[K][K] = (real_t)1.0 - sum;
}
}
protected:
bool prepare( mesh_t& _m )
{
_m.add_property( vp_pos_ );
_m.add_property( ep_pos_ );
return true;
}
bool cleanup( mesh_t& _m )
{
_m.remove_property( vp_pos_ );
_m.remove_property( ep_pos_ );
return true;
}
bool subdivide( mesh_t& _m, size_t _n)
{
typename mesh_t::FaceIter fit, f_end;
typename mesh_t::EdgeIter eit, e_end;
typename mesh_t::VertexIter vit;
// Do _n subdivisions
for (size_t i=0; i < _n; ++i)
{
// This is an interpolating scheme, old vertices remain the same.
typename mesh_t::VertexIter initialVerticesEnd = _m.vertices_end();
for ( vit = _m.vertices_begin(); vit != initialVerticesEnd; ++vit)
_m.property( vp_pos_, vit.handle() ) = _m.point(vit.handle());
// Compute position for new vertices and store them in the edge property
for (eit=_m.edges_begin(); eit != _m.edges_end(); ++eit)
compute_midpoint( _m, eit.handle() );
// Split each edge at midpoint and store precomputed positions (stored in
// edge property ep_pos_) in the vertex property vp_pos_;
// Attention! Creating new edges, hence make sure the loop ends correctly.
e_end = _m.edges_end();
for (eit=_m.edges_begin(); eit != e_end; ++eit)
split_edge(_m, eit.handle() );
// Commit changes in topology and reconsitute consistency
// Attention! Creating new faces, hence make sure the loop ends correctly.
f_end = _m.faces_end();
for (fit = _m.faces_begin(); fit != f_end; ++fit)
split_face(_m, fit.handle() );
// Commit changes in geometry
for ( vit = /*initialVerticesEnd;*/_m.vertices_begin();
vit != _m.vertices_end(); ++vit)
_m.set_point(vit, _m.property( vp_pos_, vit ) );
#if defined(_DEBUG) || defined(DEBUG)
// Now we have an consistent mesh!
assert( OpenMesh::Utils::MeshCheckerT<mesh_t>(_m).check() );
#endif
}
return true;
}
private: // topological modifiers
void split_face(mesh_t& _m, const typename mesh_t::FaceHandle& _fh)
{
typename mesh_t::HalfedgeHandle
heh1(_m.halfedge_handle(_fh)),
heh2(_m.next_halfedge_handle(_m.next_halfedge_handle(heh1))),
heh3(_m.next_halfedge_handle(_m.next_halfedge_handle(heh2)));
// Cutting off every corner of the 6_gon
corner_cutting( _m, heh1 );
corner_cutting( _m, heh2 );
corner_cutting( _m, heh3 );
}
void corner_cutting(mesh_t& _m, const typename mesh_t::HalfedgeHandle& _he)
{
// Define Halfedge Handles
typename mesh_t::HalfedgeHandle
heh1(_he),
heh5(heh1),
heh6(_m.next_halfedge_handle(heh1));
// Cycle around the polygon to find correct Halfedge
for (; _m.next_halfedge_handle(_m.next_halfedge_handle(heh5)) != heh1;
heh5 = _m.next_halfedge_handle(heh5))
{}
typename mesh_t::VertexHandle
vh1 = _m.to_vertex_handle(heh1),
vh2 = _m.to_vertex_handle(heh5);
typename mesh_t::HalfedgeHandle
heh2(_m.next_halfedge_handle(heh5)),
heh3(_m.new_edge( vh1, vh2)),
heh4(_m.opposite_halfedge_handle(heh3));
/* Intermediate result
*
* *
* 5 /|\
* /_ \
* vh2> * *
* /|\3 |\
* /_ \|4 \
* *----\*----\*
* 1 ^ 6
* vh1 (adjust_outgoing halfedge!)
*/
// Old and new Face
typename mesh_t::FaceHandle fh_old(_m.face_handle(heh6));
typename mesh_t::FaceHandle fh_new(_m.new_face());
// Re-Set Handles around old Face
_m.set_next_halfedge_handle(heh4, heh6);
_m.set_next_halfedge_handle(heh5, heh4);
_m.set_face_handle(heh4, fh_old);
_m.set_face_handle(heh5, fh_old);
_m.set_face_handle(heh6, fh_old);
_m.set_halfedge_handle(fh_old, heh4);
// Re-Set Handles around new Face
_m.set_next_halfedge_handle(heh1, heh3);
_m.set_next_halfedge_handle(heh3, heh2);
_m.set_face_handle(heh1, fh_new);
_m.set_face_handle(heh2, fh_new);
_m.set_face_handle(heh3, fh_new);
_m.set_halfedge_handle(fh_new, heh1);
}
void split_edge(mesh_t& _m, const typename mesh_t::EdgeHandle& _eh)
{
typename mesh_t::HalfedgeHandle
heh = _m.halfedge_handle(_eh, 0),
opp_heh = _m.halfedge_handle(_eh, 1);
typename mesh_t::HalfedgeHandle new_heh, opp_new_heh, t_heh;
typename mesh_t::VertexHandle vh;
typename mesh_t::VertexHandle vh1(_m.to_vertex_handle(heh));
typename mesh_t::Point zero(0,0,0);
// new vertex
vh = _m.new_vertex( zero );
// memorize position, will be set later
_m.property( vp_pos_, vh ) = _m.property( ep_pos_, _eh );
// Re-link mesh entities
if (_m.is_boundary(_eh))
{
for (t_heh = heh;
_m.next_halfedge_handle(t_heh) != opp_heh;
t_heh = _m.opposite_halfedge_handle(_m.next_halfedge_handle(t_heh)))
{}
}
else
{
for (t_heh = _m.next_halfedge_handle(opp_heh);
_m.next_halfedge_handle(t_heh) != opp_heh;
t_heh = _m.next_halfedge_handle(t_heh) )
{}
}
new_heh = _m.new_edge(vh, vh1);
opp_new_heh = _m.opposite_halfedge_handle(new_heh);
_m.set_vertex_handle( heh, vh );
_m.set_next_halfedge_handle(t_heh, opp_new_heh);
_m.set_next_halfedge_handle(new_heh, _m.next_halfedge_handle(heh));
_m.set_next_halfedge_handle(heh, new_heh);
_m.set_next_halfedge_handle(opp_new_heh, opp_heh);
if (_m.face_handle(opp_heh).is_valid())
{
_m.set_face_handle(opp_new_heh, _m.face_handle(opp_heh));
_m.set_halfedge_handle(_m.face_handle(opp_new_heh), opp_new_heh);
}
_m.set_face_handle( new_heh, _m.face_handle(heh) );
_m.set_halfedge_handle( vh, new_heh);
_m.set_halfedge_handle( _m.face_handle(heh), heh );
_m.set_halfedge_handle( vh1, opp_new_heh );
// Never forget this, when playing with the topology
_m.adjust_outgoing_halfedge( vh );
_m.adjust_outgoing_halfedge( vh1 );
}
private: // geometry helper
void compute_midpoint(mesh_t& _m, const typename mesh_t::EdgeHandle& _eh)
{
typename mesh_t::HalfedgeHandle heh, opp_heh;
heh = _m.halfedge_handle( _eh, 0);
opp_heh = _m.halfedge_handle( _eh, 1);
typename mesh_t::Point pos(0,0,0);
typename mesh_t::VertexHandle a_0(_m.to_vertex_handle(heh));
typename mesh_t::VertexHandle a_1(_m.to_vertex_handle(opp_heh));
// boundary edge: 4-point scheme
if (_m.is_boundary(_eh) )
{
pos = _m.point(a_0);
pos += _m.point(a_1);
pos *= 9.0/16;
typename mesh_t::Point tpos;
if(_m.is_boundary(heh))
{
tpos = _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
tpos += _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh))));
}
else
{
assert(_m.is_boundary(opp_heh));
tpos = _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(opp_heh)));
tpos += _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(opp_heh))));
}
tpos *= -1.0/16;
pos += tpos;
}
else
{
int valence_a_0 = _m.valence(a_0);
int valence_a_1 = _m.valence(a_1);
assert(valence_a_0>2);
assert(valence_a_1>2);
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
{
real_t alpha = real_t(1.0/2);
real_t beta = real_t(1.0/8);
real_t gamma = real_t(-1.0/16);
//get points
typename mesh_t::VertexHandle b_0, b_1, c_0, c_1, c_2, c_3;
typename mesh_t::HalfedgeHandle t_he;
t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(heh));
b_0 = _m.to_vertex_handle(t_he);
if(!_m.is_boundary(_m.opposite_halfedge_handle(t_he)))
{
t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
c_0 = _m.to_vertex_handle(t_he);
}
t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh));
b_1 = _m.to_vertex_handle(t_he);
if(!_m.is_boundary(t_he))
{
t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
c_1 = _m.to_vertex_handle(t_he);
}
t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(opp_heh));
assert(b_1.idx()==_m.to_vertex_handle(t_he).idx());
if(!_m.is_boundary(_m.opposite_halfedge_handle(t_he)))
{
t_he = _m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he));
c_2 = _m.to_vertex_handle(t_he);
}
t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(opp_heh));
assert(b_0==_m.to_vertex_handle(t_he));
if(!_m.is_boundary(t_he))
{
t_he = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(t_he));
c_3 = _m.to_vertex_handle(t_he);
}
//compute position.
//a0,a1,b0,b1 must exist.
assert(a_0.is_valid());
assert(a_1.is_valid());
assert(b_0.is_valid());
assert(b_1.is_valid());
//The other vertices may be created from symmetry is they are on the other side of the boundary.
pos = _m.point(a_0);
pos += _m.point(a_1);
pos *= alpha;
typename mesh_t::Point tpos ( _m.point(b_0) );
tpos += _m.point(b_1);
tpos *= beta;
pos += tpos;
typename mesh_t::Point pc_0, pc_1, pc_2, pc_3;
if(c_0.is_valid())
pc_0 = _m.point(c_0);
else //create the point by symmetry
{
pc_0 = _m.point(a_1) + _m.point(b_0) - _m.point(a_0);
}
if(c_1.is_valid())
pc_1 = _m.point(c_1);
else //create the point by symmetry
{
pc_1 = _m.point(a_1) + _m.point(b_1) - _m.point(a_0);
}
if(c_2.is_valid())
pc_2 = _m.point(c_2);
else //create the point by symmetry
{
pc_2 = _m.point(a_0) + _m.point(b_1) - _m.point(a_1);
}
if(c_3.is_valid())
pc_3 = _m.point(c_3);
else //create the point by symmetry
{
pc_3 = _m.point(a_0) + _m.point(b_0) - _m.point(a_1);
}
tpos = pc_0;
tpos += pc_1;
tpos += pc_2;
tpos += pc_3;
tpos *= gamma;
pos += tpos;
}
else //at least one endpoint is [irregular and not in boundary]
{
double normFactor = 0.0;
if(valence_a_0!=6 && !_m.is_boundary(a_0))
{
assert((int)weights[valence_a_0].size()==valence_a_0+1);
typename mesh_t::HalfedgeHandle t_he = opp_heh;
for(int i = 0; i < valence_a_0 ; t_he=_m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he)), ++i)
{
pos += weights[valence_a_0][i] * _m.point(_m.to_vertex_handle(t_he));
}
assert(t_he==opp_heh);
//add irregular vertex:
pos += weights[valence_a_0][valence_a_0] * _m.point(a_0);
++normFactor;
}
if(valence_a_1!=6 && !_m.is_boundary(a_1))
{
assert((int)weights[valence_a_1].size()==valence_a_1+1);
typename mesh_t::HalfedgeHandle t_he = heh;
for(int i = 0; i < valence_a_1 ; t_he=_m.next_halfedge_handle(_m.opposite_halfedge_handle(t_he)), ++i)
{
pos += weights[valence_a_1][i] * _m.point(_m.to_vertex_handle(t_he));
}
assert(t_he==heh);
//add irregular vertex:
pos += weights[valence_a_1][valence_a_1] * _m.point(a_1);
++normFactor;
}
assert(normFactor>0.1); //normFactor should be 1 or 2
//if both vertices are irregular, average positions:
pos /= normFactor;
}
}
_m.property( ep_pos_, _eh ) = pos;
}
private: // data
OpenMesh::VPropHandleT< typename mesh_t::Point > vp_pos_;
OpenMesh::EPropHandleT< typename mesh_t::Point > ep_pos_;
weights_t weights;
};
} // END_NS_UNIFORM
} // END_NS_SUBDIVIDER
} // END_NS_OPENMESH
#endif

603
src/OpenMesh/Tools/Subdivider/Uniform/Sqrt3InterpolatingSubdividerLabsikGreinerT.hh Normal file
View File

@@ -0,0 +1,603 @@
/*===========================================================================*\
* *
* OpenMesh *
* Copyright (C) 2001-2010 by Computer Graphics Group, RWTH Aachen *
* www.openmesh.org *
* *
*---------------------------------------------------------------------------*
* This file is part of OpenMesh. *
* *
* OpenMesh is free software: you can redistribute it and/or modify *
* it under the terms of the GNU Lesser General Public License as *
* published by the Free Software Foundation, either version 3 of *
* the License, or (at your option) any later version with the *
* following exceptions: *
* *
* If other files instantiate templates or use macros *
* or inline functions from this file, or you compile this file and *
* link it with other files to produce an executable, this file does *
* not by itself cause the resulting executable to be covered by the *
* GNU Lesser General Public License. This exception does not however *
* invalidate any other reasons why the executable file might be *
* covered by the GNU Lesser General Public License. *
* *
* OpenMesh is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU Lesser General Public License for more details. *
* *
* You should have received a copy of the GNU LesserGeneral Public *
* License along with OpenMesh. If not, *
* see <http://www.gnu.org/licenses/>. *
* *
\*==========================================================================*/
/*==========================================================================*\
* *
* $Revision: 410 $ *
* $Date: 2010-06-17 12:45:58 +0200 (Do, 17. Jun 2010) $ *
* *
\*==========================================================================*/
/** \file Sqrt3InterpolatingSubdividerLabsikGreinerT.hh
Interpolating Labsik Greiner Subdivider as described in "interpolating sqrt(3) subdivision" Labsik & Greiner, 2000
clement.courbet@ecp.fr
*/
//=============================================================================
//
// CLASS InterpolatingSqrt3LGT
//
//=============================================================================
#ifndef OPENMESH_SUBDIVIDER_UNIFORM_INTERP_SQRT3T_LABSIK_GREINER_HH
#define OPENMESH_SUBDIVIDER_UNIFORM_INTERP_SQRT3T_LABSIK_GREINER_HH
//== INCLUDES =================================================================
#include <OpenMesh/Core/Mesh/Handles.hh>
#include <OpenMesh/Core/System/config.hh>
#include <OpenMesh/Tools/Subdivider/Uniform/SubdividerT.hh>
#if defined(_DEBUG) || defined(DEBUG)
// Makes life lot easier, when playing/messing around with low-level topology
// changing methods of OpenMesh
# include <OpenMesh/Tools/Utils/MeshCheckerT.hh>
# define ASSERT_CONSISTENCY( T, m ) \
assert(OpenMesh::Utils::MeshCheckerT<T>(m).check())
#else
# define ASSERT_CONSISTENCY( T, m )
#endif
// -------------------- STL
#include <vector>
#if defined(OM_CC_MIPS)
# include <math.h>
#else
# include <cmath>
#endif
//#define MIRROR_TRIANGLES
//#define MIN_NORM
//== NAMESPACE ================================================================
namespace OpenMesh { // BEGIN_NS_OPENMESH
namespace Subdivider { // BEGIN_NS_DECIMATER
namespace Uniform { // BEGIN_NS_UNIFORM
//== CLASS DEFINITION =========================================================
/** %Uniform Interpolating Sqrt3 subdivision algorithm
*
*/
template <typename MeshType, typename RealType = float>
class InterpolatingSqrt3LGT : public SubdividerT< MeshType, RealType >
{
public:
typedef RealType real_t;
typedef MeshType mesh_t;
typedef SubdividerT< mesh_t, real_t > parent_t;
typedef std::vector< std::vector<real_t> > weights_t;
public:
InterpolatingSqrt3LGT(void) : parent_t()
{ init_weights(); }
InterpolatingSqrt3LGT(MeshType &_m) : parent_t(_m)
{ init_weights(); }
virtual ~InterpolatingSqrt3LGT() {}
public:
const char *name() const { return "Uniform Interpolating Sqrt3"; }
/// Pre-compute weights
void init_weights(size_t _max_valence=50)
{
weights_.resize(_max_valence);
weights_[3].resize(4);
weights_[3][0] = +4.0/27;
weights_[3][1] = -5.0/27;
weights_[3][2] = +4.0/27;
weights_[3][3] = +8.0/9;
weights_[4].resize(5);
weights_[4][0] = +2.0/9;
weights_[4][1] = -1.0/9;
weights_[4][2] = -1.0/9;
weights_[4][3] = +2.0/9;
weights_[4][4] = +7.0/9 ;
for(unsigned int K=5; K<_max_valence; ++K)
{
weights_[K].resize(K+1);
double aH = 2.0*cos(M_PI/K)/3.0;
weights_[K][K] = 1.0 - aH*aH;
for(unsigned int i=0; i<K; ++i)
{
weights_[K][i] = (aH*aH + 2.0*aH*cos(2.0*i*M_PI/K + M_PI/K) + 2.0*aH*aH*cos(4.0*i*M_PI/K + 2.0*M_PI/K))/K;
}
}
//just to be sure:
weights_[6].resize(0);
}
protected:
bool prepare( MeshType& _m )
{
_m.request_edge_status();
_m.add_property( fp_pos_ );
_m.add_property( ep_nv_ );
_m.add_property( mp_gen_ );
_m.property( mp_gen_ ) = 0;
return _m.has_edge_status()
&& ep_nv_.is_valid() && mp_gen_.is_valid();
}
bool cleanup( MeshType& _m )
{
_m.release_edge_status();
_m.remove_property( fp_pos_ );
_m.remove_property( ep_nv_ );
_m.remove_property( mp_gen_ );
return true;
}
bool subdivide( MeshType& _m, size_t _n )
{
typename MeshType::VertexIter vit;
typename MeshType::VertexVertexIter vvit;
typename MeshType::EdgeIter eit;
typename MeshType::FaceIter fit;
typename MeshType::FaceVertexIter fvit;
typename MeshType::FaceHalfedgeIter fheit;
typename MeshType::VertexHandle vh;
typename MeshType::HalfedgeHandle heh;
typename MeshType::Point pos(0,0,0), zero(0,0,0);
size_t &gen = _m.property( mp_gen_ );
for (size_t l=0; l<_n; ++l)
{
// tag existing edges
for (eit=_m.edges_begin(); eit != _m.edges_end();++eit)
{
_m.status( eit ).set_tagged( true );
if ( (gen%2) && _m.is_boundary(eit) )
compute_new_boundary_points( _m, eit ); // *) creates new vertices
}
// insert new vertices, and store pos in vp_pos_
typename MeshType::FaceIter fend = _m.faces_end();
for (fit = _m.faces_begin();fit != fend; ++fit)
{
if (_m.is_boundary(fit))
{
if(gen%2)
_m.property(fp_pos_, fit.handle()).invalidate();
else
{
//find the interior boundary halfedge
for( heh = _m.halfedge_handle(fit.handle()); !_m.is_boundary( _m.opposite_halfedge_handle(heh) ); heh = _m.next_halfedge_handle(heh) )
;
assert(_m.is_boundary( _m.opposite_halfedge_handle(heh) ));
pos = zero;
//check for two boundaries case:
if( _m.is_boundary(_m.next_halfedge_handle(heh)) || _m.is_boundary(_m.prev_halfedge_handle(heh)) )
{
if(_m.is_boundary(_m.prev_halfedge_handle(heh)))
heh = _m.prev_halfedge_handle(heh); //ensure that the boundary halfedges are heh and heh->next
//check for three boundaries case:
if(_m.is_boundary(_m.next_halfedge_handle(_m.next_halfedge_handle(heh))))
{
//three boundaries, use COG of triangle
pos += real_t(1.0/3) * _m.point(_m.to_vertex_handle(heh));
pos += real_t(1.0/3) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
pos += real_t(1.0/3) * _m.point(_m.to_vertex_handle(_m.prev_halfedge_handle(heh)));
}
else
{
#ifdef MIRROR_TRIANGLES
//two boundaries, mirror two triangles
pos += real_t(2.0/9) * _m.point(_m.to_vertex_handle(heh));
pos += real_t(4.0/9) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
pos += real_t(4.0/9) * _m.point(_m.to_vertex_handle(_m.prev_halfedge_handle(heh)));
pos += real_t(-1.0/9) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh)))));
#else
pos += real_t(7.0/24) * _m.point(_m.to_vertex_handle(heh));
pos += real_t(3.0/8) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
pos += real_t(3.0/8) * _m.point(_m.to_vertex_handle(_m.prev_halfedge_handle(heh)));
pos += real_t(-1.0/24) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh)))));
#endif
}
}
else
{
vh = _m.to_vertex_handle(_m.next_halfedge_handle(heh));
//check last vertex regularity
if((_m.valence(vh) == 6) || _m.is_boundary(vh))
{
#ifdef MIRROR_TRIANGLES
//use regular rule and mirror one triangle
pos += real_t(5.0/9) * _m.point(vh);
pos += real_t(3.0/9) * _m.point(_m.to_vertex_handle(heh));
pos += real_t(3.0/9) * _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(heh)));
pos += real_t(-1.0/9) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(_m.opposite_halfedge_handle(_m.next_halfedge_handle(heh)))));
pos += real_t(-1.0/9) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh)))));
#else
#ifdef MIN_NORM
pos += real_t(1.0/9) * _m.point(vh);
pos += real_t(1.0/3) * _m.point(_m.to_vertex_handle(heh));
pos += real_t(1.0/3) * _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(heh)));
pos += real_t(1.0/9) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(_m.opposite_halfedge_handle(_m.next_halfedge_handle(heh)))));
pos += real_t(1.0/9) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh)))));
#else
pos += real_t(1.0/2) * _m.point(vh);
pos += real_t(1.0/3) * _m.point(_m.to_vertex_handle(heh));
pos += real_t(1.0/3) * _m.point(_m.to_vertex_handle(_m.opposite_halfedge_handle(heh)));
pos += real_t(-1.0/12) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(_m.opposite_halfedge_handle(_m.next_halfedge_handle(heh)))));
pos += real_t(-1.0/12) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(_m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh)))));
#endif
#endif
}
else
{
//irregular setting, use usual irregular rule
unsigned int K = _m.valence(vh);
pos += weights_[K][K]*_m.point(vh);
heh = _m.opposite_halfedge_handle( _m.next_halfedge_handle(heh) );
for(unsigned int i = 0; i<K; ++i, heh = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh)) )
{
pos += weights_[K][i]*_m.point(_m.to_vertex_handle(heh));
}
}
}
vh = _m.add_vertex( pos );
_m.property(fp_pos_, fit.handle()) = vh;
}
}
else
{
pos = zero;
int nOrdinary = 0;
//check number of extraordinary vertices
for(fvit = _m.fv_iter( fit ); fvit; ++fvit)
if( (_m.valence(fvit.handle())) == 6 || _m.is_boundary(fvit.handle()) )
++nOrdinary;
if(nOrdinary==3)
{
for(fheit = _m.fh_iter( fit ); fheit; ++fheit)
{
//one ring vertex has weight 32/81
heh = fheit.handle();
assert(_m.to_vertex_handle(heh).is_valid());
pos += real_t(32.0/81) * _m.point(_m.to_vertex_handle(heh));
//tip vertex has weight -1/81
heh = _m.opposite_halfedge_handle(heh);
assert(heh.is_valid());
assert(_m.next_halfedge_handle(heh).is_valid());
assert(_m.to_vertex_handle(_m.next_halfedge_handle(heh)).is_valid());
pos -= real_t(1.0/81) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
//outer vertices have weight -2/81
heh = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh));
assert(heh.is_valid());
assert(_m.next_halfedge_handle(heh).is_valid());
assert(_m.to_vertex_handle(_m.next_halfedge_handle(heh)).is_valid());
pos -= real_t(2.0/81) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
heh = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh));
assert(heh.is_valid());
assert(_m.next_halfedge_handle(heh).is_valid());
assert(_m.to_vertex_handle(_m.next_halfedge_handle(heh)).is_valid());
pos -= real_t(2.0/81) * _m.point(_m.to_vertex_handle(_m.next_halfedge_handle(heh)));
}
}
else
{
//only use irregular vertices:
for(fheit = _m.fh_iter( fit ); fheit; ++fheit)
{
vh = _m.to_vertex_handle(fheit);
if( (_m.valence(vh) != 6) && (!_m.is_boundary(vh)) )
{
unsigned int K = _m.valence(vh);
pos += weights_[K][K]*_m.point(vh);
heh = _m.opposite_halfedge_handle( fheit.handle() );
for(unsigned int i = 0; i<K; ++i, heh = _m.opposite_halfedge_handle(_m.prev_halfedge_handle(heh)) )
{
pos += weights_[K][i]*_m.point(_m.to_vertex_handle(heh));
}
}
}
pos *= real_t(1.0/(3-nOrdinary));
}
vh = _m.add_vertex( pos );
_m.property(fp_pos_, fit.handle()) = vh;
}
}
//split faces
for (fit = _m.faces_begin();fit != fend; ++fit)
{
if ( _m.is_boundary(fit) && (gen%2))
{
boundary_split( _m, fit );
}
else
{
assert(_m.property(fp_pos_, fit.handle()).is_valid());
_m.split( fit, _m.property(fp_pos_, fit.handle()) );
}
}
// flip old edges
for (eit=_m.edges_begin(); eit != _m.edges_end(); ++eit)
if ( _m.status( eit ).tagged() && !_m.is_boundary( eit ) )
_m.flip(eit);
// Now we have an consistent mesh!
ASSERT_CONSISTENCY( MeshType, _m );
// increase generation by one
++gen;
}
return true;
}
private:
// Pre-compute location of new boundary points for odd generations
// and store them in the edge property ep_nv_;
void compute_new_boundary_points( MeshType& _m,
const typename MeshType::EdgeHandle& _eh)
{
assert( _m.is_boundary(_eh) );
typename MeshType::HalfedgeHandle heh;
typename MeshType::VertexHandle vh1, vh2, vh3, vh4, vhl, vhr;
typename MeshType::Point zero(0,0,0), P1, P2, P3, P4;
/*
// *---------*---------*
// / \ / \ / \
// / \ / \ / \
// / \ / \ / \
// / \ / \ / \
// *---------*--#---#--*---------*
//
// ^ ^ ^ ^ ^ ^
// P1 P2 pl pr P3 P4
*/
// get halfedge pointing from P3 to P2 (outer boundary halfedge)
heh = _m.halfedge_handle(_eh,
_m.is_boundary(_m.halfedge_handle(_eh,1)));
assert( _m.is_boundary( _m.next_halfedge_handle( heh ) ) );
assert( _m.is_boundary( _m.prev_halfedge_handle( heh ) ) );
vh1 = _m.to_vertex_handle( _m.next_halfedge_handle( heh ) );
vh2 = _m.to_vertex_handle( heh );
vh3 = _m.from_vertex_handle( heh );
vh4 = _m.from_vertex_handle( _m.prev_halfedge_handle( heh ));
P1 = _m.point(vh1);
P2 = _m.point(vh2);
P3 = _m.point(vh3);
P4 = _m.point(vh4);
vhl = _m.add_vertex(real_t(-5.0/81)*P1 + real_t(20.0/27)*P2 + real_t(10.0/27)*P3 + real_t(-4.0/81)*P4);
vhr = _m.add_vertex(real_t(-5.0/81)*P4 + real_t(20.0/27)*P3 + real_t(10.0/27)*P2 + real_t(-4.0/81)*P1);
_m.property(ep_nv_, _eh).first = vhl;
_m.property(ep_nv_, _eh).second = vhr;
}
void boundary_split( MeshType& _m, const typename MeshType::FaceHandle& _fh )
{
assert( _m.is_boundary(_fh) );
typename MeshType::VertexHandle vhl, vhr;
typename MeshType::FaceEdgeIter fe_it;
typename MeshType::HalfedgeHandle heh;
// find boundary edge
for( fe_it=_m.fe_iter( _fh ); fe_it && !_m.is_boundary( fe_it ); ++fe_it ) {};
// use precomputed, already inserted but not linked vertices
vhl = _m.property(ep_nv_, fe_it).first;
vhr = _m.property(ep_nv_, fe_it).second;
/*
// *---------*---------*
// / \ / \ / \
// / \ / \ / \
// / \ / \ / \
// / \ / \ / \
// *---------*--#---#--*---------*
//
// ^ ^ ^ ^ ^ ^
// P1 P2 pl pr P3 P4
*/
// get halfedge pointing from P2 to P3 (inner boundary halfedge)
heh = _m.halfedge_handle(fe_it,
_m.is_boundary(_m.halfedge_handle(fe_it,0)));
typename MeshType::HalfedgeHandle pl_P3;
// split P2->P3 (heh) in P2->pl (heh) and pl->P3
boundary_split( _m, heh, vhl ); // split edge
pl_P3 = _m.next_halfedge_handle( heh ); // store next halfedge handle
boundary_split( _m, heh ); // split face
// split pl->P3 in pl->pr and pr->P3
boundary_split( _m, pl_P3, vhr );
boundary_split( _m, pl_P3 );
assert( _m.is_boundary( vhl ) && _m.halfedge_handle(vhl).is_valid() );
assert( _m.is_boundary( vhr ) && _m.halfedge_handle(vhr).is_valid() );
}
void boundary_split(MeshType& _m,
const typename MeshType::HalfedgeHandle& _heh,
const typename MeshType::VertexHandle& _vh)
{
assert( _m.is_boundary( _m.edge_handle(_heh) ) );
typename MeshType::HalfedgeHandle
heh(_heh),
opp_heh( _m.opposite_halfedge_handle(_heh) ),
new_heh, opp_new_heh;
typename MeshType::VertexHandle to_vh(_m.to_vertex_handle(heh));
typename MeshType::HalfedgeHandle t_heh;
/*
* P5
* *
* /|\
* / \
* / \
* / \
* / \
* /_ heh new \
* *-----\*-----\*\-----*
* ^ ^ t_heh
* _vh to_vh
*
* P1 P2 P3 P4
*/
// Re-Setting Handles
// find halfedge point from P4 to P3
for(t_heh = heh;
_m.next_halfedge_handle(t_heh) != opp_heh;
t_heh = _m.opposite_halfedge_handle(_m.next_halfedge_handle(t_heh)))
{}
assert( _m.is_boundary( t_heh ) );
new_heh = _m.new_edge( _vh, to_vh );
opp_new_heh = _m.opposite_halfedge_handle(new_heh);
// update halfedge connectivity
_m.set_next_halfedge_handle(t_heh, opp_new_heh); // P4-P3 -> P3-P2
_m.set_next_halfedge_handle(new_heh, _m.next_halfedge_handle(heh)); // P2-P3 -> P3-P5
_m.set_next_halfedge_handle(heh, new_heh); // P1-P2 -> P2-P3
_m.set_next_halfedge_handle(opp_new_heh, opp_heh); // P3-P2 -> P2-P1
// both opposite halfedges point to same face
_m.set_face_handle(opp_new_heh, _m.face_handle(opp_heh));
// let heh finally point to new inserted vertex
_m.set_vertex_handle(heh, _vh);
// let heh and new_heh point to same face
_m.set_face_handle(new_heh, _m.face_handle(heh));
// let opp_new_heh be the new outgoing halfedge for to_vh
// (replaces for opp_heh)
_m.set_halfedge_handle( to_vh, opp_new_heh );
// let opp_heh be the outgoing halfedge for _vh
_m.set_halfedge_handle( _vh, opp_heh );
}
void boundary_split( MeshType& _m,
const typename MeshType::HalfedgeHandle& _heh)
{
assert( _m.is_boundary( _m.opposite_halfedge_handle( _heh ) ) );
typename MeshType::HalfedgeHandle
heh(_heh),
n_heh(_m.next_halfedge_handle(heh));
typename MeshType::VertexHandle
to_vh(_m.to_vertex_handle(heh));
typename MeshType::HalfedgeHandle
heh2(_m.new_edge(to_vh,
_m.to_vertex_handle(_m.next_halfedge_handle(n_heh)))),
heh3(_m.opposite_halfedge_handle(heh2));
typename MeshType::FaceHandle
new_fh(_m.new_face()),
fh(_m.face_handle(heh));
// Relink (half)edges
_m.set_face_handle(heh, new_fh);
_m.set_face_handle(heh2, new_fh);
_m.set_next_halfedge_handle(heh2, _m.next_halfedge_handle(_m.next_halfedge_handle(n_heh)));
_m.set_next_halfedge_handle(heh, heh2);
_m.set_face_handle( _m.next_halfedge_handle(heh2), new_fh);
_m.set_next_halfedge_handle(heh3, n_heh);
_m.set_next_halfedge_handle(_m.next_halfedge_handle(n_heh), heh3);
_m.set_face_handle(heh3, fh);
_m.set_halfedge_handle( fh, n_heh);
_m.set_halfedge_handle(new_fh, heh);
}
private:
weights_t weights_;
OpenMesh::FPropHandleT< typename MeshType::VertexHandle > fp_pos_;
OpenMesh::EPropHandleT< std::pair< typename MeshType::VertexHandle,
typename MeshType::VertexHandle> > ep_nv_;
OpenMesh::MPropHandleT< size_t > mp_gen_;
};
//=============================================================================
} // END_NS_UNIFORM
} // END_NS_SUBDIVIDER
} // END_NS_OPENMESH
//=============================================================================
#endif // OPENMESH_SUBDIVIDER_UNIFORM_SQRT3T_HH
//=============================================================================