IsotropicRemesherT_impl.hh 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();

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