Developer Documentation
 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 >
105 {
106  // create Triangle BSP
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
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 >
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 >
321 {
322  _mesh.update_normals();
323 
324  //add checked property
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
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 
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...
int targetValence(MeshT &_mesh, const typename MeshT::VertexHandle &_vh)
returns 4 for boundary vertices and 6 otherwise
void reserve(size_t _n)
Reserve memory for _n entries.
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
void remesh(MeshT &_mesh, const double _targetEdgeLength)
do the remeshing
void push_back(Handle _h)
Add a handle to the BSP.
void splitLongEdges(MeshT &_mesh, const double _maxEdgeLength)
performs edge splits until all edges are shorter than the threshold
void build(unsigned int _max_handles, unsigned int _max_depth)