IsotropicRemesherT_impl.hh

/*===========================================================================*\
*                                                                            *
*                              OpenFlipper                                   *
 *           Copyright (c) 2001-2015, RWTH-Aachen University                 *
 *           Department of Computer Graphics and Multimedia                  *
 *                          All rights reserved.                             *
 *                            www.openflipper.org                            *
 *                                                                           *
 *---------------------------------------------------------------------------*
 * This file is part of OpenFlipper.                                         *
 *---------------------------------------------------------------------------*
 *                                                                           *
 * Redistribution and use in source and binary forms, with or without        *
 * modification, are permitted provided that the following conditions        *
 * are met:                                                                  *
 *                                                                           *
 * 1. Redistributions of source code must retain the above copyright notice, *
 *    this list of conditions and the following disclaimer.                  *
 *                                                                           *
 * 2. Redistributions in binary form must reproduce the above copyright      *
 *    notice, this list of conditions and the following disclaimer in the    *
 *    documentation and/or other materials provided with the distribution.   *
 *                                                                           *
 * 3. Neither the name of the copyright holder nor the names of its          *
 *    contributors may be used to endorse or promote products derived from   *
 *    this software without specific prior written permission.               *
 *                                                                           *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS       *
 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A           *
 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,  *
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,       *
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR        *
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF    *
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING      *
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS        *
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.              *
*                                                                            *
\*===========================================================================*/
 
 
 
 
#define ISOTROPICREMESHER_C
 
#include "IsotropicRemesherT.hh"
 
#include <OpenMesh/Core/Mesh/PolyConnectivity.hh>
 
#include <ACG/Geometry/Algorithms.hh>
 
// -------------------- BSP
#include <ACG/Geometry/bsp/TriangleBSPT.hh>
 
template< class MeshT >
void IsotropicRemesher< MeshT >::remesh( MeshT& _mesh, const double _targetEdgeLength )
{
  const double low  = (4.0 / 5.0) * _targetEdgeLength;
  const double high = (4.0 / 3.0) * _targetEdgeLength;
 
  MeshT meshCopy = _mesh;
  OpenMeshTriangleBSPT< MeshT >* triangleBSP = getTriangleBSP(meshCopy);
 
  for (int i=0; i < 10; i++){
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 0);
// std::cerr << "Iteration = " << i << std::endl;
 
    splitLongEdges(_mesh, high);
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 2);
    
// std::cerr << "collapse" << std::endl;
    collapseShortEdges(_mesh, low, high);
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 4);
    
// std::cerr << "equal" << std::endl;
    equalizeValences(_mesh);
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 6);
    
// std::cerr << "relax" << std::endl;
    tangentialRelaxation(_mesh);
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 8);
    
// std::cerr << "project" << std::endl;
    projectToSurface(_mesh, meshCopy, triangleBSP);
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 10);
  }
 
  delete triangleBSP;
}
 
template< class MeshT >
OpenMeshTriangleBSPT< MeshT >* IsotropicRemesher< MeshT >::getTriangleBSP(MeshT& _mesh)
{
  // create Triangle BSP
  OpenMeshTriangleBSPT< MeshT >* triangle_bsp = new OpenMeshTriangleBSPT< MeshT >( _mesh );
 
  // build Triangle BSP
  triangle_bsp->reserve(_mesh.n_faces());
 
  for (auto f_it : _mesh.faces())
      triangle_bsp->push_back( f_it );
 
  triangle_bsp->build(10, 100);
 
  // return pointer to triangle bsp
  return triangle_bsp;
}
 
template< class MeshT >
void IsotropicRemesher< MeshT >::splitLongEdges( MeshT& _mesh, const double _maxEdgeLength )
{
 
  const double _maxEdgeLengthSqr = _maxEdgeLength * _maxEdgeLength;
 
  //iterate over all edges
  for (auto e_it : _mesh.edges()) {
    const typename OpenMesh::SmartHalfedgeHandle & hh = e_it.h0();
 
    const typename MeshT::VertexHandle & v0 = hh.from();
    const typename MeshT::VertexHandle & v1 = hh.to();
 
    typename MeshT::Point vec = _mesh.point(v1) - _mesh.point(v0);
 
    //edge to long?
    if ( vec.sqrnorm() > _maxEdgeLengthSqr ){
 
      const typename MeshT::Point midPoint = _mesh.point(v0) + ( 0.5 * vec );
 
      //split at midpoint
      typename MeshT::VertexHandle vh = _mesh.add_vertex( midPoint );
      
      bool hadFeature = _mesh.status(e_it).feature();
      
      _mesh.split(e_it, vh);
      
      if ( hadFeature ){
        
          for (auto voh_it : _mesh.voh_range(vh))
              if ( voh_it.to() == v0 || voh_it.to() == v1 )
                  _mesh.status( voh_it.edge() ).set_feature( true );
      }
    }
  }
}
 
template< class MeshT >
void IsotropicRemesher< MeshT >::collapseShortEdges( MeshT& _mesh, const double _minEdgeLength, const double _maxEdgeLength )
{
  const double _minEdgeLengthSqr = _minEdgeLength * _minEdgeLength;
  const double _maxEdgeLengthSqr = _maxEdgeLength * _maxEdgeLength;
 
  //add checked property
  OpenMesh::EPropHandleT< bool > checked;
  if ( !_mesh.get_property_handle(checked, "Checked Property") )
    _mesh.add_property(checked,"Checked Property" );
 
  //init property
  for (auto e_it : _mesh.edges())
    _mesh.property(checked, e_it) = false;
 
 
  bool finished = false;
 
  while( !finished ){
 
    finished = true;
 
    for (auto e_it : _mesh.edges()) {
    
        if ( _mesh.property(checked, e_it) )
          continue;
      
        _mesh.property(checked, e_it) = true;
 
        const typename OpenMesh::SmartHalfedgeHandle & hh = e_it.h0();
 
        const typename MeshT::VertexHandle & v0 = hh.from();
        const typename MeshT::VertexHandle & v1 = hh.to();
 
        const typename MeshT::Point vec = _mesh.point(v1) - _mesh.point(v0);
 
        const double edgeLength = vec.sqrnorm();
 
        // edge too short but don't try to collapse edges that have length 0
        if ( (edgeLength < _minEdgeLengthSqr) && (edgeLength > DBL_EPSILON) ){
    
          //check if the collapse is ok
          const typename MeshT::Point & B = _mesh.point(v1);
 
          bool collapse_ok = true;
    
          for (auto voh_it : _mesh.voh_range(v0))
              if ( (( B - _mesh.point( voh_it.to() ) ).sqrnorm() > _maxEdgeLengthSqr )
                   || ( _mesh.status( voh_it.edge() ).feature())
                   || ( voh_it.edge().is_boundary() ) ){
                collapse_ok = false;
                break;
              }
 
          if( collapse_ok && _mesh.is_collapse_ok(hh) ) {
            _mesh.collapse( hh );
 
            finished = false;
          }
 
        } 
    }
 
  }
 
  _mesh.remove_property(checked);
 
  _mesh.garbage_collection();
}
 
template< class MeshT >
void IsotropicRemesher< MeshT >::equalizeValences( MeshT& _mesh )
{
 
  for (auto e_it : _mesh.edges()) {
 
    if ( !_mesh.is_flip_ok(e_it) ) continue;
    if ( _mesh.status( e_it ).feature() ) continue;
 
    const typename OpenMesh::SmartHalfedgeHandle & h0 = e_it.h0();
    const typename OpenMesh::SmartHalfedgeHandle & h1 = e_it.h1();
 
    if (h0.is_valid() && h1.is_valid())
 
      if (h0.face().is_valid() && h1.face().is_valid()){
        //get vertices of corresponding faces
        //const typename MeshT::VertexHandle & a = h0.to();
        const typename OpenMesh::SmartVertexHandle & a = h0.to();
        const typename OpenMesh::SmartVertexHandle & b = h1.to();
        const typename OpenMesh::SmartVertexHandle & c = h0.next().to();
        const typename OpenMesh::SmartVertexHandle & d = h1.next().to();
 
        const int deviation_pre =  abs((int)(a.valence() - targetValence(_mesh, a)))
                                  +abs((int)(b.valence() - targetValence(_mesh, b)))
                                  +abs((int)(c.valence() - targetValence(_mesh, c)))
                                  +abs((int)(d.valence() - targetValence(_mesh, d)));
        _mesh.flip(e_it);
 
        const int deviation_post = abs((int)(a.valence() - targetValence(_mesh, a)))
                                  +abs((int)(b.valence() - targetValence(_mesh, b)))
                                  +abs((int)(c.valence() - targetValence(_mesh, c)))
                                  +abs((int)(d.valence() - targetValence(_mesh, d)));
 
        if (deviation_pre <= deviation_post)
          _mesh.flip(e_it);
      }
  }
}
 
template< class MeshT > 
inline
int IsotropicRemesher< MeshT >::targetValence( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
 
  if (isBoundary(_mesh,_vh))
    return 4;
  else
    return 6;
}
 
template< class MeshT > 
inline
bool IsotropicRemesher< MeshT >::isBoundary( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
 
  for (auto voh_it : _mesh.voh_range(_vh))
    if ( voh_it.edge().is_boundary() )
      return true;
 
  return false;
}
 
template< class MeshT > 
inline
bool IsotropicRemesher< MeshT >::isFeature( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
 
  for (auto voh_it : _mesh.voh_range(_vh))
    if ( _mesh.status( voh_it.edge()).feature() )
      return true;
 
  return false;
}
 
template< class MeshT >
void IsotropicRemesher< MeshT >::tangentialRelaxation( MeshT& _mesh )
{
  _mesh.update_normals();
 
  //add checked property
  OpenMesh::VPropHandleT< typename MeshT::Point > q;
  if ( !_mesh.get_property_handle(q, "q Property") )
    _mesh.add_property(q,"q Property" );
 
  //first compute barycenters
  for (auto v_it : _mesh.vertices()){
 
    typename MeshT::Point tmp(0.0, 0.0, 0.0);
    uint N = 0;
 
    for (auto vv_it : v_it.vertices()){
      tmp += _mesh.point(vv_it);
      N++;
    }
 
    if (N > 0)
      tmp /= (double) N;
 
    _mesh.property(q, v_it) = tmp;
  }
 
  //move to new position
  for (auto v_it : _mesh.vertices()){
    if ( !v_it.is_boundary() && !v_it.feature() )
      _mesh.set_point(v_it,  _mesh.property(q, v_it) + (_mesh.normal(v_it)| (_mesh.point(v_it) - _mesh.property(q, v_it) ) ) * _mesh.normal(v_it));
  }
 
  _mesh.remove_property(q);
}
 
template <class MeshT>
template <class SpatialSearchT>
typename MeshT::Point
IsotropicRemesher< MeshT >::findNearestPoint(const MeshT&                   _mesh,
                                          const typename MeshT::Point&   _point,
                                          typename MeshT::FaceHandle&    _fh,
                                          SpatialSearchT*                _ssearch,
                                          double*                        _dbest)
{
 
  typename MeshT::Point  p_best = _mesh.point(_mesh.vertex_handle(0));
  typename MeshT::Scalar d_best = (_point-p_best).sqrnorm();
 
  typename MeshT::FaceHandle fh_best;
 
  if( _ssearch == 0 )
  {
    // exhaustive search
    for (auto cf_it : _mesh.faces()) {
      typename MeshT::ConstFaceVertexIter cfv_it = _mesh.cfv_iter(cf_it);
 
      const typename MeshT::Point& pt0 = _mesh.point(   *cfv_it);
      const typename MeshT::Point& pt1 = _mesh.point( *(++cfv_it));
      const typename MeshT::Point& pt2 = _mesh.point( *(++cfv_it));
 
      typename MeshT::Point ptn;
 
      typename MeshT::Scalar d = ACG::Geometry::distPointTriangleSquared( _point,
                     pt0,
                     pt1,
                     pt2,
                     ptn );
 
      if( d < d_best && d >= 0.0)
      {
        d_best = d;
        p_best = ptn;
 
        fh_best = cf_it;
      }
    }
 
    // return facehandle
    _fh = fh_best;
 
    // return distance
    if(_dbest)
      *_dbest = sqrt(d_best);
 
 
    return p_best;
  }
  else
  {
    typename MeshT::FaceHandle     fh = _ssearch->nearest(_point).handle;
    typename MeshT::CFVIter        fv_it = _mesh.cfv_iter(fh);
 
    const typename MeshT::Point&   pt0 = _mesh.point( *(  fv_it));
    const typename MeshT::Point&   pt1 = _mesh.point( *(++fv_it));
    const typename MeshT::Point&   pt2 = _mesh.point( *(++fv_it));
 
    // project
    d_best = ACG::Geometry::distPointTriangleSquared(_point, pt0, pt1, pt2, p_best);
 
    // return facehandle
    _fh = fh;
 
    // return distance
    if(_dbest)
      *_dbest = sqrt(d_best);
 
    return p_best;
  }
}
 
 
template< class MeshT >
template< class SpatialSearchT >
void IsotropicRemesher< MeshT >::projectToSurface( MeshT& _mesh, MeshT& _original, SpatialSearchT* _ssearch )
{
 
  for (auto v_it : _mesh.vertices()){
 
    if ( v_it.is_boundary() ) continue;
    if ( v_it.feature() ) continue;
 
    typename MeshT::Point p = _mesh.point(v_it);
    typename MeshT::FaceHandle fhNear;
    double distance;
 
    typename MeshT::Point pNear = findNearestPoint(_original, p, fhNear, _ssearch, &distance);
 
    _mesh.set_point(v_it, pNear);
  }
}