Commit ce620e8e authored by Jan Möbius's avatar Jan Möbius

git-svn-id: http://www.openflipper.org/svnrepo/OpenFlipper/branches/Free@5085 383ad7c9-94d9-4d36-a494-682f7c89f535
parent 1ca729c4
#== SYSTEM PART -- DON'T TOUCH ==============================================
include $(ACGMAKE)/Config
#==============================================================================
SUBDIRS :=
PACKAGES := openmesh2 acg2
PROJ_LIBS :=
MODULES := cxx
#== SYSTEM PART -- DON'T TOUCH ==============================================
include $(ACGMAKE)/Rules
#==============================================================================
//=============================================================================
//
// OpenFlipper
// Copyright (C) 2008 by Computer Graphics Group, RWTH Aachen
// www.openflipper.org
//
//-----------------------------------------------------------------------------
//
// License
//
// OpenFlipper 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.
//
// OpenFlipper 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 Lesser General Public License
// along with OpenFlipper. If not, see <http://www.gnu.org/licenses/>.
//
//-----------------------------------------------------------------------------
//
// $Revision$
// $Author$
// $Date$
//
//=============================================================================
//=============================================================================
//
// IMPLEMENTATION
//
//=============================================================================
#define CURVATURE_C
//== INCLUDES =================================================================
// #include "MeshFunctions.hh"
// #include <ACG/Geometry/Algorithms.hh>
// #include "Math_Tools/GeometryFunctions.hh"
#include "Math_Tools/Math_Tools.hh"
#include "Math_Tools/GeometryFunctions.hh"
#include <iostream>
#include <OpenMesh/Core/Geometry/MathDefs.hh>
#include <math.h>
//== NAMESPACES ===============================================================
namespace curvature {
//== IMPLEMENTATION ==========================================================
/*! compute consistent dirscrete gaussian curvature (vertex is a small sphere patch, edges are small cylinders)
*/
template< typename MeshT >
double
gauss_curvature(MeshT& _mesh, const typename MeshT::VertexHandle& _vh) {
if (_mesh.status(_vh).deleted())
return 0.0;
double gauss_curv( 2.0*M_PI);
const typename MeshT::Point p0 = _mesh.point(_vh);
typename MeshT::CVOHIter voh_it( _mesh.cvoh_iter(_vh));
typename MeshT::CVOHIter n_voh_it = voh_it;
if ( ! voh_it.handle().is_valid() )
return 0.0;
// move to next
++n_voh_it;
for(; voh_it; ++voh_it, ++n_voh_it)
{
typename MeshT::Point p1 = _mesh.point(_mesh.to_vertex_handle( voh_it));
typename MeshT::Point p2 = _mesh.point(_mesh.to_vertex_handle( n_voh_it));
gauss_curv -= acos(OpenMesh::sane_aarg( ((p1-p0).normalize() | (p2-p0).normalize()) ));
}
return gauss_curv;
}
template<class MeshT, class VectorT, class REALT>
void discrete_mean_curv_op( const MeshT& _m,
const typename MeshT::VertexHandle& _vh,
VectorT& _n,
REALT& _area )
{
_n = VectorT(0,0,0);
_area = 0.0;
typename MeshT::ConstVertexOHalfedgeIter voh_it = _m.cvoh_iter(_vh);
if ( ! voh_it.handle().is_valid() )
return;
for(; voh_it; ++voh_it)
{
if ( _m.is_boundary( _m.edge_handle( voh_it.handle() ) ) )
continue;
const typename MeshT::Point p0 = _m.point( _vh );
const typename MeshT::Point p1 = _m.point( _m.to_vertex_handle( voh_it.handle()));
// const typename MeshT::Point p2 = _m.point( _m.to_vertex_handle( _m.next_halfedge_handle( voh_it.handle())));
const typename MeshT::Point p2 = _m.point( _m.from_vertex_handle( _m.prev_halfedge_handle( voh_it.handle())));
const typename MeshT::Point p3 = _m.point( _m.to_vertex_handle( _m.next_halfedge_handle( _m.opposite_halfedge_handle(voh_it.handle()))));
const REALT alpha = acos( OpenMesh::sane_aarg((p0-p2).normalize() | (p1-p2).normalize()) );
const REALT beta = acos( OpenMesh::sane_aarg((p0-p3).normalize() | (p1-p3).normalize()) );
REALT cotw = 0.0;
if ( !OpenMesh::is_eq(alpha,M_PI/2.0) )
cotw += (REALT(1.0))/tan(alpha);
if ( !OpenMesh::is_eq(beta,M_PI/2.0) )
cotw += (REALT(1.0))/tan(beta);
#ifdef WIN32
if ( OpenMesh::is_zero(cotw) )
continue;
#else
if ( OpenMesh::is_zero(cotw) || isinf(cotw) )
continue;
#endif
// calculate area
const int obt = GeometryFunctions::is_obtuse(p0,p1,p2);
if(obt == 0)
{
REALT gamma = acos( OpenMesh::sane_aarg((p0-p1).normalize() | (p2-p1).normalize()) );
REALT tmp = 0.0;
if ( !OpenMesh::is_eq(alpha,M_PI/2.0) )
tmp += (p0-p1).sqrnorm()*1.0/tan(alpha);
if ( !OpenMesh::is_eq(gamma,M_PI/2.0) )
tmp += (p0-p2).sqrnorm()*1.0/tan(gamma);
#ifdef WIN32
if ( OpenMesh::is_zero(tmp) )
continue;
#else
if ( OpenMesh::is_zero(tmp) || isinf(tmp) )
continue;
#endif
_area += 0.125*( tmp );
}
else
{
if(obt == 1)
{
_area += ((p1-p0) % (p2-p0)).norm() * 0.5 * 0.5;
}
else
{
_area += ((p1-p0) % (p2-p0)).norm() * 0.5 * 0.25;
}
}
_n += ((p0-p1)*cotw);
// error handling
//if(_area < 0) std::cerr << "error: triangle area < 0\n";
// if(isnan(_area))
// {
// REALT gamma = acos( ((p0-p1).normalize() | (p2-p1).normalize()) );
/*
std::cerr << "***************************\n";
std::cerr << "error : trianlge area = nan\n";
std::cerr << "alpha : " << alpha << std::endl;
std::cerr << "beta : " << beta << std::endl;
std::cerr << "gamma : " << gamma << std::endl;
std::cerr << "cotw : " << cotw << std::endl;
std::cerr << "normal: " << _n << std::endl;
std::cerr << "p0 : " << p0 << std::endl;
std::cerr << "p1 : " << p1 << std::endl;
std::cerr << "p2 : " << p2 << std::endl;
std::cerr << "p3 : " << p3 << std::endl;
std::cerr << "***************************\n";
*/
// }
}
_n /= 4.0*_area;
}
//=============================================================================
} // curvature Namespace
//=============================================================================
//=============================================================================
//
// OpenFlipper
// Copyright (C) 2008 by Computer Graphics Group, RWTH Aachen
// www.openflipper.org
//
//-----------------------------------------------------------------------------
//
// License
//
// OpenFlipper 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.
//
// OpenFlipper 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 Lesser General Public License
// along with OpenFlipper. If not, see <http://www.gnu.org/licenses/>.
//
//-----------------------------------------------------------------------------
//
// $Revision$
// $Author$
// $Date$
//
//=============================================================================
//=============================================================================
//
//
//=============================================================================
#ifndef CURVATURE_HH
#define CURVATURE_HH
/*! \file Curvature.hh
\brief Functions for calculating curvatures
*/
//== INCLUDES =================================================================
#include <vector>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
namespace curvature {
//== DEFINITIONS ==============================================================
/*! compute consistent dirscrete gaussian curvature (vertex is a small sphere patch, edges are small cylinders)
*/
template< typename MeshT >
double
gauss_curvature(MeshT& _mesh, const typename MeshT::VertexHandle& _vh);
/**
* Mean Curvature Normal Operator
* warning: if mean curvature < 0 _n points to the inside
* warning: if mean curvature = 0 -> no normal information
@param _n mean_curvature(vit)*n(vit)
@param _area global vertex area
*/
template<class MeshT, class VectorT, class REALT>
void discrete_mean_curv_op( const MeshT& _m,
const typename MeshT::VertexHandle& _vh,
VectorT& _n,
REALT& _area );
//=============================================================================
}
//=============================================================================
#if defined(INCLUDE_TEMPLATES) && !defined(CURVATURE_C)
#define CURVATURE_TEMPLATES
#include "Curvature.cc"
#endif
//=============================================================================
#endif // CURVATURE_HH defined
//=============================================================================
//=============================================================================
//
// OpenFlipper
// Copyright (C) 2008 by Computer Graphics Group, RWTH Aachen
// www.openflipper.org
//
//-----------------------------------------------------------------------------
//
// License
//
// OpenFlipper 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.
//
// OpenFlipper 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 Lesser General Public License
// along with OpenFlipper. If not, see <http://www.gnu.org/licenses/>.
//
//-----------------------------------------------------------------------------
//
// $Revision$
// $Author$
// $Date$
//
//=============================================================================
//=============================================================================
//
// IMPLEMENTATION
//
//=============================================================================
#define MESHFUNCTIONS_C
//== INCLUDES =================================================================
#include "MeshFunctions.hh"
#include <ACG/Geometry/Algorithms.hh>
#include "Math_Tools/GeometryFunctions.hh"
#include "Math_Tools/Math_Tools.hh"
#include <set>
#include <iostream>
#include <OpenMesh/Core/Geometry/MathDefs.hh>
//== NAMESPACES ===============================================================
namespace MeshFunctions {
//== IMPLEMENTATION ==========================================================
template < typename MeshT , typename VectorT >
bool get_boundary(MeshT& _mesh,
typename MeshT::VertexHandle _vh,
std::vector< std::pair< VectorT , typename MeshT::VertexHandle > >& _boundary)
{
_boundary.clear();
typename MeshT::VertexHandle last = _vh;
const typename MeshT::VertexHandle begin = _vh;
std::set< typename MeshT::VertexHandle > set;
//walk along boundary
do {
//insert current Vertex (which is on boundary to set (needed to detect loops)
set.insert( _vh );
// store the current vertex in the boundary list
_boundary.push_back(std::pair< VectorT , typename MeshT::VertexHandle > ( (VectorT)_mesh.point(_vh) , _vh ) );
// iterate over all outgoing halfedges of the current vertex to find next one
for (typename MeshT::VertexOHalfedgeIter vohe_it(_mesh,_vh); vohe_it ; ++vohe_it) {
//if vertex is on the boundary use it as the next one (if its not the one, we are comming from)
if ( _mesh.is_boundary(vohe_it) && ( _mesh.to_vertex_handle(vohe_it) != last ) ) {
last = _vh;
_vh = _mesh.to_vertex_handle(vohe_it);
break;
}
}
// stop if we find a vertex which is already in the list (can also detect partial loops)
} while ( set.count( _vh ) == 0 );
if ( begin != _vh ) {
std::cout << "Warning in ( GeometryFunctions.cc get_boundary ) : boundary loop may be partial ( start- != endpoint ) " << std::endl;
}
return true;
}
template < typename MeshT , typename VectorT >
bool get_boundary(MeshT& _mesh,
std::vector< std::pair< VectorT , typename MeshT::VertexHandle > >& _boundary)
{
typename MeshT::VertexHandle vh;
// Search for one Vertex on boundary
bool found = false;
for (typename MeshT::HalfedgeIter he_it=_mesh.halfedges_begin() ; he_it != _mesh.halfedges_end() ; ++he_it) {
if ( _mesh.is_boundary(he_it) ) {
vh = _mesh.from_vertex_handle(he_it);
found = true;
break;
}
}
if ( found )
return get_boundary(_mesh , vh , _boundary);
else {
std::cerr << "Did not find Mesh boundary!! ( GeometryFunctions.cc get_boundary )" << std::endl;
return false;
}
}
template < typename MeshT , typename VectorT >
void smooth_boundary(MeshT& _mesh,
typename MeshT::VertexHandle _vh)
{
std::vector< std::pair< VectorT , typename MeshT::VertexHandle > > boundary;
//get the boundary
get_boundary (_mesh , _vh , boundary);
for (uint i = 1 ; i <= boundary.size() ; ++i ) {
_mesh.point( boundary[ i % boundary.size() ].second ) = ( boundary[ ( i - 1) % boundary.size() ].first
+ boundary[ ( i ) % boundary.size() ].first * 2.0
+ boundary[ ( i + 1) % boundary.size() ].first ) * 0.25;
}
}
template < typename MeshT >
bool neighbour(const MeshT& _mesh ,
const typename MeshT::FaceHandle _fh1 ,
const typename MeshT::FaceHandle _fh2 )
{
for ( typename MeshT::FaceFaceIter ff_it(_mesh,_fh1) ; ff_it ; ++ff_it)
if (ff_it.handle() == _fh2)
return true;
return false;
}
template <class MeshT , typename VectorT >
bool
cut_face(const VectorT& _porigin,
const VectorT& _pnormal,
const MeshT& _mesh,
const typename MeshT::FaceHandle& _fh)
{
typename MeshT::ConstFaceVertexIter fv_it = _mesh.cfv_iter(_fh);
const VectorT p0 = _mesh.point( fv_it);
const VectorT p1 = _mesh.point(++fv_it);
const VectorT p2 = _mesh.point(++fv_it);
unsigned int npos(0), nneg(0);
if( GeometryFunctions::dist_plane< VectorT , double >( _porigin, _pnormal, p0) > 0 ) ++npos; else ++nneg;
if( GeometryFunctions::dist_plane< VectorT , double >( _porigin, _pnormal, p1) > 0 ) ++npos; else ++nneg;
if( GeometryFunctions::dist_plane< VectorT , double >( _porigin, _pnormal, p2) > 0 ) ++npos; else ++nneg;
if( npos && nneg ) return true;
else return false;
}
template < typename MeshT >
double
calc_area( const MeshT& _mesh)
{
double area = 0.0;
for ( typename MeshT::ConstFaceIter f_it = _mesh.faces_begin() ; f_it != _mesh.faces_end() ; ++f_it) {
typename MeshT::ConstFaceVertexIter fv_it = _mesh.cfv_iter(f_it);
const ACG::Vec3d vertex1 = (ACG::Vec3d)_mesh.point( fv_it );
const ACG::Vec3d vertex2 = (ACG::Vec3d)_mesh.point( ++fv_it );
const ACG::Vec3d vertex3 = (ACG::Vec3d)_mesh.point( ++fv_it );
area += ((vertex1 - vertex2) % (vertex3-vertex2)).norm();
}
return area;
}
template < typename MeshT >
double
calc_angle_around( const MeshT& _mesh , const typename MeshT::VertexHandle _vh)
{
double angle = 0.0;
const typename MeshT::Point p0 = _mesh.point(_vh);
typename MeshT::ConstVertexOHalfedgeIter voh_it(_mesh,_vh);
typename MeshT::ConstVertexOHalfedgeIter nx_voh_it = voh_it;
++nx_voh_it;
for ( ; voh_it; ++voh_it , ++nx_voh_it) {
const typename MeshT::Point edge_1 = MathTools::sane_normalized( _mesh.point(_mesh.to_vertex_handle(voh_it)) - p0);
const typename MeshT::Point edge_2 = MathTools::sane_normalized( (_mesh.point(_mesh.to_vertex_handle(nx_voh_it))) - p0);
angle += acos(OpenMesh::sane_aarg(edge_1 | edge_2));
}
return angle;
}
//=============================================================================
} // MeshFunctions Namespace
//=============================================================================
//=============================================================================
//
// OpenFlipper
// Copyright (C) 2008 by Computer Graphics Group, RWTH Aachen
// www.openflipper.org
//
//-----------------------------------------------------------------------------
//
// License
//
// OpenFlipper 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.
//
// OpenFlipper 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 Lesser General Public License
// along with OpenFlipper. If not, see <http://www.gnu.org/licenses/>.
//
//-----------------------------------------------------------------------------
//
// $Revision$
// $Author$
// $Date$
//
//=============================================================================
//=============================================================================
//
//
//=============================================================================
#ifndef MESHFUNCTIONS_HH
#define MESHFUNCTIONS_HH
/*! \file MeshFunctions.hh
\brief Functions for modifying a Mesh
General file with template functions doing modifications on a given Mesh
(e.g smooth, reposition,... )
*/
//== INCLUDES =================================================================
#include <vector>
//== FORWARDDECLARATIONS ======================================================
//== NAMESPACES ===============================================================
/// Namespace providing different Mesh editing functions
namespace MeshFunctions {
//== DEFINITIONS ==============================================================
/**
get one boundary at the vertex
@param _mesh Mesh
@param _vh Vertex handle of one vertex of the boundary
@param _boundary Coords and vertex handles of the boundary
*/
template < typename MeshT , typename VectorT >
bool get_boundary(MeshT& _mesh,
typename MeshT::VertexHandle _vh,
std::vector< std::pair< VectorT , typename MeshT::VertexHandle > >& _boundary);
/**
get boundary of a mesh (assumes that there is only one boundary!!
@param _mesh Mesh
@param _boundary Coords and vertex handles of the boundary
@return Found a boundary?
*/
template < typename MeshT , typename VectorT >
bool get_boundary(MeshT& _mesh,
std::vector< std::pair< VectorT ,
typename MeshT::VertexHandle > >& _boundary);
/**
function to smooth a boundary of a given mesh by moving each vertex to the
mean Position of its adjazent verticies
@param _mesh Mesh to work on
@param _vh Vertex handle on the boundary
*/
template < typename MeshT , typename VectorT >
void smooth_boundary(MeshT& _mesh ,
typename MeshT::VertexHandle _vh);
/**
Checks for two faces if they are adjazent
@param _mesh Mesh
@param _fh1 First Face
@param _fh2 Second Face
*/
template < typename MeshT >
bool neighbour(const MeshT& _mesh ,
const typename MeshT::FaceHandle _fh1 ,