IsotropicRemesherT.cc 15.5 KB
Newer Older
Jan Möbius's avatar
Jan Möbius committed
1 2 3
/*===========================================================================*\
*                                                                            *
*                              OpenFlipper                                   *
Martin Schultz's avatar
Martin Schultz committed
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

Jan Möbius's avatar
 
Jan Möbius committed
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>
Jan Möbius's avatar
 
Jan Möbius committed
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);
Jan Möbius's avatar
 
Jan Möbius committed
67
// std::cerr << "Iteration = " << i << std::endl;
Anne Kathrein's avatar
Anne Kathrein committed
68 69

    splitLongEdges(_mesh, high);
70 71 72
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 2);
    
Jan Möbius's avatar
 
Jan Möbius committed
73 74
// std::cerr << "collapse" << std::endl;
    collapseShortEdges(_mesh, low, high);
75 76 77
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 4);
    
Jan Möbius's avatar
 
Jan Möbius committed
78 79
// std::cerr << "equal" << std::endl;
    equalizeValences(_mesh);
80 81 82
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 6);
    
Jan Möbius's avatar
 
Jan Möbius committed
83 84
// std::cerr << "relax" << std::endl;
    tangentialRelaxation(_mesh);
85 86 87
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 8);
    
Jan Möbius's avatar
 
Jan Möbius committed
88 89
// std::cerr << "project" << std::endl;
    projectToSurface(_mesh, meshCopy, triangleBSP);
90 91
    if (prgEmt_)
      prgEmt_->sendProgressSignal(10*i + 10);
Jan Möbius's avatar
 
Jan Möbius committed
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 );
Jan Möbius's avatar
 
Jan Möbius committed
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 );
Jan Möbius's avatar
 
Jan Möbius committed
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
      }
Jan Möbius's avatar
 
Jan Möbius committed
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;
Jan Möbius's avatar
 
Jan Möbius committed
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) )
Jan Möbius's avatar
 
Jan Möbius committed
189 190
          continue;
      
Jan Möbius's avatar
Jan Möbius committed
191
        _mesh.property(checked, *e_it) = true;
Jan Möbius's avatar
 
Jan Möbius committed
192

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

Jan Möbius's avatar
 
Jan Möbius committed
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) ){
Jan Möbius's avatar
 
Jan Möbius committed
204 205 206
    
          //check if the collapse is ok
          const typename MeshT::Point & B = _mesh.point(v1);
207

Jan Möbius's avatar
 
Jan Möbius committed
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 ) ) )){
Jan Möbius's avatar
 
Jan Möbius committed
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

Jan Möbius's avatar
 
Jan Möbius committed
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
    
Jan Möbius's avatar
 
Jan Möbius committed
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 );
Jan Möbius's avatar
 
Jan Möbius committed
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);
Jan Möbius's avatar
 
Jan Möbius committed
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);
Jan Möbius's avatar
 
Jan Möbius committed
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 ) ) )
Jan Möbius's avatar
 
Jan Möbius committed
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;
}

Jan Möbius's avatar
 
Jan Möbius committed
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);
Jan Möbius's avatar
 
Jan Möbius committed
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;
Jan Möbius's avatar
 
Jan Möbius committed
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));
Jan Möbius's avatar
 
Jan Möbius committed
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);
Jan Möbius's avatar
 
Jan Möbius committed
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));
Jan Möbius's avatar
 
Jan Möbius committed
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)
Jan Möbius's avatar
 
Jan Möbius committed
391 392 393 394
      {
        d_best = d;
        p_best = ptn;

Jan Möbius's avatar
Jan Möbius committed
395
        fh_best = *cf_it;
Jan Möbius's avatar
 
Jan Möbius committed
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));
Jan Möbius's avatar
 
Jan Möbius committed
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;
Jan Möbius's avatar
 
Jan Möbius committed
445

Jan Möbius's avatar
Jan Möbius committed
446
    typename MeshT::Point p = _mesh.point(*v_it);
Jan Möbius's avatar
 
Jan Möbius committed
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);
Jan Möbius's avatar
 
Jan Möbius committed
453 454 455
  }
}