Developer Documentation

IsotropicRemesherT_impl.hh
1 /*===========================================================================*\
2 * *
3 * OpenFlipper *
4 * Copyright (c) 2001-2015, RWTH-Aachen University *
5 * Department of Computer Graphics and Multimedia *
6 * All rights reserved. *
7 * www.openflipper.org *
8 * *
9 *---------------------------------------------------------------------------*
10 * This file is part of OpenFlipper. *
11 *---------------------------------------------------------------------------*
12 * *
13 * Redistribution and use in source and binary forms, with or without *
14 * modification, are permitted provided that the following conditions *
15 * are met: *
16 * *
17 * 1. Redistributions of source code must retain the above copyright notice, *
18 * this list of conditions and the following disclaimer. *
19 * *
20 * 2. Redistributions in binary form must reproduce the above copyright *
21 * notice, this list of conditions and the following disclaimer in the *
22 * documentation and/or other materials provided with the distribution. *
23 * *
24 * 3. Neither the name of the copyright holder nor the names of its *
25 * contributors may be used to endorse or promote products derived from *
26 * this software without specific prior written permission. *
27 * *
28 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS *
29 * "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED *
30 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A *
31 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER *
32 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, *
33 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, *
34 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR *
35 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF *
36 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING *
37 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS *
38 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. *
39 * *
40 \*===========================================================================*/
41
42
43
44
45 #define ISOTROPICREMESHER_C
46
47 #include "IsotropicRemesherT.hh"
48
49 #include <ACG/Geometry/Algorithms.hh>
50
51 // -------------------- BSP
52 #include <ACG/Geometry/bsp/TriangleBSPT.hh>
53
55 template< class MeshT >
56 void IsotropicRemesher< MeshT >::remesh( MeshT& _mesh, const double _targetEdgeLength )
57 {
58 const double low = (4.0 / 5.0) * _targetEdgeLength;
59 const double high = (4.0 / 3.0) * _targetEdgeLength;
60
61 MeshT meshCopy = _mesh;
62 OpenMeshTriangleBSPT< MeshT >* triangleBSP = getTriangleBSP(meshCopy);
63
64 for (int i=0; i < 10; i++){
65 if (prgEmt_)
66 prgEmt_->sendProgressSignal(10*i + 0);
67 // std::cerr << "Iteration = " << i << std::endl;
68
69 splitLongEdges(_mesh, high);
70 if (prgEmt_)
71 prgEmt_->sendProgressSignal(10*i + 2);
72
73 // std::cerr << "collapse" << std::endl;
74 collapseShortEdges(_mesh, low, high);
75 if (prgEmt_)
76 prgEmt_->sendProgressSignal(10*i + 4);
77
78 // std::cerr << "equal" << std::endl;
79 equalizeValences(_mesh);
80 if (prgEmt_)
81 prgEmt_->sendProgressSignal(10*i + 6);
82
83 // std::cerr << "relax" << std::endl;
84 tangentialRelaxation(_mesh);
85 if (prgEmt_)
86 prgEmt_->sendProgressSignal(10*i + 8);
87
88 // std::cerr << "project" << std::endl;
89 projectToSurface(_mesh, meshCopy, triangleBSP);
90 if (prgEmt_)
91 prgEmt_->sendProgressSignal(10*i + 10);
92 }
93
94 delete triangleBSP;
95 }
96
97 template< class MeshT >
98 OpenMeshTriangleBSPT< MeshT >* IsotropicRemesher< MeshT >::getTriangleBSP(MeshT& _mesh)
99 {
100 // create Triangle BSP
101 OpenMeshTriangleBSPT< MeshT >* triangle_bsp = new OpenMeshTriangleBSPT< MeshT >( _mesh );
102
103 // build Triangle BSP
104 triangle_bsp->reserve(_mesh.n_faces());
105
106 typename MeshT::FIter f_it = _mesh.faces_begin();
107 typename MeshT::FIter f_end = _mesh.faces_end();
108
109 for (; f_it!=f_end; ++f_it)
110 triangle_bsp->push_back( *f_it );
111
112 triangle_bsp->build(10, 100);
113
114 // return pointer to triangle bsp
115 return triangle_bsp;
116 }
117
119 template< class MeshT >
120 void IsotropicRemesher< MeshT >::splitLongEdges( MeshT& _mesh, const double _maxEdgeLength )
121 {
122
123 const double _maxEdgeLengthSqr = _maxEdgeLength * _maxEdgeLength;
124
125 typename MeshT::EdgeIter e_it;
126 typename MeshT::EdgeIter e_end = _mesh.edges_end();
127
128 //iterate over all edges
129 for (e_it = _mesh.edges_begin(); e_it != e_end; ++e_it){
130 const typename MeshT::HalfedgeHandle & hh = _mesh.halfedge_handle( *e_it, 0 );
131
132 const typename MeshT::VertexHandle & v0 = _mesh.from_vertex_handle(hh);
133 const typename MeshT::VertexHandle & v1 = _mesh.to_vertex_handle(hh);
134
135 typename MeshT::Point vec = _mesh.point(v1) - _mesh.point(v0);
136
137 //edge to long?
138 if ( vec.sqrnorm() > _maxEdgeLengthSqr ){
139
140 const typename MeshT::Point midPoint = _mesh.point(v0) + ( 0.5 * vec );
141
142 //split at midpoint
143 typename MeshT::VertexHandle vh = _mesh.add_vertex( midPoint );
144
145 bool hadFeature = _mesh.status(*e_it).feature();
146
147 _mesh.split(*e_it, vh);
148
149 if ( hadFeature ){
150
151 typename MeshT::VOHIter vh_it;
152 for (vh_it = _mesh.voh_iter(vh); vh_it.is_valid(); ++vh_it)
153 if ( _mesh.to_vertex_handle(*vh_it) == v0 || _mesh.to_vertex_handle(*vh_it) == v1 )
154 _mesh.status( _mesh.edge_handle( *vh_it ) ).set_feature( true );
155 }
156 }
157 }
158 }
159
161 template< class MeshT >
162 void IsotropicRemesher< MeshT >::collapseShortEdges( MeshT& _mesh, const double _minEdgeLength, const double _maxEdgeLength )
163 {
164 const double _minEdgeLengthSqr = _minEdgeLength * _minEdgeLength;
165 const double _maxEdgeLengthSqr = _maxEdgeLength * _maxEdgeLength;
166
167 //add checked property
168 OpenMesh::EPropHandleT< bool > checked;
169 if ( !_mesh.get_property_handle(checked, "Checked Property") )
170 _mesh.add_property(checked,"Checked Property" );
171
172 //init property
173 typename MeshT::ConstEdgeIter e_it;
174 typename MeshT::ConstEdgeIter e_end = _mesh.edges_end();
175
176 for (e_it = _mesh.edges_begin(); e_it != e_end; ++e_it)
177 _mesh.property(checked, *e_it) = false;
178
179
180 bool finished = false;
181
182 while( !finished ){
183
184 finished = true;
185
186 for (e_it = _mesh.edges_begin(); e_it != _mesh.edges_end() ; ++e_it){
187
188 if ( _mesh.property(checked, *e_it) )
189 continue;
190
191 _mesh.property(checked, *e_it) = true;
192
193 const typename MeshT::HalfedgeHandle & hh = _mesh.halfedge_handle( *e_it, 0 );
194
195 const typename MeshT::VertexHandle & v0 = _mesh.from_vertex_handle(hh);
196 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 const double edgeLength = vec.sqrnorm();
201
202 // edge too short but don't try to collapse edges that have length 0
203 if ( (edgeLength < _minEdgeLengthSqr) && (edgeLength > DBL_EPSILON) ){
204
205 //check if the collapse is ok
206 const typename MeshT::Point & B = _mesh.point(v1);
207
208 bool collapse_ok = true;
209
210 for(typename MeshT::VOHIter vh_it = _mesh.voh_iter(v0); vh_it.is_valid(); ++vh_it)
211 if ( (( B - _mesh.point( _mesh.to_vertex_handle(*vh_it) ) ).sqrnorm() > _maxEdgeLengthSqr )
212 || ( _mesh.status( _mesh.edge_handle( *vh_it ) ).feature())
213 || ( _mesh.is_boundary( _mesh.edge_handle( *vh_it ) ) )){
214 collapse_ok = false;
215 break;
216 }
217
218 if( collapse_ok && _mesh.is_collapse_ok(hh) ) {
219 _mesh.collapse( hh );
220
221 finished = false;
222 }
223
224 }
225 }
226
227 }
228
229 _mesh.remove_property(checked);
230
231 _mesh.garbage_collection();
232 }
233
234 template< class MeshT >
235 void IsotropicRemesher< MeshT >::equalizeValences( MeshT& _mesh )
236 {
237
238 typename MeshT::EdgeIter e_it;
239 typename MeshT::EdgeIter e_end = _mesh.edges_end();
240
241 for (e_it = _mesh.edges_sbegin(); e_it != e_end; ++e_it){
242
243 if ( !_mesh.is_flip_ok(*e_it) ) continue;
244 if ( _mesh.status( *e_it ).feature() ) continue;
245
246
247 const typename MeshT::HalfedgeHandle & h0 = _mesh.halfedge_handle( *e_it, 0 );
248 const typename MeshT::HalfedgeHandle & h1 = _mesh.halfedge_handle( *e_it, 1 );
249
250 if (h0.is_valid() && h1.is_valid())
251
252 if (_mesh.face_handle(h0).is_valid() && _mesh.face_handle(h1).is_valid()){
253 //get vertices of corresponding faces
254 const typename MeshT::VertexHandle & a = _mesh.to_vertex_handle(h0);
255 const typename MeshT::VertexHandle & b = _mesh.to_vertex_handle(h1);
256 const typename MeshT::VertexHandle & c = _mesh.to_vertex_handle(_mesh.next_halfedge_handle(h0));
257 const typename MeshT::VertexHandle & d = _mesh.to_vertex_handle(_mesh.next_halfedge_handle(h1));
258
259 const int deviation_pre = abs((int)(_mesh.valence(a) - targetValence(_mesh, a)))
260 +abs((int)(_mesh.valence(b) - targetValence(_mesh, b)))
261 +abs((int)(_mesh.valence(c) - targetValence(_mesh, c)))
262 +abs((int)(_mesh.valence(d) - targetValence(_mesh, d)));
263 _mesh.flip(*e_it);
264
265 const int deviation_post = abs((int)(_mesh.valence(a) - targetValence(_mesh, a)))
266 +abs((int)(_mesh.valence(b) - targetValence(_mesh, b)))
267 +abs((int)(_mesh.valence(c) - targetValence(_mesh, c)))
268 +abs((int)(_mesh.valence(d) - targetValence(_mesh, d)));
269
270 if (deviation_pre <= deviation_post)
271 _mesh.flip(*e_it);
272 }
273 }
274 }
275
277 template< class MeshT >
278 inline
279 int IsotropicRemesher< MeshT >::targetValence( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
280
281 if (isBoundary(_mesh,_vh))
282 return 4;
283 else
284 return 6;
285 }
286
287 template< class MeshT >
288 inline
289 bool IsotropicRemesher< MeshT >::isBoundary( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
290
291 typename MeshT::ConstVertexOHalfedgeIter voh_it;
292
293 for (voh_it = _mesh.voh_iter( _vh ); voh_it.is_valid(); ++voh_it )
294 if ( _mesh.is_boundary( _mesh.edge_handle( *voh_it ) ) )
295 return true;
296
297 return false;
298 }
299
300 template< class MeshT >
301 inline
302 bool IsotropicRemesher< MeshT >::isFeature( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
303
304 typename MeshT::ConstVertexOHalfedgeIter voh_it;
305
306 for (voh_it = _mesh.voh_iter( _vh ); voh_it.is_valid(); ++voh_it )
307 if ( _mesh.status( _mesh.edge_handle( *voh_it ) ).feature() )
308 return true;
309
310 return false;
311 }
312
313 template< class MeshT >
314 void IsotropicRemesher< MeshT >::tangentialRelaxation( MeshT& _mesh )
315 {
316 _mesh.update_normals();
317
318 //add checked property
319 OpenMesh::VPropHandleT< typename MeshT::Point > q;
320 if ( !_mesh.get_property_handle(q, "q Property") )
321 _mesh.add_property(q,"q Property" );
322
323 typename MeshT::ConstVertexIter v_it;
324 typename MeshT::ConstVertexIter v_end = _mesh.vertices_end();
325
326 //first compute barycenters
327 for (v_it = _mesh.vertices_sbegin(); v_it != v_end; ++v_it){
328
329 typename MeshT::Point tmp(0.0, 0.0, 0.0);
330 uint N = 0;
331
332 typename MeshT::VVIter vv_it;
333 for (vv_it = _mesh.vv_iter(*v_it); vv_it.is_valid(); ++vv_it){
334 tmp += _mesh.point(*vv_it);
335 N++;
336 }
337
338 if (N > 0)
339 tmp /= (double) N;
340
341 _mesh.property(q, *v_it) = tmp;
342 }
343
344 //move to new position
345 for (v_it = _mesh.vertices_sbegin(); v_it != v_end; ++v_it){
346 if ( !isBoundary(_mesh, *v_it) && !isFeature(_mesh, *v_it) )
347 _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 _mesh.remove_property(q);
351 }
352
353 template <class MeshT>
354 template <class SpatialSearchT>
355 typename MeshT::Point
356 IsotropicRemesher< MeshT >::findNearestPoint(const MeshT& _mesh,
357 const typename MeshT::Point& _point,
358 typename MeshT::FaceHandle& _fh,
359 SpatialSearchT* _ssearch,
360 double* _dbest)
361 {
362
363 typename MeshT::Point p_best = _mesh.point(_mesh.vertex_handle(0));
364 typename MeshT::Scalar d_best = (_point-p_best).sqrnorm();
365
366 typename MeshT::FaceHandle fh_best;
367
368 if( _ssearch == 0 )
369 {
370 // exhaustive search
371 typename MeshT::ConstFaceIter cf_it = _mesh.faces_begin();
372 typename MeshT::ConstFaceIter cf_end = _mesh.faces_end();
373
374 for(; cf_it != cf_end; ++cf_it)
375 {
376 typename MeshT::ConstFaceVertexIter cfv_it = _mesh.cfv_iter(*cf_it);
377
378 const typename MeshT::Point& pt0 = _mesh.point( *cfv_it);
379 const typename MeshT::Point& pt1 = _mesh.point( *(++cfv_it));
380 const typename MeshT::Point& pt2 = _mesh.point( *(++cfv_it));
381
382 typename MeshT::Point ptn;
383
384 typename MeshT::Scalar d = ACG::Geometry::distPointTriangleSquared( _point,
385 pt0,
386 pt1,
387 pt2,
388 ptn );
389
390 if( d < d_best && d >= 0.0)
391 {
392 d_best = d;
393 p_best = ptn;
394
395 fh_best = *cf_it;
396 }
397 }
398
399 // return facehandle
400 _fh = fh_best;
401
402 // return distance
403 if(_dbest)
404 *_dbest = sqrt(d_best);
405
406
407 return p_best;
408 }
409 else
410 {
411 typename MeshT::FaceHandle fh = _ssearch->nearest(_point).handle;
412 typename MeshT::CFVIter fv_it = _mesh.cfv_iter(fh);
413
414 const typename MeshT::Point& pt0 = _mesh.point( *( fv_it));
415 const typename MeshT::Point& pt1 = _mesh.point( *(++fv_it));
416 const typename MeshT::Point& pt2 = _mesh.point( *(++fv_it));
417
418 // project
419 d_best = ACG::Geometry::distPointTriangleSquared(_point, pt0, pt1, pt2, p_best);
420
421 // return facehandle
422 _fh = fh;
423
424 // return distance
425 if(_dbest)
426 *_dbest = sqrt(d_best);
427
428 return p_best;
429 }
430 }
431
432
433 template< class MeshT >
434 template< class SpatialSearchT >
435 void IsotropicRemesher< MeshT >::projectToSurface( MeshT& _mesh, MeshT& _original, SpatialSearchT* _ssearch )
436 {
437
438 typename MeshT::VertexIter v_it;
439 typename MeshT::VertexIter v_end = _mesh.vertices_end();
440
441 for (v_it = _mesh.vertices_begin(); v_it != v_end; ++v_it){
442
443 if (isBoundary(_mesh, *v_it)) continue;
444 if ( isFeature(_mesh, *v_it)) continue;
445
446 typename MeshT::Point p = _mesh.point(*v_it);
447 typename MeshT::FaceHandle fhNear;
448 double distance;
449
450 typename MeshT::Point pNear = findNearestPoint(_original, p, fhNear, _ssearch, &distance);
451
452 _mesh.set_point(*v_it, pNear);
453 }
454 }
455
OpenMesh::VPropHandleT
Definition: Property.hh:484
IsotropicRemesher::remesh
void remesh(MeshT &_mesh, const double _targetEdgeLength)
do the remeshing
Definition: IsotropicRemesherT_impl.hh:56
OpenMesh::EPropHandleT
Definition: Property.hh:514
OpenMesh::VectorT::sqrnorm
Scalar sqrnorm(const VectorT< Scalar, DIM > &_v)
Definition: MeshNode2T_impl.hh:283
TriangleBSPCoreT< OpenMeshTriangleBSPTraits< Mesh > >::build
void build(unsigned int _max_handles, unsigned int _max_depth)
Definition: TriangleBSPCoreT_impl.hh:63
IsotropicRemesher
Definition: IsotropicRemesherT.hh:50
TriangleBSPCoreT< OpenMeshTriangleBSPTraits< Mesh > >::push_back
void push_back(Handle _h)
Add a handle to the BSP.
Definition: TriangleBSPCoreT.hh:97
IsotropicRemesher::targetValence
int targetValence(MeshT &_mesh, const typename MeshT::VertexHandle &_vh)
returns 4 for boundary vertices and 6 otherwise
Definition: IsotropicRemesherT_impl.hh:279
IsotropicRemesher::splitLongEdges
void splitLongEdges(MeshT &_mesh, const double _maxEdgeLength)
performs edge splits until all edges are shorter than the threshold
Definition: IsotropicRemesherT_impl.hh:120
IsotropicRemesher::collapseShortEdges
void collapseShortEdges(MeshT &_mesh, const double _minEdgeLength, const double _maxEdgeLength)
collapse edges shorter than minEdgeLength if collapsing doesn&#39;t result in new edge longer than maxEdg...
Definition: IsotropicRemesherT_impl.hh:162
TriangleBSPCoreT< OpenMeshTriangleBSPTraits< Mesh > >::reserve
void reserve(size_t _n)
Reserve memory for _n entries.
Definition: TriangleBSPCoreT.hh:95
OpenMeshTriangleBSPT
Definition: TriangleBSPT.hh:211