DGtal  1.5.beta
SurfaceMesh.ih
1 /**
2  * This program is free software: you can redistribute it and/or modify
3  * it under the terms of the GNU Lesser General Public License as
4  * published by the Free Software Foundation, either version 3 of the
5  * License, or (at your option) any later version.
6  *
7  * This program is distributed in the hope that it will be useful,
8  * but WITHOUT ANY WARRANTY; without even the implied warranty of
9  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
10  * GNU General Public License for more details.
11  *
12  * You should have received a copy of the GNU General Public License
13  * along with this program. If not, see <http://www.gnu.org/licenses/>.
14  *
15  **/
16 
17 /**
18  * @file SurfaceMesh.ih
19  * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20  * Laboratory of Mathematics (CNRS, UMR 5127), University of Savoie, France
21  *
22  * @date 2020/02/18
23  *
24  * Implementation of inline methods defined in SurfaceMesh.h
25  *
26  * This file is part of the DGtal library.
27  */
28 
29 
30 //////////////////////////////////////////////////////////////////////////////
31 #include <cstdlib>
32 #include <limits>
33 //////////////////////////////////////////////////////////////////////////////
34 
35 ///////////////////////////////////////////////////////////////////////////////
36 // IMPLEMENTATION of inline methods.
37 ///////////////////////////////////////////////////////////////////////////////
38 
39 ///////////////////////////////////////////////////////////////////////////////
40 // ----------------------- Standard services ------------------------------
41 
42 //-----------------------------------------------------------------------------
43 template <typename TRealPoint, typename TRealVector>
44 template <typename RealPointIterator, typename VerticesIterator>
45 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
46 SurfaceMesh( RealPointIterator itPos, RealPointIterator itPosEnd,
47  VerticesIterator itVertices, VerticesIterator itVerticesEnd )
48 {
49  bool ok = init( itPos, itPosEnd, itVertices, itVerticesEnd );
50  if ( !ok ) clear();
51 }
52 
53 //-----------------------------------------------------------------------------
54 template <typename TRealPoint, typename TRealVector>
55 template <typename RealPointIterator, typename VerticesIterator>
56 bool
57 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
58 init( RealPointIterator itPos, RealPointIterator itPosEnd,
59  VerticesIterator itVertices, VerticesIterator itVerticesEnd )
60 {
61  clear();
62  myPositions = std::vector< RealPoint >( itPos, itPosEnd );
63  myIncidentFaces.resize( myPositions.size() );
64  Index f = 0; // current face index
65  bool ok = true;
66  for ( ; itVertices != itVerticesEnd; ++itVertices, ++f )
67  {
68  Vertices f_vtcs;
69  for ( auto it = itVertices->begin(), itE = itVertices->end(); it != itE; ++it )
70  {
71  Index vtx = *it;
72  if ( vtx >= nbVertices() )
73  {
74  trace.warning() << "[SurfaceMesh::init] Invalid vtx "
75  << vtx << " at face " << f
76  << " since #V=" << nbVertices()
77  << ". Ignoring vertex." << std::endl;
78  ok = false;
79  }
80  else
81  {
82  myIncidentFaces[ vtx ].push_back( f );
83  f_vtcs.push_back( vtx );
84  }
85  }
86  myIncidentVertices.push_back( f_vtcs );
87  }
88  computeNeighbors();
89  computeEdges();
90  return ok;
91 }
92 
93 //-----------------------------------------------------------------------------
94 template <typename TRealPoint, typename TRealVector>
95 void
96 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
97 clear()
98 {
99  myIncidentVertices.clear();
100  myIncidentFaces.clear();
101  myPositions.clear();
102  myVertexNormals.clear();
103  myFaceNormals.clear();
104  myNeighborFaces.clear();
105  myNeighborVertices.clear();
106  myEdgeVertices.clear();
107  myVertexPairEdge.clear();
108  myEdgeFaces.clear();
109  myEdgeRightFaces.clear();
110  myEdgeLeftFaces.clear();
111 }
112 
113 //-----------------------------------------------------------------------------
114 template <typename TRealPoint, typename TRealVector>
115 template <typename RealVectorIterator>
116 bool
117 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
118 setVertexNormals( RealVectorIterator itN, RealVectorIterator itNEnd )
119 {
120  myVertexNormals = std::vector< RealVector >( itN, itNEnd );
121  return myVertexNormals.size() == myPositions.size();
122 }
123 
124 //-----------------------------------------------------------------------------
125 template <typename TRealPoint, typename TRealVector>
126 template <typename RealVectorIterator>
127 bool
128 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
129 setFaceNormals( RealVectorIterator itN, RealVectorIterator itNEnd )
130 {
131  myFaceNormals = std::vector< RealVector >( itN, itNEnd );
132  return myFaceNormals.size() == myIncidentVertices.size();
133 }
134 
135 //-----------------------------------------------------------------------------
136 template <typename TRealPoint, typename TRealVector>
137 void
138 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
139 computeFaceNormalsFromPositions()
140 {
141  myFaceNormals.resize( myIncidentVertices.size() );
142  for ( Face f = 0; f < myFaceNormals.size(); f++ )
143  computeFaceNormalFromPositions( f );
144 }
145 
146 //-----------------------------------------------------------------------------
147 template <typename TRealPoint, typename TRealVector>
148 void
149 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
150 computeFaceNormalFromPositions( const Face f )
151 {
152  RealPoint p; // barycenter
153  RealVector n; // normal
154  const Vertices& vtcs = incidentVertices( f );
155  // compute barycenter
156  for ( auto idx : vtcs ) p += myPositions[ idx ];
157  p /= vtcs.size();
158  // compute normal as sum of triangle normal vectors.
159  for ( Index i = 0; i < vtcs.size(); ++i )
160  {
161  const Index j = vtcs[ i ];
162  const Index nj = vtcs[ (i+1) % vtcs.size() ];
163  n += ( myPositions[ j ] - p ).crossProduct( myPositions[ nj ] - p );
164  }
165  auto n_norm = n.norm();
166  myFaceNormals[ f ] = n_norm != 0.0 ? n / n_norm : n;
167 }
168 
169 //-----------------------------------------------------------------------------
170 template <typename TRealPoint, typename TRealVector>
171 void
172 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
173 computeFaceNormalsFromVertexNormals()
174 {
175  if ( myVertexNormals.empty() ) return;
176  myFaceNormals.resize( myIncidentVertices.size() );
177  Index f = 0;
178  for ( auto face : myIncidentVertices )
179  {
180  RealVector n; // normal
181  for ( auto idx : face ) n += myVertexNormals[ idx ];
182  auto n_norm = n.norm();
183  myFaceNormals[ f ] = n_norm != 0.0 ? n / n_norm : n;
184  f++;
185  }
186 }
187 //-----------------------------------------------------------------------------
188 template <typename TRealPoint, typename TRealVector>
189 void
190 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
191 computeVertexNormalsFromFaceNormals()
192 {
193  if ( myFaceNormals.empty() ) return;
194  myVertexNormals.resize( myIncidentFaces.size() );
195  Index v = 0;
196  for ( auto vertex : myIncidentFaces )
197  {
198  RealVector n; // normal
199  for ( auto idx : vertex ) n += myFaceNormals[ idx ];
200  auto n_norm = n.norm();
201  myVertexNormals[ v ] = n_norm != 0.0 ? n / n_norm : n;
202  v++;
203  }
204 }
205 //-----------------------------------------------------------------------------
206 template <typename TRealPoint, typename TRealVector>
207 void
208 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
209 computeVertexNormalsFromFaceNormalsWithMaxWeights()
210 {
211  if ( myFaceNormals.empty() ) return;
212  myVertexNormals.resize( myIncidentFaces.size() );
213  Index v = 0;
214  for ( auto incident_faces : myIncidentFaces )
215  {
216  RealVector n; // normal
217  const auto weights = getMaxWeights( v );
218  Index i = 0;
219  for ( auto idx_f : incident_faces ) n += weights[ i++ ] * myFaceNormals[ idx_f ];
220  auto n_norm = n.norm();
221  myVertexNormals[ v ] = n_norm != 0.0 ? n / n_norm : n;
222  v++;
223  }
224 }
225 
226 //-----------------------------------------------------------------------------
227 template <typename TRealPoint, typename TRealVector>
228 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalars
229 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
230 getMaxWeights( Index v ) const
231 {
232  Scalars weights;
233  const auto & neighbors = myNeighborVertices[ v ];
234  const RealPoint x = myPositions[ v ];
235  for ( auto idx_f : myIncidentFaces[ v ] )
236  {
237  // Find adjacent vertices to v
238  std::vector< Index > adj_vertices;
239  for ( auto idx_v : myIncidentVertices[ idx_f ] )
240  {
241  auto it = std::find( neighbors.cbegin(), neighbors.cend(), idx_v );
242  if ( it != neighbors.cend() ) adj_vertices.push_back( *it );
243  }
244  if ( adj_vertices.size() != 2 )
245  {
246  trace.warning() << "[SurfaceMesh::getMaxWeights] "
247  << adj_vertices.size() << " adjacent vertices to vertex "
248  << v << " on face" << idx_f << "." << std::endl;
249  for ( auto a : adj_vertices ) std::cerr << " " << a;
250  std::cerr << std::endl;
251  }
252  if (adj_vertices.size() >= 2 )
253  {
254  const Scalar area = faceArea( idx_f );
255  const Scalar l1 = ( myPositions[ adj_vertices[ 0 ] ] - x ).squaredNorm();
256  const Scalar l2 = ( myPositions[ adj_vertices[ 1 ] ] - x ).squaredNorm();
257  const Scalar l1l2 = l1 * l2;
258  weights.push_back( l1l2 != 0 ? fabs( area ) / l1l2 : 0.0 );
259  }
260  else weights.push_back( 0.0 );
261  }
262  return weights;
263 }
264 
265 //-----------------------------------------------------------------------------
266 template <typename TRealPoint, typename TRealVector>
267 template <typename AnyRing>
268 std::vector<AnyRing>
269 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
270 computeFaceValuesFromVertexValues( const std::vector<AnyRing>& vvalues ) const
271 {
272  ASSERT( vvalues.size() == nbVertices() );
273  std::vector<AnyRing> fvalues( nbFaces() );
274  Index f = 0;
275  for ( auto face : myIncidentVertices )
276  {
277  AnyRing n = NumberTraits<AnyRing>::ZERO;
278  for ( auto idx : face ) n += vvalues[ idx ];
279  fvalues[ f++ ] = n / face.size();
280  }
281  return fvalues;
282 }
283 
284 template <typename TRealPoint, typename TRealVector>
285 template <typename AnyRing>
286 std::vector<AnyRing>
287 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
288 computeVertexValuesFromFaceValues( const std::vector<AnyRing>& fvalues ) const
289 {
290  ASSERT( fvalues.size() == nbFaces() );
291  std::vector<AnyRing> vvalues( nbVertices() );
292  Index v = 0;
293  for ( auto vertex : myIncidentFaces )
294  {
295  AnyRing n = NumberTraits<AnyRing>::ZERO;
296  for ( auto idx : vertex ) n += fvalues[ idx ];
297  vvalues[ v++ ] = n / vertex.size();
298  }
299  return vvalues;
300 }
301 
302 //-----------------------------------------------------------------------------
303 template <typename TRealPoint, typename TRealVector>
304 std::vector<TRealVector>
305 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
306 computeFaceUnitVectorsFromVertexUnitVectors
307 ( const std::vector<RealVector>& vuvectors ) const
308 {
309  ASSERT( vuvectors.size() == nbVertices() );
310  std::vector<RealVector> fuvectors( nbFaces() );
311  Index f = 0;
312  for ( auto face : myIncidentVertices )
313  {
314  RealVector n;
315  for ( auto idx : face ) n += vuvectors[ idx ];
316  const auto n_norm = n.norm();
317  fuvectors[ f++ ] = n_norm != 0.0 ? n / n_norm : n;
318  }
319  return fuvectors;
320 }
321 
322 //-----------------------------------------------------------------------------
323 template <typename TRealPoint, typename TRealVector>
324 std::vector<TRealVector>
325 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
326 computeVertexUnitVectorsFromFaceUnitVectors
327 ( const std::vector<RealVector>& fuvectors ) const
328 {
329  ASSERT( fuvectors.size() == nbFaces() );
330  std::vector<RealVector> vuvectors( nbVertices() );
331  Index v = 0;
332  for ( auto vertex : myIncidentFaces )
333  {
334  RealVector n;
335  for ( auto idx : vertex ) n += fuvectors[ idx ];
336  const auto n_norm = n.norm();
337  vuvectors[ v++ ] = n_norm != 0.0 ? n / n_norm : n;
338  }
339  return vuvectors;
340 }
341 
342 
343 //-----------------------------------------------------------------------------
344 template <typename TRealPoint, typename TRealVector>
345 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edge
346 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
347 makeEdge( Vertex i, Vertex j ) const
348 {
349  VertexPair vp = i < j ? std::make_pair( i,j ) : std::make_pair( j,i );
350  const auto it = myVertexPairEdge.find( vp );
351  if ( it == myVertexPairEdge.cend() ) return nbEdges();
352  return it->second;
353 }
354 
355 //-----------------------------------------------------------------------------
356 template <typename TRealPoint, typename TRealVector>
357 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
358 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
359 averageEdgeLength() const
360 {
361  double lengths = 0.0;
362  for ( Edge e = 0; e < nbEdges(); ++e )
363  {
364  auto vtcs = edgeVertices( e );
365  const RealPoint p = myPositions[ vtcs.first ];
366  const RealVector pq = myPositions[ vtcs.second ] - p;
367  lengths += pq.norm();
368  }
369  lengths /= nbEdges();
370  return lengths;
371 }
372 
373 //-----------------------------------------------------------------------------
374 template <typename TRealPoint, typename TRealVector>
375 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
376 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
377 localWindow( Face f ) const
378 {
379  const RealPoint x = faceCentroid( f );
380  Scalar local_length = 0.0;
381  for ( auto v : incidentVertices( f ) )
382  local_length += ( myPositions[ v ] - x ).norm();
383  return local_length / incidentVertices( f ).size();
384 }
385 
386 //-----------------------------------------------------------------------------
387 template <typename TRealPoint, typename TRealVector>
388 void
389 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
390 perturbateWithUniformRandomNoise( Scalar p )
391 {
392  for ( auto& x : myPositions )
393  {
394  RealVector d( rand01()*2.0 - 1.0, rand01()*2.0 - 1.0, rand01()*2.0 - 1.0 );
395  d = d.getNormalized();
396  Scalar l = rand01() * p;
397  x += l * d;
398  }
399 }
400 
401 //-----------------------------------------------------------------------------
402 template <typename TRealPoint, typename TRealVector>
403 void
404 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
405 perturbateWithAdaptiveUniformRandomNoise( Scalar p )
406 {
407  Scalars local_average_lengths( nbVertices() );
408  for ( Index v = 0; v < nbVertices(); ++v )
409  {
410  const RealPoint x = myPositions[ v ];
411  Scalar local_length = 0.0;
412  for ( auto nv : myNeighborVertices[ v ] )
413  local_length += ( myPositions[ nv ] - x ).norm();
414  local_average_lengths[ v ] = local_length / myNeighborVertices[ v ].size();
415  }
416  for ( Index v = 0; v < nbVertices(); ++v )
417  {
418  RealVector d( rand01()*2.0 - 1.0, rand01()*2.0 - 1.0, rand01()*2.0 - 1.0 );
419  d = d.getNormalized();
420  Scalar l = rand01() * p * local_average_lengths[ v ];
421  myPositions[ v ] += l * d;
422  }
423 }
424 
425 //-----------------------------------------------------------------------------
426 template <typename TRealPoint, typename TRealVector>
427 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::RealPoint
428 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
429 faceCentroid( Index f ) const
430 {
431  RealPoint c;
432  for ( auto v : myIncidentVertices[ f ] )
433  c += myPositions[ v ];
434  return c / myIncidentVertices[ f ].size();
435 }
436 
437 //-----------------------------------------------------------------------------
438 template <typename TRealPoint, typename TRealVector>
439 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::RealPoint
440 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
441 edgeCentroid( Index e ) const
442 {
443  const auto& vtcs = edgeVertices( e );
444  RealPoint c = myPositions[ vtcs.first ] + myPositions[ vtcs.second ];
445  return c / 2.0;
446 }
447 
448 //-----------------------------------------------------------------------------
449 template <typename TRealPoint, typename TRealVector>
450 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
451 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
452 faceArea( Index f ) const
453 {
454  Scalar area = 0.0;
455  const auto & inc_vtcs = myIncidentVertices[ f ];
456  RealPoint p = myPositions[ inc_vtcs.back() ];
457  const Index m = inc_vtcs.size() - 2;
458  for ( Index i = 0; i < m; ++i )
459  area += ( myPositions[ inc_vtcs[ i ] ] - p )
460  .crossProduct( myPositions[ inc_vtcs[ i+1 ] ] - p ).norm();
461  return area / 2.0;
462 }
463 
464 //-----------------------------------------------------------------------------
465 template <typename TRealPoint, typename TRealVector>
466 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces
467 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
468 computeFacesInclusionsInBall( Scalar r, Index f ) const
469 {
470  return computeFacesInclusionsInBall( r, f, faceCentroid( f ) );
471 }
472 
473 //-----------------------------------------------------------------------------
474 template <typename TRealPoint, typename TRealVector>
475 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces
476 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
477 computeFacesInclusionsInBall( Scalar r, Index f, RealPoint p ) const
478 {
479  WeightedFaces result;
480  if ( r < 0.000001 )
481  {
482  result.push_back( std::make_pair( f, 0.000001 ) );
483  return result;
484  }
485  std::unordered_set< Index > marked;
486  std::queue< Index > active;
487  active.push( f );
488  marked.insert( f );
489  while ( ! active.empty() )
490  {
491  Index current = active.front();
492  active.pop();
493  Scalar weight = faceInclusionRatio( p, r, current );
494  if ( weight > 0.0 )
495  {
496  result.push_back( std::make_pair( current, weight ) );
497  auto neighbors = myNeighborFaces[ current ];
498  for ( auto n : neighbors )
499  if ( marked.find( n ) == marked.end() )
500  {
501  active.push( n );
502  marked.insert( n );
503  }
504  }
505  }
506  return result;
507 }
508 
509 //-----------------------------------------------------------------------------
510 template <typename TRealPoint, typename TRealVector>
511 std::tuple
512 < typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Vertices,
513  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedEdges,
514  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces >
515 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
516 computeCellsInclusionsInBall( Scalar r, Index f ) const
517 {
518  return computeCellsInclusionsInBall( r, f, faceCentroid( f ) );
519 }
520 
521 //-----------------------------------------------------------------------------
522 template <typename TRealPoint, typename TRealVector>
523 std::tuple
524 < typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Vertices,
525  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedEdges,
526  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces >
527 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
528 computeCellsInclusionsInBall( Scalar r, Index f, RealPoint p ) const
529 {
530  Vertices result_v;
531  std::set<Vertex> set_v;
532  WeightedEdges result_e;
533  WeightedFaces result_f;
534  if ( r < 0.000001 )
535  {
536  result_f.push_back( std::make_pair( f, 0.000001 ) );
537  return std::make_tuple( result_v, result_e, result_f );
538  }
539  std::set< Index > marked;
540  std::queue< Index > active;
541  active.push( f );
542  marked.insert( f );
543  while ( ! active.empty() )
544  {
545  Index current = active.front();
546  active.pop();
547  Scalar fweight = faceInclusionRatio( p, r, current );
548  if ( fweight > 0.0 )
549  {
550  result_f.push_back( std::make_pair( current, fweight ) );
551  // Taking care of faces, and the breadth-first traversal
552  const auto& neighbors = myNeighborFaces[ current ];
553  for ( auto n : neighbors )
554  if ( marked.find( n ) == marked.end() )
555  {
556  active.push( n );
557  marked.insert( n );
558  }
559  // Taking care of edges and vertices
560  const auto& inc_v = myIncidentVertices[ current ];
561  for ( Size i = 0; i < inc_v.size(); ++i )
562  {
563  const Vertex vi = inc_v[ i ];
564  const Vertex vn = inc_v[ (i+1) % inc_v.size() ];
565  if ( vertexInclusionRatio( p, r, vi ) > 0.0 )
566  set_v.insert( vi );
567  if ( vn < vi ) continue; // edges are ordered pairs
568  const Edge e_ij = makeEdge( vi, vn );
569  if ( e_ij >= nbEdges() ) {
570  trace.error() << "bad edge " << vi << " " << vn << std::endl;
571  continue;
572  }
573  Scalar eweight = edgeInclusionRatio( p, r, e_ij );
574  if ( eweight > 0.0 )
575  result_e.push_back( std::make_pair( e_ij, eweight ) );
576  }
577  }
578  }
579  result_v = Vertices( set_v.cbegin(), set_v.cend() );
580  return std::make_tuple( result_v, result_e, result_f );
581 }
582 
583 //-----------------------------------------------------------------------------
584 template <typename TRealPoint, typename TRealVector>
585 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
586 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
587 faceInclusionRatio( RealPoint p, Scalar r, Index f ) const
588 {
589  const auto vertices = myIncidentVertices[ f ];
590  const RealPoint b = faceCentroid( f );
591  Scalar d_min = ( b - p ).norm();
592  Scalar d_max = d_min;
593  for ( auto v : vertices )
594  {
595  Scalar d = ( myPositions[ v ] - p ).norm();
596  d_max = std::max( d_max, d );
597  d_min = std::min( d_min, d );
598  }
599  if ( d_max <= r ) return 1.0;
600  else if ( r <= d_min ) return 0.0;
601  return ( r - d_min ) / ( d_max - d_min );
602 }
603 
604 //-----------------------------------------------------------------------------
605 template <typename TRealPoint, typename TRealVector>
606 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
607 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
608 edgeInclusionRatio( RealPoint p, Scalar r, Index e ) const
609 {
610  const auto vertices = edgeVertices( e );
611  const RealPoint b = edgeCentroid( e );
612  const Scalar d0 = ( myPositions[ vertices.first ] - p ).norm();
613  const Scalar d1 = ( myPositions[ vertices.second ] - p ).norm();
614  Scalar d_min = ( b - p ).norm();
615  Scalar d_max = d_min;
616  d_max = std::max( d_max, std::max( d0, d1 ) );
617  d_min = std::min( d_min, std::min( d0, d1 ) );
618  if ( d_max <= r ) return 1.0;
619  else if ( r <= d_min ) return 0.0;
620  return ( r - d_min ) / ( d_max - d_min );
621 }
622 
623 //-----------------------------------------------------------------------------
624 template <typename TRealPoint, typename TRealVector>
625 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
626 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
627 vertexInclusionRatio( RealPoint p, Scalar r, Index v ) const
628 {
629  const RealPoint b = myPositions[ v ];
630  return ( ( b - p ).norm() <= r ) ? 1.0 : 0.0;
631 }
632 
633 
634 ///////////////////////////////////////////////////////////////////////////////
635 // Interface - public :
636 
637 //-----------------------------------------------------------------------------
638 template <typename TRealPoint, typename TRealVector>
639 void
640 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
641 selfDisplay ( std::ostream & out ) const
642 {
643  out << "[SurfaceMesh" << ( isValid() ? " (OK)" : " (KO)" )
644  << " #V=" << myPositions.size()
645  << " #VN=" << myVertexNormals.size()
646  << " #E=" << myEdgeVertices.size()
647  << " #F=" << myIncidentVertices.size()
648  << " #FN=" << myFaceNormals.size();
649  double nb_nf = 0.0;
650  double nb_nv = 0.0;
651  double nb_nfe = 0.0;
652  for ( auto nf : myNeighborFaces ) nb_nf += nf.size();
653  for ( auto nv : myNeighborVertices ) nb_nv += nv.size();
654  for ( auto nfe : myEdgeFaces ) nb_nfe += nfe.size();
655  nb_nf /= nbFaces();
656  nb_nv /= nbVertices();
657  nb_nfe /= nbEdges();
658  out << " E[IF]=" << nb_nf << " E[IV]=" << nb_nv << " E[IFE]=" << nb_nfe;
659  out << "]";
660 }
661 
662 //-----------------------------------------------------------------------------
663 template <typename TRealPoint, typename TRealVector>
664 bool
665 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
666 isValid() const
667 {
668  return myPositions.size() == myIncidentFaces.size()
669  && ( myVertexNormals.size() == 0
670  || ( myVertexNormals.size() == myPositions.size() ) )
671  && ( myFaceNormals.size() == 0
672  || ( myFaceNormals.size() == myIncidentVertices.size() ) );
673 }
674 
675 
676 //-----------------------------------------------------------------------------
677 template <typename TRealPoint, typename TRealVector>
678 void
679 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
680 computeNeighbors()
681 {
682  myNeighborFaces .resize( nbFaces() );
683  myNeighborVertices.resize( nbVertices() );
684  std::vector< std::set< Index > > tmp( nbVertices() );
685  // For each vertex, computes its neighboring vertices
686  for ( auto incident_vertices : allIncidentVertices() )
687  {
688  const Size nb_iv = incident_vertices.size();
689  for ( Size k = 0; k < nb_iv; ++k )
690  {
691  tmp[ incident_vertices[ k ] ].insert( incident_vertices[ (k+1)%nb_iv ] );
692  tmp[ incident_vertices[ (k+1)%nb_iv ] ].insert( incident_vertices[ k ] );
693  }
694  }
695  // For each vertex, computes its neighboring vertices
696  for ( Index idx_v = 0; idx_v < nbVertices(); ++idx_v )
697  myNeighborVertices[ idx_v ] = Vertices( tmp[ idx_v ].cbegin(), tmp[ idx_v ].cend() );
698 
699  // For each face, computes its neighboring faces
700  Index idx_f = 0;
701  for ( auto incident_vertices : allIncidentVertices() )
702  {
703  std::set< Index > neighbor_faces_set;
704  std::sort( incident_vertices.begin(), incident_vertices.end() );
705  for ( auto idx_v : incident_vertices )
706  {
707  const auto & incident_faces = incidentFaces( idx_v );
708  for ( auto inc_f : incident_faces )
709  {
710  // Keep only faces incident to two vertices of f.
711  auto incident_vertices2 = incidentVertices( inc_f );
712  std::sort( incident_vertices2.begin(), incident_vertices2.end() );
713  Vertices common;
714  std::set_intersection( incident_vertices.cbegin(), incident_vertices.cend(),
715  incident_vertices2.cbegin(), incident_vertices2.cend(),
716  std::back_inserter( common ) );
717  if ( common.size() == 2 )
718  neighbor_faces_set.insert( inc_f );
719  }
720  }
721  neighbor_faces_set.erase( idx_f );
722  Faces neighbor_faces( neighbor_faces_set.begin(), neighbor_faces_set.end() );
723  myNeighborFaces[ idx_f++ ] = neighbor_faces;
724  }
725 }
726 
727 //-----------------------------------------------------------------------------
728 template <typename TRealPoint, typename TRealVector>
729 void
730 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
731 computeEdges()
732 {
733  std::map< VertexPair, std::vector<Face> > edge2face_right;
734  std::map< VertexPair, std::vector<Face> > edge2face_left;
735  std::set< VertexPair > edges;
736  Index idx_f = 0;
737  for ( auto incident_vertices : allIncidentVertices() )
738  {
739  const Size n = incident_vertices.size();
740  for ( Size i = 0; i < n; i++ )
741  {
742  VertexPair e = std::make_pair( incident_vertices[ i ],
743  incident_vertices[ (i+1) % n ] );
744  if ( e.first < e.second ) {
745  edge2face_left[ e ].push_back( idx_f );
746  } else {
747  std::swap( e.first, e.second );
748  edge2face_right[ e ].push_back( idx_f );
749  }
750  edges.insert( e );
751  }
752  idx_f++;
753  }
754  const Size nbe = edges.size();
755  myVertexPairEdge.clear();
756  myEdgeVertices.resize ( nbe );
757  myEdgeFaces.resize ( nbe );
758  myEdgeRightFaces.resize( nbe );
759  myEdgeLeftFaces.resize ( nbe );
760  Index idx_e = 0;
761  for ( auto e : edges )
762  {
763  myEdgeVertices [ idx_e ] = e;
764  myVertexPairEdge[ e ] = idx_e;
765  myEdgeLeftFaces [ idx_e ] = edge2face_left [ e ];
766  myEdgeRightFaces[ idx_e ] = edge2face_right[ e ];
767  myEdgeFaces [ idx_e ] = myEdgeRightFaces[ idx_e ];
768  myEdgeFaces [ idx_e ].insert( myEdgeFaces[ idx_e ].end(),
769  myEdgeLeftFaces[ idx_e ].cbegin(),
770  myEdgeLeftFaces[ idx_e ].cend() );
771  idx_e++;
772  }
773 }
774 
775 //-----------------------------------------------------------------------------
776 template <typename TRealPoint, typename TRealVector>
777 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
778 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
779 computeManifoldBoundaryEdges() const
780 {
781  Edges edges;
782  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
783  {
784  if ( myEdgeFaces[ e ].size() == 1 )
785  edges.push_back( e );
786  }
787  return edges;
788 }
789 
790 //-----------------------------------------------------------------------------
791 template <typename TRealPoint, typename TRealVector>
792 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
793 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
794 computeManifoldInnerEdges() const
795 {
796  Edges edges;
797  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
798  {
799  if ( myEdgeFaces[ e ].size() == 2 )
800  edges.push_back( e );
801  }
802  return edges;
803 }
804 
805 //-----------------------------------------------------------------------------
806 template <typename TRealPoint, typename TRealVector>
807 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
808 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
809 computeManifoldInnerConsistentEdges() const
810 {
811  Edges edges;
812  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
813  {
814  if ( ( myEdgeRightFaces[ e ].size() == 1 )
815  && ( myEdgeLeftFaces[ e ].size() == 1 ) )
816  edges.push_back( e );
817  }
818  return edges;
819 }
820 
821 //-----------------------------------------------------------------------------
822 template <typename TRealPoint, typename TRealVector>
823 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
824 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
825 computeManifoldInnerUnconsistentEdges() const
826 {
827  Edges edges;
828  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
829  {
830  if ( ( myEdgeRightFaces[ e ].size() == 2 )
831  || ( myEdgeLeftFaces[ e ].size() == 2 ) )
832  edges.push_back( e );
833  }
834  return edges;
835 }
836 
837 //-----------------------------------------------------------------------------
838 template <typename TRealPoint, typename TRealVector>
839 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
840 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
841 computeNonManifoldEdges() const
842 {
843  Edges edges;
844  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
845  {
846  if ( myEdgeFaces[ e ].size() >= 3 )
847  edges.push_back( e );
848  }
849  return edges;
850 }
851 
852 
853 //-----------------------------------------------------------------------------
854 template <typename TRealPoint, typename TRealVector>
855 bool
856 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
857 isFlippable( const Edge e ) const
858 {
859  /// An edge is (topologically) flippable iff: (1) it does not lie
860  /// on the boundary, (2) it is bordered by two triangles, one that
861  /// to its right, one to its left, (3) the two other vertices of
862  /// the quad are not already neighbors.
863 
864  // edge is (i,j) with i<j
865 
866  // (1) the edge must be bordered by two faces, one on its left, one on its right.
867  const Faces& rfaces = edgeRightFaces( e );
868  if ( rfaces.size() != 1 ) return false; //< not one face to the right
869  const Faces& lfaces = edgeLeftFaces ( e );
870  if ( lfaces.size() != 1 ) return false; //< not one face to the left
871 
872  // (2) both faces must be triangles
873  const Face rf = rfaces.front(); //< some `(..., j, i, ... )` since faces are ccw.
874  const Face lf = lfaces.front(); //< some `(..., i, j, ... )` since faces are ccw.
875  const Vertices& rvtx = incidentVertices( rf );
876  if ( rvtx.size() != 3 ) return false; //< right face is not a triangle
877  const Vertices& lvtx = incidentVertices( lf );
878  if ( lvtx.size() != 3 ) return false; //< left face is not a triangle
879 
880  // (3) the two other vertices of the quad are not already neighbors.
881  Vertex i, j;
882  std::tie( i, j ) = edgeVertices( e );
883  const auto ir = ( rvtx[ 0 ] == i ) ? 0 : ( ( rvtx[ 1 ] == i ) ? 1 : 2 );
884  const auto il = ( lvtx[ 0 ] == i ) ? 0 : ( ( lvtx[ 1 ] == i ) ? 1 : 2 );
885  const Vertex k = rvtx[ (ir + 1) % 3 ]; // right triangle is (j,i,k)
886  const Vertex l = lvtx[ (il + 2) % 3 ]; // left triangle is (i,j,l).
887  if ( ! ( k != i && k != j && l != i && l != j && k != l ) )
888  {
889  /// May happen if we have identical triangles (i,j,k) and (j,i,k)
890  if ( k == l ) return false;
891  /// Otherwise there is a problem.
892  trace.error() << "[SurfaceMesh::isFlippable] Invalid neighbors to edge: "
893  << "(i,j,k,l) == (" << i << "," << j << "," << k << "," << l << ")"
894  << std::endl
895  << "right=(" << rvtx[ 0 ] << "," << rvtx[ 1 ] << "," << rvtx[ 2 ] << ")" << std::endl
896  << "left =(" << lvtx[ 0 ] << "," << lvtx[ 1 ] << "," << lvtx[ 2 ] << ")" << std::endl;
897  return false;
898  }
899  const Vertices& Nk = neighborVertices( k );
900  const auto itl = std::find( Nk.cbegin(), Nk.cend(), l );
901  return itl == Nk.cend();
902 }
903 
904 //-----------------------------------------------------------------------------
905 template <typename TRealPoint, typename TRealVector>
906 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::VertexPair
907 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
908 otherDiagonal( const Edge e ) const
909 {
910  // only valid if `isFlippable( e )` is true.
911  const Faces& rfaces = edgeRightFaces( e );
912  const Faces& lfaces = edgeLeftFaces ( e );
913  const Face rf = rfaces.front(); //< some `(..., j, i, ... )` since faces are ccw.
914  const Face lf = lfaces.front(); //< some `(..., i, j, ... )` since faces are ccw.
915  const Vertices& rvtx = incidentVertices( rf );
916  const Vertices& lvtx = incidentVertices( lf );
917  Vertex i, j;
918  std::tie( i, j ) = edgeVertices( e );
919  const auto ir = ( rvtx[ 0 ] == i ) ? 0 : ( ( rvtx[ 1 ] == i ) ? 1 : 2 );
920  const auto il = ( lvtx[ 0 ] == i ) ? 0 : ( ( lvtx[ 1 ] == i ) ? 1 : 2 );
921  const Vertex k = rvtx[ (ir + 1) % 3 ]; // right triangle is (j,i,k)
922  const Vertex l = lvtx[ (il + 2) % 3 ]; // left triangle is (i,j,l).
923  return ( k < l ) ? std::make_pair( k, l ) : std::make_pair( l, k );
924 }
925 
926 
927 //-----------------------------------------------------------------------------
928 template <typename TRealPoint, typename TRealVector>
929 void
930 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
931 flip( const Edge e, bool recompute_face_normals )
932 {
933  /*
934  l l
935  / \ /|\
936  / \ / | \
937  / \ / | \
938  / lf \ / | \
939  / \ / | \
940  i --- e --- j ==> i lf e rf j if k < l otherwise rf and lf are swapped
941  \ / \ | /
942  \ rf / \ | /
943  \ / \ | /
944  \ / \ | /
945  \ / \|/
946  k k
947  */
948 
949  // (1) We must collect all information: right and left face, vertices k and l
950  const Face rf = edgeRightFaces( e ).front(); //< some `(..., j, i, ... )` since faces are ccw.
951  const Face lf = edgeLeftFaces ( e ).front(); //< some `(..., i, j, ... )` since faces are ccw.
952  Vertices& rvtx = myIncidentVertices[ rf ];
953  Vertices& lvtx = myIncidentVertices[ lf ];
954  Vertex i, j;
955  std::tie( i, j ) = edgeVertices( e );
956  const auto ir = ( rvtx[ 0 ] == i ) ? 0 : ( ( rvtx[ 1 ] == i ) ? 1 : 2 );
957  const auto il = ( lvtx[ 0 ] == i ) ? 0 : ( ( lvtx[ 1 ] == i ) ? 1 : 2 );
958  const Vertex k = rvtx[ (ir + 1) % 3 ]; // right triangle is (j,i,k)
959  const Vertex l = lvtx[ (il + 2) % 3 ]; // left triangle is (i,j,l).
960 
961  // (2) we must update all arrays.
962  // std::vector< Faces > myNeighborFaces; //< not done
963  if ( k < l )
964  {
965  /*
966  l l
967  / \ /|\
968  / \ / | \
969  / \ / | \
970  / lf \ / | \
971  / \ / | \
972  i --- e --- j ==> i lf e rf j ( k < l )
973  \ / \ | /
974  \ rf / \ | /
975  \ / \ | /
976  \ / \ | /
977  \ / \|/
978  k k
979  */
980  // e=(k,l) with rf as right face et lf as left face
981  rvtx[ 0 ] = l; rvtx[ 1 ] = k; rvtx[ 2 ] = j;
982  lvtx[ 0 ] = k; lvtx[ 1 ] = l; lvtx[ 2 ] = i;
983  VertexPair kl = std::make_pair( k, l );
984  myEdgeVertices [ e ] = kl;
985  myVertexPairEdge[ kl ] = e;
986  removeIndex ( myIncidentFaces[ i ], rf );
987  removeIndex ( myIncidentFaces[ j ], lf );
988  addIndex ( myIncidentFaces[ k ], lf );
989  addIndex ( myIncidentFaces[ l ], rf );
990  removeIndex ( myNeighborVertices[ i ], j );
991  removeIndex ( myNeighborVertices[ j ], i );
992  addIndex ( myNeighborVertices[ k ], l );
993  addIndex ( myNeighborVertices[ l ], k );
994  // No need to update myEdgeFaces, myEdgeRightFaces and myEdgeLeftFaces for edge e.
995  const auto e_ik = makeEdge( i, k );
996  const bool e_ik_left = i < k;
997  replaceIndex( myEdgeFaces[ e_ik ], rf, lf );
998  if ( e_ik_left ) myEdgeLeftFaces [ e_ik ][ 0 ] = lf;
999  else myEdgeRightFaces[ e_ik ][ 0 ] = lf;
1000  // nothing to change for e_kj (rf is still the incident face)
1001  const auto e_jl = makeEdge( j, l );
1002  const bool e_jl_left = j < l;
1003  replaceIndex( myEdgeFaces[ e_jl ], lf, rf );
1004  if ( e_jl_left ) myEdgeLeftFaces [ e_jl ][ 0 ] = rf;
1005  else myEdgeRightFaces[ e_jl ][ 0 ] = rf;
1006  // nothing to change for e_li (lf is still the incident face)
1007  // vertex normals are not updated.
1008  // face normals are recomputed if asked for.
1009  if ( recompute_face_normals )
1010  {
1011  computeFaceNormalFromPositions( rf );
1012  computeFaceNormalFromPositions( lf );
1013  }
1014  }
1015  else // k > l
1016  {
1017  /*
1018  l l
1019  / \ /|\
1020  / \ / | \
1021  / \ / | \
1022  / lf \ / | \
1023  / \ / | \
1024  i --- e --- j ==> i rf e lf j (k > l)
1025  \ / \ | /
1026  \ rf / \ | /
1027  \ / \ | /
1028  \ / \ | /
1029  \ / \|/
1030  k k
1031  */
1032  // e=(l,k) with rf as right face et lf as left face
1033  rvtx[ 0 ] = i; rvtx[ 1 ] = k; rvtx[ 2 ] = l;
1034  lvtx[ 0 ] = j; lvtx[ 1 ] = l; lvtx[ 2 ] = k;
1035  VertexPair lk = std::make_pair( l, k );
1036  myEdgeVertices [ e ] = lk;
1037  myVertexPairEdge[ lk ] = e;
1038  removeIndex ( myIncidentFaces[ i ], lf );
1039  removeIndex ( myIncidentFaces[ j ], rf );
1040  addIndex ( myIncidentFaces[ k ], lf );
1041  addIndex ( myIncidentFaces[ l ], rf );
1042  removeIndex ( myNeighborVertices[ i ], j );
1043  removeIndex ( myNeighborVertices[ j ], i );
1044  addIndex ( myNeighborVertices[ k ], l );
1045  addIndex ( myNeighborVertices[ l ], k );
1046  // No need to update myEdgeFaces, myEdgeRightFaces and myEdgeLeftFaces for edge e.
1047  const auto e_kj = makeEdge( k, j );
1048  const bool e_kj_left = k < j;
1049  replaceIndex( myEdgeFaces[ e_kj ], rf, lf );
1050  if ( e_kj_left ) myEdgeLeftFaces [ e_kj ][ 0 ] = lf;
1051  else myEdgeRightFaces[ e_kj ][ 0 ] = lf;
1052  // nothing to change for e_jl (lf is still the incident face)
1053  const auto e_li = makeEdge( l, i );
1054  const bool e_li_left = l < i;
1055  replaceIndex( myEdgeFaces[ e_li ], lf, rf );
1056  if ( e_li_left ) myEdgeLeftFaces [ e_li ][ 0 ] = rf;
1057  else myEdgeRightFaces[ e_li ][ 0 ] = rf;
1058  // nothing to change for e_ik (rf is still the incident face)
1059  // vertex normals are not updated.
1060  // face normals are recomputed if asked for.
1061  if ( recompute_face_normals )
1062  {
1063  computeFaceNormalFromPositions( rf );
1064  computeFaceNormalFromPositions( lf );
1065  }
1066  }
1067 }
1068 
1069 ///////////////////////////////////////////////////////////////////////////////
1070 // Implementation of inline functions //
1071 
1072 //-----------------------------------------------------------------------------
1073 template <typename TRealPoint, typename TRealVector>
1074 std::ostream&
1075 DGtal::operator<< ( std::ostream & out,
1076  const SurfaceMesh<TRealPoint, TRealVector> & object )
1077 {
1078  object.selfDisplay( out );
1079  return out;
1080 }
1081 
1082 
1083 ///////////////////////////////////////////////////////////////////////////////
1084 ///////////////////////////////////////////////////////////////////////////////