IsotropicRemesherT.cc 15.5 KB
Newer Older
Jan Möbius's avatar
Jan Möbius committed
1 2 3
/*===========================================================================*\
*                                                                            *
*                              OpenFlipper                                   *
4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38
 *           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.              *
Jan Möbius's avatar
Jan Möbius committed
39 40 41
*                                                                            *
\*===========================================================================*/

Jan Möbius's avatar
Jan Möbius committed
42

Jan Möbius's avatar
Jan Möbius committed
43

44 45 46 47 48 49 50 51

#define ISOTROPICREMESHER_C

#include "IsotropicRemesherT.hh"

#include <ACG/Geometry/Algorithms.hh>

// -------------------- BSP
52
#include <ACG/Geometry/bsp/TriangleBSPT.hh>
53 54 55 56 57 58 59 60 61 62 63 64

/// do the remeshing
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++){
65 66
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 0);
67
// std::cerr << "Iteration = " << i << std::endl;
68 69

    splitLongEdges(_mesh, high);
70 71 72
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 2);
    
73 74
// std::cerr << "collapse" << std::endl;
    collapseShortEdges(_mesh, low, high);
75 76 77
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 4);
    
78 79
// std::cerr << "equal" << std::endl;
    equalizeValences(_mesh);
80 81 82
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 6);
    
83 84
// std::cerr << "relax" << std::endl;
    tangentialRelaxation(_mesh);
85 86 87
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 8);
    
88 89
// std::cerr << "project" << std::endl;
    projectToSurface(_mesh, meshCopy, triangleBSP);
90 91
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 10);
92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109
  }

  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());

  typename MeshT::FIter f_it  = _mesh.faces_begin();
  typename MeshT::FIter f_end = _mesh.faces_end();

  for (; f_it!=f_end; ++f_it)
Jan Möbius's avatar
Jan Möbius committed
110
    triangle_bsp->push_back( *f_it );
111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129

  triangle_bsp->build(10, 100);

  // return pointer to triangle bsp
  return triangle_bsp;
}

/// performs edge splits until all edges are shorter than the threshold
template< class MeshT >
void IsotropicRemesher< MeshT >::splitLongEdges( MeshT& _mesh, const double _maxEdgeLength )
{

  const double _maxEdgeLengthSqr = _maxEdgeLength * _maxEdgeLength;

  typename MeshT::EdgeIter e_it;
  typename MeshT::EdgeIter e_end = _mesh.edges_end();

  //iterate over all edges
  for (e_it = _mesh.edges_begin(); e_it != e_end; ++e_it){
Jan Möbius's avatar
Jan Möbius committed
130
    const typename MeshT::HalfedgeHandle & hh = _mesh.halfedge_handle( *e_it, 0 );
131 132 133 134 135 136 137 138 139 140 141 142

    const typename MeshT::VertexHandle & v0 = _mesh.from_vertex_handle(hh);
    const typename MeshT::VertexHandle & v1 = _mesh.to_vertex_handle(hh);

    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
143 144
      typename MeshT::VertexHandle vh = _mesh.add_vertex( midPoint );
      
Jan Möbius's avatar
Jan Möbius committed
145
      bool hadFeature = _mesh.status(*e_it).feature();
146
      
Jan Möbius's avatar
Jan Möbius committed
147
      _mesh.split(*e_it, vh);
148 149 150 151
      
      if ( hadFeature ){
        
        typename MeshT::VOHIter vh_it;
Jan Möbius's avatar
Jan Möbius committed
152 153 154
        for (vh_it = _mesh.voh_iter(vh); vh_it.is_valid(); ++vh_it)
          if ( _mesh.to_vertex_handle(*vh_it) == v0 || _mesh.to_vertex_handle(*vh_it) == v1 )
            _mesh.status( _mesh.edge_handle( *vh_it ) ).set_feature( true );
155
      }
156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176
    }
  }
}

/// collapse edges shorter than minEdgeLength if collapsing doesn't result in new edge longer than maxEdgeLength
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
  typename MeshT::ConstEdgeIter e_it;
  typename MeshT::ConstEdgeIter e_end = _mesh.edges_end();

  for (e_it = _mesh.edges_begin(); e_it != e_end; ++e_it)
Jan Möbius's avatar
Jan Möbius committed
177
    _mesh.property(checked, *e_it) = false;
178 179 180 181 182 183 184 185 186 187


  bool finished = false;

  while( !finished ){

    finished = true;

    for (e_it = _mesh.edges_begin(); e_it != _mesh.edges_end() ; ++e_it){
    
Jan Möbius's avatar
Jan Möbius committed
188
        if ( _mesh.property(checked, *e_it) )
189 190
          continue;
      
Jan Möbius's avatar
Jan Möbius committed
191
        _mesh.property(checked, *e_it) = true;
192

Jan Möbius's avatar
Jan Möbius committed
193
        const typename MeshT::HalfedgeHandle & hh = _mesh.halfedge_handle( *e_it, 0 );
194 195 196
      
        const typename MeshT::VertexHandle & v0 = _mesh.from_vertex_handle(hh);
        const typename MeshT::VertexHandle & v1 = _mesh.to_vertex_handle(hh);
197

198
        const typename MeshT::Point vec = _mesh.point(v1) - _mesh.point(v0);
199 200 201

        const double edgeLength = vec.sqrnorm();

202 203
        // edge too short but don't try to collapse edges that have length 0
        if ( (edgeLength < _minEdgeLengthSqr) && (edgeLength > DBL_EPSILON) ){
204 205 206
    
          //check if the collapse is ok
          const typename MeshT::Point & B = _mesh.point(v1);
207

208 209
          bool collapse_ok = true;
    
Jan Möbius's avatar
Jan Möbius committed
210 211 212 213
          for(typename MeshT::VOHIter vh_it = _mesh.voh_iter(v0); vh_it.is_valid(); ++vh_it)
            if ( (( B - _mesh.point( _mesh.to_vertex_handle(*vh_it) ) ).sqrnorm() > _maxEdgeLengthSqr )
                 || ( _mesh.status( _mesh.edge_handle( *vh_it ) ).feature())
                 || ( _mesh.is_boundary( _mesh.edge_handle( *vh_it ) ) )){
214 215 216 217 218 219 220 221 222
              collapse_ok = false;
              break;
            }
  
          if( collapse_ok && _mesh.is_collapse_ok(hh) ) {
            _mesh.collapse( hh );

            finished = false;
          }
223

224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242
        } 
    }

  }

  _mesh.remove_property(checked);

  _mesh.garbage_collection();
}

template< class MeshT >
void IsotropicRemesher< MeshT >::equalizeValences( MeshT& _mesh )
{

  typename MeshT::EdgeIter e_it;
  typename MeshT::EdgeIter e_end = _mesh.edges_end();

  for (e_it = _mesh.edges_sbegin(); e_it != e_end; ++e_it){

Jan Möbius's avatar
Jan Möbius committed
243 244
    if ( !_mesh.is_flip_ok(*e_it) ) continue;
    if ( _mesh.status( *e_it ).feature() ) continue;
245
    
246

Jan Möbius's avatar
Jan Möbius committed
247 248
    const typename MeshT::HalfedgeHandle & h0 = _mesh.halfedge_handle( *e_it, 0 );
    const typename MeshT::HalfedgeHandle & h1 = _mesh.halfedge_handle( *e_it, 1 );
249 250 251 252 253 254 255 256 257 258 259 260 261 262

    if (h0.is_valid() && h1.is_valid())

      if (_mesh.face_handle(h0).is_valid() && _mesh.face_handle(h1).is_valid()){
        //get vertices of corresponding faces
        const typename MeshT::VertexHandle & a = _mesh.to_vertex_handle(h0);
        const typename MeshT::VertexHandle & b = _mesh.to_vertex_handle(h1);
        const typename MeshT::VertexHandle & c = _mesh.to_vertex_handle(_mesh.next_halfedge_handle(h0));
        const typename MeshT::VertexHandle & d = _mesh.to_vertex_handle(_mesh.next_halfedge_handle(h1));

        const int deviation_pre =  abs((int)(_mesh.valence(a) - targetValence(_mesh, a)))
                                  +abs((int)(_mesh.valence(b) - targetValence(_mesh, b)))
                                  +abs((int)(_mesh.valence(c) - targetValence(_mesh, c)))
                                  +abs((int)(_mesh.valence(d) - targetValence(_mesh, d)));
Jan Möbius's avatar
Jan Möbius committed
263
        _mesh.flip(*e_it);
264 265 266 267 268 269 270

        const int deviation_post = abs((int)(_mesh.valence(a) - targetValence(_mesh, a)))
                                  +abs((int)(_mesh.valence(b) - targetValence(_mesh, b)))
                                  +abs((int)(_mesh.valence(c) - targetValence(_mesh, c)))
                                  +abs((int)(_mesh.valence(d) - targetValence(_mesh, d)));

        if (deviation_pre <= deviation_post)
Jan Möbius's avatar
Jan Möbius committed
271
          _mesh.flip(*e_it);
272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
      }
  }
}

///returns 4 for boundary vertices and 6 otherwise
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 ){

  typename MeshT::ConstVertexOHalfedgeIter voh_it;

Jan Möbius's avatar
Jan Möbius committed
293 294
  for (voh_it = _mesh.voh_iter( _vh ); voh_it.is_valid(); ++voh_it )
    if ( _mesh.is_boundary( _mesh.edge_handle( *voh_it ) ) )
295 296 297 298 299
      return true;

  return false;
}

300 301 302 303 304 305
template< class MeshT > 
inline
bool IsotropicRemesher< MeshT >::isFeature( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){

  typename MeshT::ConstVertexOHalfedgeIter voh_it;

Jan Möbius's avatar
Jan Möbius committed
306 307
  for (voh_it = _mesh.voh_iter( _vh ); voh_it.is_valid(); ++voh_it )
    if ( _mesh.status( _mesh.edge_handle( *voh_it ) ).feature() )
308 309 310 311 312
      return true;

  return false;
}

313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332
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" );

  typename MeshT::ConstVertexIter v_it;
  typename MeshT::ConstVertexIter v_end = _mesh.vertices_end();

  //first compute barycenters
  for (v_it = _mesh.vertices_sbegin(); v_it != v_end; ++v_it){

    typename MeshT::Point tmp(0.0, 0.0, 0.0);
    uint N = 0;

    typename MeshT::VVIter vv_it;
Jan Möbius's avatar
Jan Möbius committed
333 334
    for (vv_it = _mesh.vv_iter(*v_it); vv_it.is_valid(); ++vv_it){
      tmp += _mesh.point(*vv_it);
335 336 337 338 339 340
      N++;
    }

    if (N > 0)
      tmp /= (double) N;

Jan Möbius's avatar
Jan Möbius committed
341
    _mesh.property(q, *v_it) = tmp;
342 343 344 345
  }

  //move to new position
  for (v_it = _mesh.vertices_sbegin(); v_it != v_end; ++v_it){
Jan Möbius's avatar
Jan Möbius committed
346 347
    if ( !isBoundary(_mesh, *v_it) && !isFeature(_mesh, *v_it) )
      _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));
348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375
  }

  _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
    typename MeshT::ConstFaceIter cf_it  = _mesh.faces_begin();
    typename MeshT::ConstFaceIter cf_end = _mesh.faces_end();

    for(; cf_it != cf_end; ++cf_it)
    {
Jan Möbius's avatar
Jan Möbius committed
376
      typename MeshT::ConstFaceVertexIter cfv_it = _mesh.cfv_iter(*cf_it);
377

Jan Möbius's avatar
Jan Möbius committed
378 379 380
      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));
381 382 383 384 385 386 387 388 389

      typename MeshT::Point ptn;

      typename MeshT::Scalar d = ACG::Geometry::distPointTriangleSquared( _point,
                     pt0,
                     pt1,
                     pt2,
                     ptn );

Anne Kathrein's avatar
Anne Kathrein committed
390
      if( d < d_best && d >= 0.0)
391 392 393 394
      {
        d_best = d;
        p_best = ptn;

Jan Möbius's avatar
Jan Möbius committed
395
        fh_best = *cf_it;
396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413
      }
    }

    // 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);

Jan Möbius's avatar
Jan Möbius committed
414 415 416
    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));
417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442

    // 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 )
{

  typename MeshT::VertexIter v_it;
  typename MeshT::VertexIter v_end = _mesh.vertices_end();

  for (v_it = _mesh.vertices_begin(); v_it != v_end; ++v_it){

Jan Möbius's avatar
Jan Möbius committed
443 444
    if (isBoundary(_mesh, *v_it)) continue;
    if ( isFeature(_mesh, *v_it)) continue;
445

Jan Möbius's avatar
Jan Möbius committed
446
    typename MeshT::Point p = _mesh.point(*v_it);
447 448 449 450 451
    typename MeshT::FaceHandle fhNear;
    double distance;

    typename MeshT::Point pNear = findNearestPoint(_original, p, fhNear, _ssearch, &distance);

Jan Möbius's avatar
Jan Möbius committed
452
    _mesh.set_point(*v_it, pNear);
453 454 455
  }
}