Developer Documentation

- Main Page
- Related Pages
- Modules
- Namespaces
- Classes
- Files

- File List
- File Members

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Modules Pages
IsotropicRemesherT.cc
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 * $Revision$ *
45 * $LastChangedBy$ *
46 * $Date$ *
47 * *
48 \*===========================================================================*/
49
50
51 #define ISOTROPICREMESHER_C
52
53 #include "IsotropicRemesherT.hh"
54
55 #include <ACG/Geometry/Algorithms.hh>
56
57 // -------------------- BSP
58 #include <ACG/Geometry/bsp/TriangleBSPT.hh>
59
61 template< class MeshT >
62 void IsotropicRemesher< MeshT >::remesh( MeshT& _mesh, const double _targetEdgeLength )
63 {
64 const double low = (4.0 / 5.0) * _targetEdgeLength;
65 const double high = (4.0 / 3.0) * _targetEdgeLength;
66
67 MeshT meshCopy = _mesh;
68 OpenMeshTriangleBSPT< MeshT >* triangleBSP = getTriangleBSP(meshCopy);
69
70 for (int i=0; i < 10; i++){
71 if (prgEmt_)
72 prgEmt_->sendProgressSignal(10*i + 0);
73 // std::cerr << "Iteration = " << i << std::endl;
74
75 splitLongEdges(_mesh, high);
76 if (prgEmt_)
77 prgEmt_->sendProgressSignal(10*i + 2);
78
79 // std::cerr << "collapse" << std::endl;
80 collapseShortEdges(_mesh, low, high);
81 if (prgEmt_)
82 prgEmt_->sendProgressSignal(10*i + 4);
83
84 // std::cerr << "equal" << std::endl;
85 equalizeValences(_mesh);
86 if (prgEmt_)
87 prgEmt_->sendProgressSignal(10*i + 6);
88
89 // std::cerr << "relax" << std::endl;
90 tangentialRelaxation(_mesh);
91 if (prgEmt_)
92 prgEmt_->sendProgressSignal(10*i + 8);
93
94 // std::cerr << "project" << std::endl;
95 projectToSurface(_mesh, meshCopy, triangleBSP);
96 if (prgEmt_)
97 prgEmt_->sendProgressSignal(10*i + 10);
98 }
99
100 delete triangleBSP;
101 }
102
103 template< class MeshT >
104 OpenMeshTriangleBSPT< MeshT >* IsotropicRemesher< MeshT >::getTriangleBSP(MeshT& _mesh)
105 {
106 // create Triangle BSP
107 OpenMeshTriangleBSPT< MeshT >* triangle_bsp = new OpenMeshTriangleBSPT< MeshT >( _mesh );
108
109 // build Triangle BSP
110 triangle_bsp->reserve(_mesh.n_faces());
111
112 typename MeshT::FIter f_it = _mesh.faces_begin();
113 typename MeshT::FIter f_end = _mesh.faces_end();
114
115 for (; f_it!=f_end; ++f_it)
116 triangle_bsp->push_back( *f_it );
117
118 triangle_bsp->build(10, 100);
119
120 // return pointer to triangle bsp
121 return triangle_bsp;
122 }
123
125 template< class MeshT >
126 void IsotropicRemesher< MeshT >::splitLongEdges( MeshT& _mesh, const double _maxEdgeLength )
127 {
128
129 const double _maxEdgeLengthSqr = _maxEdgeLength * _maxEdgeLength;
130
131 typename MeshT::EdgeIter e_it;
132 typename MeshT::EdgeIter e_end = _mesh.edges_end();
133
134 //iterate over all edges
135 for (e_it = _mesh.edges_begin(); e_it != e_end; ++e_it){
136 const typename MeshT::HalfedgeHandle & hh = _mesh.halfedge_handle( *e_it, 0 );
137
138 const typename MeshT::VertexHandle & v0 = _mesh.from_vertex_handle(hh);
139 const typename MeshT::VertexHandle & v1 = _mesh.to_vertex_handle(hh);
140
141 typename MeshT::Point vec = _mesh.point(v1) - _mesh.point(v0);
142
143 //edge to long?
144 if ( vec.sqrnorm() > _maxEdgeLengthSqr ){
145
146 const typename MeshT::Point midPoint = _mesh.point(v0) + ( 0.5 * vec );
147
148 //split at midpoint
149 typename MeshT::VertexHandle vh = _mesh.add_vertex( midPoint );
150
151 bool hadFeature = _mesh.status(*e_it).feature();
152
153 _mesh.split(*e_it, vh);
154
155 if ( hadFeature ){
156
157 typename MeshT::VOHIter vh_it;
158 for (vh_it = _mesh.voh_iter(vh); vh_it.is_valid(); ++vh_it)
159 if ( _mesh.to_vertex_handle(*vh_it) == v0 || _mesh.to_vertex_handle(*vh_it) == v1 )
160 _mesh.status( _mesh.edge_handle( *vh_it ) ).set_feature( true );
161 }
162 }
163 }
164 }
165
167 template< class MeshT >
168 void IsotropicRemesher< MeshT >::collapseShortEdges( MeshT& _mesh, const double _minEdgeLength, const double _maxEdgeLength )
169 {
170 const double _minEdgeLengthSqr = _minEdgeLength * _minEdgeLength;
171 const double _maxEdgeLengthSqr = _maxEdgeLength * _maxEdgeLength;
172
173 //add checked property
174 OpenMesh::EPropHandleT< bool > checked;
175 if ( !_mesh.get_property_handle(checked, "Checked Property") )
176 _mesh.add_property(checked,"Checked Property" );
177
178 //init property
179 typename MeshT::ConstEdgeIter e_it;
180 typename MeshT::ConstEdgeIter e_end = _mesh.edges_end();
181
182 for (e_it = _mesh.edges_begin(); e_it != e_end; ++e_it)
183 _mesh.property(checked, *e_it) = false;
184
185
186 bool finished = false;
187
188 while( !finished ){
189
190 finished = true;
191
192 for (e_it = _mesh.edges_begin(); e_it != _mesh.edges_end() ; ++e_it){
193
194 if ( _mesh.property(checked, *e_it) )
195 continue;
196
197 _mesh.property(checked, *e_it) = true;
198
199 const typename MeshT::HalfedgeHandle & hh = _mesh.halfedge_handle( *e_it, 0 );
200
201 const typename MeshT::VertexHandle & v0 = _mesh.from_vertex_handle(hh);
202 const typename MeshT::VertexHandle & v1 = _mesh.to_vertex_handle(hh);
203
204 const typename MeshT::Point vec = _mesh.point(v1) - _mesh.point(v0);
205
206 const double edgeLength = vec.sqrnorm();
207
208 // edge too short but don't try to collapse edges that have length 0
209 if ( (edgeLength < _minEdgeLengthSqr) && (edgeLength > DBL_EPSILON) ){
210
211 //check if the collapse is ok
212 const typename MeshT::Point & B = _mesh.point(v1);
213
214 bool collapse_ok = true;
215
216 for(typename MeshT::VOHIter vh_it = _mesh.voh_iter(v0); vh_it.is_valid(); ++vh_it)
217 if ( (( B - _mesh.point( _mesh.to_vertex_handle(*vh_it) ) ).sqrnorm() > _maxEdgeLengthSqr )
218 || ( _mesh.status( _mesh.edge_handle( *vh_it ) ).feature())
219 || ( _mesh.is_boundary( _mesh.edge_handle( *vh_it ) ) )){
220 collapse_ok = false;
221 break;
222 }
223
224 if( collapse_ok && _mesh.is_collapse_ok(hh) ) {
225 _mesh.collapse( hh );
226
227 finished = false;
228 }
229
230 }
231 }
232
233 }
234
235 _mesh.remove_property(checked);
236
237 _mesh.garbage_collection();
238 }
239
240 template< class MeshT >
241 void IsotropicRemesher< MeshT >::equalizeValences( MeshT& _mesh )
242 {
243
244 typename MeshT::EdgeIter e_it;
245 typename MeshT::EdgeIter e_end = _mesh.edges_end();
246
247 for (e_it = _mesh.edges_sbegin(); e_it != e_end; ++e_it){
248
249 if ( !_mesh.is_flip_ok(*e_it) ) continue;
250 if ( _mesh.status( *e_it ).feature() ) continue;
251
252
253 const typename MeshT::HalfedgeHandle & h0 = _mesh.halfedge_handle( *e_it, 0 );
254 const typename MeshT::HalfedgeHandle & h1 = _mesh.halfedge_handle( *e_it, 1 );
255
256 if (h0.is_valid() && h1.is_valid())
257
258 if (_mesh.face_handle(h0).is_valid() && _mesh.face_handle(h1).is_valid()){
259 //get vertices of corresponding faces
260 const typename MeshT::VertexHandle & a = _mesh.to_vertex_handle(h0);
261 const typename MeshT::VertexHandle & b = _mesh.to_vertex_handle(h1);
262 const typename MeshT::VertexHandle & c = _mesh.to_vertex_handle(_mesh.next_halfedge_handle(h0));
263 const typename MeshT::VertexHandle & d = _mesh.to_vertex_handle(_mesh.next_halfedge_handle(h1));
264
265 const int deviation_pre = 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 _mesh.flip(*e_it);
270
271 const int deviation_post = abs((int)(_mesh.valence(a) - targetValence(_mesh, a)))
272 +abs((int)(_mesh.valence(b) - targetValence(_mesh, b)))
273 +abs((int)(_mesh.valence(c) - targetValence(_mesh, c)))
274 +abs((int)(_mesh.valence(d) - targetValence(_mesh, d)));
275
276 if (deviation_pre <= deviation_post)
277 _mesh.flip(*e_it);
278 }
279 }
280 }
281
283 template< class MeshT >
284 inline
285 int IsotropicRemesher< MeshT >::targetValence( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
286
287 if (isBoundary(_mesh,_vh))
288 return 4;
289 else
290 return 6;
291 }
292
293 template< class MeshT >
294 inline
295 bool IsotropicRemesher< MeshT >::isBoundary( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
296
297 typename MeshT::ConstVertexOHalfedgeIter voh_it;
298
299 for (voh_it = _mesh.voh_iter( _vh ); voh_it.is_valid(); ++voh_it )
300 if ( _mesh.is_boundary( _mesh.edge_handle( *voh_it ) ) )
301 return true;
302
303 return false;
304 }
305
306 template< class MeshT >
307 inline
308 bool IsotropicRemesher< MeshT >::isFeature( MeshT& _mesh, const typename MeshT::VertexHandle& _vh ){
309
310 typename MeshT::ConstVertexOHalfedgeIter voh_it;
311
312 for (voh_it = _mesh.voh_iter( _vh ); voh_it.is_valid(); ++voh_it )
313 if ( _mesh.status( _mesh.edge_handle( *voh_it ) ).feature() )
314 return true;
315
316 return false;
317 }
318
319 template< class MeshT >
320 void IsotropicRemesher< MeshT >::tangentialRelaxation( MeshT& _mesh )
321 {
322 _mesh.update_normals();
323
324 //add checked property
325 OpenMesh::VPropHandleT< typename MeshT::Point > q;
326 if ( !_mesh.get_property_handle(q, "q Property") )
327 _mesh.add_property(q,"q Property" );
328
329 typename MeshT::ConstVertexIter v_it;
330 typename MeshT::ConstVertexIter v_end = _mesh.vertices_end();
331
332 //first compute barycenters
333 for (v_it = _mesh.vertices_sbegin(); v_it != v_end; ++v_it){
334
335 typename MeshT::Point tmp(0.0, 0.0, 0.0);
336 uint N = 0;
337
338 typename MeshT::VVIter vv_it;
339 for (vv_it = _mesh.vv_iter(*v_it); vv_it.is_valid(); ++vv_it){
340 tmp += _mesh.point(*vv_it);
341 N++;
342 }
343
344 if (N > 0)
345 tmp /= (double) N;
346
347 _mesh.property(q, *v_it) = tmp;
348 }
349
350 //move to new position
351 for (v_it = _mesh.vertices_sbegin(); v_it != v_end; ++v_it){
352 if ( !isBoundary(_mesh, *v_it) && !isFeature(_mesh, *v_it) )
353 _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));
354 }
355
356 _mesh.remove_property(q);
357 }
358
359 template <class MeshT>
360 template <class SpatialSearchT>
361 typename MeshT::Point
362 IsotropicRemesher< MeshT >::findNearestPoint(const MeshT& _mesh,
363 const typename MeshT::Point& _point,
364 typename MeshT::FaceHandle& _fh,
365 SpatialSearchT* _ssearch,
366 double* _dbest)
367 {
368
369 typename MeshT::Point p_best = _mesh.point(_mesh.vertex_handle(0));
370 typename MeshT::Scalar d_best = (_point-p_best).sqrnorm();
371
372 typename MeshT::FaceHandle fh_best;
373
374 if( _ssearch == 0 )
375 {
376 // exhaustive search
377 typename MeshT::ConstFaceIter cf_it = _mesh.faces_begin();
378 typename MeshT::ConstFaceIter cf_end = _mesh.faces_end();
379
380 for(; cf_it != cf_end; ++cf_it)
381 {
382 typename MeshT::ConstFaceVertexIter cfv_it = _mesh.cfv_iter(*cf_it);
383
384 const typename MeshT::Point& pt0 = _mesh.point( *cfv_it);
385 const typename MeshT::Point& pt1 = _mesh.point( *(++cfv_it));
386 const typename MeshT::Point& pt2 = _mesh.point( *(++cfv_it));
387
388 typename MeshT::Point ptn;
389
390 typename MeshT::Scalar d = ACG::Geometry::distPointTriangleSquared( _point,
391 pt0,
392 pt1,
393 pt2,
394 ptn );
395
396 if( d < d_best && d >= 0.0)
397 {
398 d_best = d;
399 p_best = ptn;
400
401 fh_best = *cf_it;
402 }
403 }
404
405 // return facehandle
406 _fh = fh_best;
407
408 // return distance
409 if(_dbest)
410 *_dbest = sqrt(d_best);
411
412
413 return p_best;
414 }
415 else
416 {
417 typename MeshT::FaceHandle fh = _ssearch->nearest(_point).handle;
418 typename MeshT::CFVIter fv_it = _mesh.cfv_iter(fh);
419
420 const typename MeshT::Point& pt0 = _mesh.point( *( fv_it));
421 const typename MeshT::Point& pt1 = _mesh.point( *(++fv_it));
422 const typename MeshT::Point& pt2 = _mesh.point( *(++fv_it));
423
424 // project
425 d_best = ACG::Geometry::distPointTriangleSquared(_point, pt0, pt1, pt2, p_best);
426
427 // return facehandle
428 _fh = fh;
429
430 // return distance
431 if(_dbest)
432 *_dbest = sqrt(d_best);
433
434 return p_best;
435 }
436 }
437
438
439 template< class MeshT >
440 template< class SpatialSearchT >
441 void IsotropicRemesher< MeshT >::projectToSurface( MeshT& _mesh, MeshT& _original, SpatialSearchT* _ssearch )
442 {
443
444 typename MeshT::VertexIter v_it;
445 typename MeshT::VertexIter v_end = _mesh.vertices_end();
446
447 for (v_it = _mesh.vertices_begin(); v_it != v_end; ++v_it){
448
449 if (isBoundary(_mesh, *v_it)) continue;
450 if ( isFeature(_mesh, *v_it)) continue;
451
452 typename MeshT::Point p = _mesh.point(*v_it);
453 typename MeshT::FaceHandle fhNear;
454 double distance;
455
456 typename MeshT::Point pNear = findNearestPoint(_original, p, fhNear, _ssearch, &distance);
457
458 _mesh.set_point(*v_it, pNear);
459 }
460 }
461
IsotropicRemesher
Definition: IsotropicRemesherT.hh:56
IsotropicRemesher::collapseShortEdges
void collapseShortEdges(MeshT &_mesh, const double _minEdgeLength, const double _maxEdgeLength)
collapse edges shorter than minEdgeLength if collapsing doesn't result in new edge longer than maxEdg...
Definition: IsotropicRemesherT.cc:168
IsotropicRemesher::targetValence
int targetValence(MeshT &_mesh, const typename MeshT::VertexHandle &_vh)
returns 4 for boundary vertices and 6 otherwise
Definition: IsotropicRemesherT.cc:285
TriangleBSPCoreT< OpenMeshTriangleBSPTraits< Mesh > >::reserve
void reserve(size_t _n)
Reserve memory for _n entries.
Definition: TriangleBSPCoreT.hh:105
OpenMesh::EPropHandleT
Definition: Property.hh:515
ACG::Geometry::distPointTriangleSquared
Vec::value_type distPointTriangleSquared(const Vec &_p, const Vec &_v0, const Vec &_v1, const Vec &_v2, Vec &_nearestPoint)
squared distance from point _p to triangle (_v0, _v1, _v2)
Definition: Algorithms.cc:461
OpenMesh::VPropHandleT
Definition: Property.hh:487
IsotropicRemesher::remesh
void remesh(MeshT &_mesh, const double _targetEdgeLength)
do the remeshing
Definition: IsotropicRemesherT.cc:62
OpenMeshTriangleBSPT
Definition: TriangleBSPT.hh:206
TriangleBSPCoreT< OpenMeshTriangleBSPTraits< Mesh > >::push_back
void push_back(Handle _h)
Add a handle to the BSP.
Definition: TriangleBSPCoreT.hh:107
IsotropicRemesher::splitLongEdges
void splitLongEdges(MeshT &_mesh, const double _maxEdgeLength)
performs edge splits until all edges are shorter than the threshold
Definition: IsotropicRemesherT.cc:126
TriangleBSPCoreT< OpenMeshTriangleBSPTraits< Mesh > >::build
void build(unsigned int _max_handles, unsigned int _max_depth)