DGtal  1.4.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  myEdgeFaces.clear();
108  myEdgeRightFaces.clear();
109  myEdgeLeftFaces.clear();
110 }
111 
112 //-----------------------------------------------------------------------------
113 template <typename TRealPoint, typename TRealVector>
114 template <typename RealVectorIterator>
115 bool
116 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
117 setVertexNormals( RealVectorIterator itN, RealVectorIterator itNEnd )
118 {
119  myVertexNormals = std::vector< RealVector >( itN, itNEnd );
120  return myVertexNormals.size() == myPositions.size();
121 }
122 
123 //-----------------------------------------------------------------------------
124 template <typename TRealPoint, typename TRealVector>
125 template <typename RealVectorIterator>
126 bool
127 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
128 setFaceNormals( RealVectorIterator itN, RealVectorIterator itNEnd )
129 {
130  myFaceNormals = std::vector< RealVector >( itN, itNEnd );
131  return myFaceNormals.size() == myIncidentVertices.size();
132 }
133 
134 //-----------------------------------------------------------------------------
135 template <typename TRealPoint, typename TRealVector>
136 void
137 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
138 computeFaceNormalsFromPositions()
139 {
140  myFaceNormals.resize( myIncidentVertices.size() );
141  Index f = 0;
142  for ( auto face : myIncidentVertices )
143  {
144  RealPoint p; // barycenter
145  RealVector n; // normal
146  // compute barycenter
147  for ( auto idx : face ) p += myPositions[ idx ];
148  p /= face.size();
149  // compute normal as sum of triangle normal vectors.
150  for ( Index i = 0; i < face.size(); ++i )
151  {
152  const Index j = face[ i ];
153  const Index nj = face[ (i+1) % face.size() ];
154  n += ( myPositions[ j ] - p ).crossProduct( myPositions[ nj ] - p );
155  }
156  auto n_norm = n.norm();
157  myFaceNormals[ f ] = n_norm != 0.0 ? n / n_norm : n;
158  f++;
159  }
160 }
161 
162 //-----------------------------------------------------------------------------
163 template <typename TRealPoint, typename TRealVector>
164 void
165 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
166 computeFaceNormalsFromVertexNormals()
167 {
168  if ( myVertexNormals.empty() ) return;
169  myFaceNormals.resize( myIncidentVertices.size() );
170  Index f = 0;
171  for ( auto face : myIncidentVertices )
172  {
173  RealVector n; // normal
174  for ( auto idx : face ) n += myVertexNormals[ idx ];
175  auto n_norm = n.norm();
176  myFaceNormals[ f ] = n_norm != 0.0 ? n / n_norm : n;
177  f++;
178  }
179 }
180 //-----------------------------------------------------------------------------
181 template <typename TRealPoint, typename TRealVector>
182 void
183 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
184 computeVertexNormalsFromFaceNormals()
185 {
186  if ( myFaceNormals.empty() ) return;
187  myVertexNormals.resize( myIncidentFaces.size() );
188  Index v = 0;
189  for ( auto vertex : myIncidentFaces )
190  {
191  RealVector n; // normal
192  for ( auto idx : vertex ) n += myFaceNormals[ idx ];
193  auto n_norm = n.norm();
194  myVertexNormals[ v ] = n_norm != 0.0 ? n / n_norm : n;
195  v++;
196  }
197 }
198 //-----------------------------------------------------------------------------
199 template <typename TRealPoint, typename TRealVector>
200 void
201 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
202 computeVertexNormalsFromFaceNormalsWithMaxWeights()
203 {
204  if ( myFaceNormals.empty() ) return;
205  myVertexNormals.resize( myIncidentFaces.size() );
206  Index v = 0;
207  for ( auto incident_faces : myIncidentFaces )
208  {
209  RealVector n; // normal
210  const auto weights = getMaxWeights( v );
211  Index i = 0;
212  for ( auto idx_f : incident_faces ) n += weights[ i++ ] * myFaceNormals[ idx_f ];
213  auto n_norm = n.norm();
214  myVertexNormals[ v ] = n_norm != 0.0 ? n / n_norm : n;
215  v++;
216  }
217 }
218 
219 //-----------------------------------------------------------------------------
220 template <typename TRealPoint, typename TRealVector>
221 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalars
222 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
223 getMaxWeights( Index v ) const
224 {
225  Scalars weights;
226  const auto & neighbors = myNeighborVertices[ v ];
227  const RealPoint x = myPositions[ v ];
228  for ( auto idx_f : myIncidentFaces[ v ] )
229  {
230  // Find adjacent vertices to v
231  std::vector< Index > adj_vertices;
232  for ( auto idx_v : myIncidentVertices[ idx_f ] )
233  {
234  auto it = std::find( neighbors.cbegin(), neighbors.cend(), idx_v );
235  if ( it != neighbors.cend() ) adj_vertices.push_back( *it );
236  }
237  if ( adj_vertices.size() != 2 )
238  {
239  trace.warning() << "[SurfaceMesh::getMaxWeights] "
240  << adj_vertices.size() << " adjacent vertices to vertex "
241  << v << " on face" << idx_f << "." << std::endl;
242  for ( auto a : adj_vertices ) std::cerr << " " << a;
243  std::cerr << std::endl;
244  }
245  if (adj_vertices.size() >= 2 )
246  {
247  const Scalar area = faceArea( idx_f );
248  const Scalar l1 = ( myPositions[ adj_vertices[ 0 ] ] - x ).squaredNorm();
249  const Scalar l2 = ( myPositions[ adj_vertices[ 1 ] ] - x ).squaredNorm();
250  const Scalar l1l2 = l1 * l2;
251  weights.push_back( l1l2 != 0 ? fabs( area ) / l1l2 : 0.0 );
252  }
253  else weights.push_back( 0.0 );
254  }
255  return weights;
256 }
257 
258 //-----------------------------------------------------------------------------
259 template <typename TRealPoint, typename TRealVector>
260 template <typename AnyRing>
261 std::vector<AnyRing>
262 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
263 computeFaceValuesFromVertexValues( const std::vector<AnyRing>& vvalues ) const
264 {
265  ASSERT( vvalues.size() == nbVertices() );
266  std::vector<AnyRing> fvalues( nbFaces() );
267  Index f = 0;
268  for ( auto face : myIncidentVertices )
269  {
270  AnyRing n = NumberTraits<AnyRing>::ZERO;
271  for ( auto idx : face ) n += vvalues[ idx ];
272  fvalues[ f++ ] = n / face.size();
273  }
274  return fvalues;
275 }
276 
277 template <typename TRealPoint, typename TRealVector>
278 template <typename AnyRing>
279 std::vector<AnyRing>
280 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
281 computeVertexValuesFromFaceValues( const std::vector<AnyRing>& fvalues ) const
282 {
283  ASSERT( fvalues.size() == nbFaces() );
284  std::vector<AnyRing> vvalues( nbVertices() );
285  Index v = 0;
286  for ( auto vertex : myIncidentFaces )
287  {
288  AnyRing n = NumberTraits<AnyRing>::ZERO;
289  for ( auto idx : vertex ) n += fvalues[ idx ];
290  vvalues[ v++ ] = n / vertex.size();
291  }
292  return vvalues;
293 }
294 
295 //-----------------------------------------------------------------------------
296 template <typename TRealPoint, typename TRealVector>
297 std::vector<TRealVector>
298 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
299 computeFaceUnitVectorsFromVertexUnitVectors
300 ( const std::vector<RealVector>& vuvectors ) const
301 {
302  ASSERT( vuvectors.size() == nbVertices() );
303  std::vector<RealVector> fuvectors( nbFaces() );
304  Index f = 0;
305  for ( auto face : myIncidentVertices )
306  {
307  RealVector n;
308  for ( auto idx : face ) n += vuvectors[ idx ];
309  const auto n_norm = n.norm();
310  fuvectors[ f++ ] = n_norm != 0.0 ? n / n_norm : n;
311  }
312  return fuvectors;
313 }
314 
315 //-----------------------------------------------------------------------------
316 template <typename TRealPoint, typename TRealVector>
317 std::vector<TRealVector>
318 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
319 computeVertexUnitVectorsFromFaceUnitVectors
320 ( const std::vector<RealVector>& fuvectors ) const
321 {
322  ASSERT( fuvectors.size() == nbFaces() );
323  std::vector<RealVector> vuvectors( nbVertices() );
324  Index v = 0;
325  for ( auto vertex : myIncidentFaces )
326  {
327  RealVector n;
328  for ( auto idx : vertex ) n += fuvectors[ idx ];
329  const auto n_norm = n.norm();
330  vuvectors[ v++ ] = n_norm != 0.0 ? n / n_norm : n;
331  }
332  return vuvectors;
333 }
334 
335 
336 //-----------------------------------------------------------------------------
337 template <typename TRealPoint, typename TRealVector>
338 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edge
339 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
340 makeEdge( Vertex i, Vertex j ) const
341 {
342  VertexPair vp = i < j ? std::make_pair( i,j ) : std::make_pair( j,i );
343  auto it = std::lower_bound( myEdgeVertices.cbegin(), myEdgeVertices.cend(), vp );
344  if ( it == myEdgeVertices.cend() || *it != vp ) return nbEdges();
345  return it - myEdgeVertices.cbegin();
346 }
347 
348 //-----------------------------------------------------------------------------
349 template <typename TRealPoint, typename TRealVector>
350 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
351 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
352 averageEdgeLength() const
353 {
354  double lengths = 0.0;
355  for ( Edge e = 0; e < nbEdges(); ++e )
356  {
357  auto vtcs = edgeVertices( e );
358  const RealPoint p = myPositions[ vtcs.first ];
359  const RealVector pq = myPositions[ vtcs.second ] - p;
360  lengths += pq.norm();
361  }
362  lengths /= nbEdges();
363  return lengths;
364 }
365 
366 //-----------------------------------------------------------------------------
367 template <typename TRealPoint, typename TRealVector>
368 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
369 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
370 localWindow( Face f ) const
371 {
372  const RealPoint x = faceCentroid( f );
373  Scalar local_length = 0.0;
374  for ( auto v : incidentVertices( f ) )
375  local_length += ( myPositions[ v ] - x ).norm();
376  return local_length / incidentVertices( f ).size();
377 }
378 
379 //-----------------------------------------------------------------------------
380 template <typename TRealPoint, typename TRealVector>
381 void
382 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
383 perturbateWithUniformRandomNoise( Scalar p )
384 {
385  for ( auto& x : myPositions )
386  {
387  RealVector d( rand01()*2.0 - 1.0, rand01()*2.0 - 1.0, rand01()*2.0 - 1.0 );
388  d = d.getNormalized();
389  Scalar l = rand01() * p;
390  x += l * d;
391  }
392 }
393 
394 //-----------------------------------------------------------------------------
395 template <typename TRealPoint, typename TRealVector>
396 void
397 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
398 perturbateWithAdaptiveUniformRandomNoise( Scalar p )
399 {
400  Scalars local_average_lengths( nbVertices() );
401  for ( Index v = 0; v < nbVertices(); ++v )
402  {
403  const RealPoint x = myPositions[ v ];
404  Scalar local_length = 0.0;
405  for ( auto nv : myNeighborVertices[ v ] )
406  local_length += ( myPositions[ nv ] - x ).norm();
407  local_average_lengths[ v ] = local_length / myNeighborVertices[ v ].size();
408  }
409  for ( Index v = 0; v < nbVertices(); ++v )
410  {
411  RealVector d( rand01()*2.0 - 1.0, rand01()*2.0 - 1.0, rand01()*2.0 - 1.0 );
412  d = d.getNormalized();
413  Scalar l = rand01() * p * local_average_lengths[ v ];
414  myPositions[ v ] += l * d;
415  }
416 }
417 
418 //-----------------------------------------------------------------------------
419 template <typename TRealPoint, typename TRealVector>
420 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::RealPoint
421 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
422 faceCentroid( Index f ) const
423 {
424  RealPoint c;
425  for ( auto v : myIncidentVertices[ f ] )
426  c += myPositions[ v ];
427  return c / myIncidentVertices[ f ].size();
428 }
429 
430 //-----------------------------------------------------------------------------
431 template <typename TRealPoint, typename TRealVector>
432 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::RealPoint
433 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
434 edgeCentroid( Index e ) const
435 {
436  const auto& vtcs = myEdgeVertices[ e ];
437  RealPoint c = myPositions[ vtcs.first ] + myPositions[ vtcs.second ];
438  return c / 2.0;
439 }
440 
441 //-----------------------------------------------------------------------------
442 template <typename TRealPoint, typename TRealVector>
443 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
444 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
445 faceArea( Index f ) const
446 {
447  Scalar area = 0.0;
448  const auto & inc_vtcs = myIncidentVertices[ f ];
449  RealPoint p = myPositions[ inc_vtcs.back() ];
450  const Index m = inc_vtcs.size() - 2;
451  for ( Index i = 0; i < m; ++i )
452  area += ( myPositions[ inc_vtcs[ i ] ] - p )
453  .crossProduct( myPositions[ inc_vtcs[ i+1 ] ] - p ).norm();
454  return area / 2.0;
455 }
456 
457 //-----------------------------------------------------------------------------
458 template <typename TRealPoint, typename TRealVector>
459 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces
460 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
461 computeFacesInclusionsInBall( Scalar r, Index f ) const
462 {
463  return computeFacesInclusionsInBall( r, f, faceCentroid( f ) );
464 }
465 
466 //-----------------------------------------------------------------------------
467 template <typename TRealPoint, typename TRealVector>
468 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces
469 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
470 computeFacesInclusionsInBall( Scalar r, Index f, RealPoint p ) const
471 {
472  WeightedFaces result;
473  if ( r < 0.000001 )
474  {
475  result.push_back( std::make_pair( f, 0.000001 ) );
476  return result;
477  }
478  std::unordered_set< Index > marked;
479  std::queue< Index > active;
480  active.push( f );
481  marked.insert( f );
482  while ( ! active.empty() )
483  {
484  Index current = active.front();
485  active.pop();
486  Scalar weight = faceInclusionRatio( p, r, current );
487  if ( weight > 0.0 )
488  {
489  result.push_back( std::make_pair( current, weight ) );
490  auto neighbors = myNeighborFaces[ current ];
491  for ( auto n : neighbors )
492  if ( marked.find( n ) == marked.end() )
493  {
494  active.push( n );
495  marked.insert( n );
496  }
497  }
498  }
499  return result;
500 }
501 
502 //-----------------------------------------------------------------------------
503 template <typename TRealPoint, typename TRealVector>
504 std::tuple
505 < typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Vertices,
506  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedEdges,
507  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces >
508 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
509 computeCellsInclusionsInBall( Scalar r, Index f ) const
510 {
511  return computeCellsInclusionsInBall( r, f, faceCentroid( f ) );
512 }
513 
514 //-----------------------------------------------------------------------------
515 template <typename TRealPoint, typename TRealVector>
516 std::tuple
517 < typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Vertices,
518  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedEdges,
519  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces >
520 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
521 computeCellsInclusionsInBall( Scalar r, Index f, RealPoint p ) const
522 {
523  Vertices result_v;
524  std::set<Vertex> set_v;
525  WeightedEdges result_e;
526  WeightedFaces result_f;
527  if ( r < 0.000001 )
528  {
529  result_f.push_back( std::make_pair( f, 0.000001 ) );
530  return std::make_tuple( result_v, result_e, result_f );
531  }
532  std::set< Index > marked;
533  std::queue< Index > active;
534  active.push( f );
535  marked.insert( f );
536  while ( ! active.empty() )
537  {
538  Index current = active.front();
539  active.pop();
540  Scalar fweight = faceInclusionRatio( p, r, current );
541  if ( fweight > 0.0 )
542  {
543  result_f.push_back( std::make_pair( current, fweight ) );
544  // Taking care of faces, and the breadth-first traversal
545  const auto& neighbors = myNeighborFaces[ current ];
546  for ( auto n : neighbors )
547  if ( marked.find( n ) == marked.end() )
548  {
549  active.push( n );
550  marked.insert( n );
551  }
552  // Taking care of edges and vertices
553  const auto& inc_v = myIncidentVertices[ current ];
554  for ( Size i = 0; i < inc_v.size(); ++i )
555  {
556  const Vertex vi = inc_v[ i ];
557  const Vertex vn = inc_v[ (i+1) % inc_v.size() ];
558  if ( vertexInclusionRatio( p, r, vi ) > 0.0 )
559  set_v.insert( vi );
560  if ( vn < vi ) continue; // edges are ordered pairs
561  const Edge e_ij = makeEdge( vi, vn );
562  if ( e_ij >= nbEdges() ) {
563  trace.error() << "bad edge " << vi << " " << vn << std::endl;
564  continue;
565  }
566  Scalar eweight = edgeInclusionRatio( p, r, e_ij );
567  if ( eweight > 0.0 )
568  result_e.push_back( std::make_pair( e_ij, eweight ) );
569  }
570  }
571  }
572  result_v = Vertices( set_v.cbegin(), set_v.cend() );
573  return std::make_tuple( result_v, result_e, result_f );
574 }
575 
576 //-----------------------------------------------------------------------------
577 template <typename TRealPoint, typename TRealVector>
578 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
579 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
580 faceInclusionRatio( RealPoint p, Scalar r, Index f ) const
581 {
582  const auto vertices = myIncidentVertices[ f ];
583  const RealPoint b = faceCentroid( f );
584  Scalar d_min = ( b - p ).norm();
585  Scalar d_max = d_min;
586  for ( auto v : vertices )
587  {
588  Scalar d = ( myPositions[ v ] - p ).norm();
589  d_max = std::max( d_max, d );
590  d_min = std::min( d_min, d );
591  }
592  if ( d_max <= r ) return 1.0;
593  else if ( r <= d_min ) return 0.0;
594  return ( r - d_min ) / ( d_max - d_min );
595 }
596 
597 //-----------------------------------------------------------------------------
598 template <typename TRealPoint, typename TRealVector>
599 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
600 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
601 edgeInclusionRatio( RealPoint p, Scalar r, Index e ) const
602 {
603  const auto vertices = myEdgeVertices[ e ];
604  const RealPoint b = edgeCentroid( e );
605  const Scalar d0 = ( myPositions[ vertices.first ] - p ).norm();
606  const Scalar d1 = ( myPositions[ vertices.second ] - p ).norm();
607  Scalar d_min = ( b - p ).norm();
608  Scalar d_max = d_min;
609  d_max = std::max( d_max, std::max( d0, d1 ) );
610  d_min = std::min( d_min, std::min( d0, d1 ) );
611  if ( d_max <= r ) return 1.0;
612  else if ( r <= d_min ) return 0.0;
613  return ( r - d_min ) / ( d_max - d_min );
614 }
615 
616 //-----------------------------------------------------------------------------
617 template <typename TRealPoint, typename TRealVector>
618 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
619 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
620 vertexInclusionRatio( RealPoint p, Scalar r, Index v ) const
621 {
622  const RealPoint b = myPositions[ v ];
623  return ( ( b - p ).norm() <= r ) ? 1.0 : 0.0;
624 }
625 
626 
627 ///////////////////////////////////////////////////////////////////////////////
628 // Interface - public :
629 
630 //-----------------------------------------------------------------------------
631 template <typename TRealPoint, typename TRealVector>
632 void
633 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
634 selfDisplay ( std::ostream & out ) const
635 {
636  out << "[SurfaceMesh" << ( isValid() ? " (OK)" : " (KO)" )
637  << " #V=" << myPositions.size()
638  << " #VN=" << myVertexNormals.size()
639  << " #E=" << myEdgeVertices.size()
640  << " #F=" << myIncidentVertices.size()
641  << " #FN=" << myFaceNormals.size();
642  double nb_nf = 0.0;
643  double nb_nv = 0.0;
644  double nb_nfe = 0.0;
645  for ( auto nf : myNeighborFaces ) nb_nf += nf.size();
646  for ( auto nv : myNeighborVertices ) nb_nv += nv.size();
647  for ( auto nfe : myEdgeFaces ) nb_nfe += nfe.size();
648  nb_nf /= nbFaces();
649  nb_nv /= nbVertices();
650  nb_nfe /= nbEdges();
651  out << " E[IF]=" << nb_nf << " E[IV]=" << nb_nv << " E[IFE]=" << nb_nfe;
652  out << "]";
653 }
654 
655 //-----------------------------------------------------------------------------
656 template <typename TRealPoint, typename TRealVector>
657 bool
658 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
659 isValid() const
660 {
661  return myPositions.size() == myIncidentFaces.size()
662  && ( myVertexNormals.size() == 0
663  || ( myVertexNormals.size() == myPositions.size() ) )
664  && ( myFaceNormals.size() == 0
665  || ( myFaceNormals.size() == myIncidentVertices.size() ) );
666 }
667 
668 
669 //-----------------------------------------------------------------------------
670 template <typename TRealPoint, typename TRealVector>
671 void
672 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
673 computeNeighbors()
674 {
675  myNeighborFaces .resize( nbFaces() );
676  myNeighborVertices.resize( nbVertices() );
677  std::vector< std::set< Index > > tmp( nbVertices() );
678  // For each vertex, computes its neighboring vertices
679  for ( auto incident_vertices : allIncidentVertices() )
680  {
681  const Size nb_iv = incident_vertices.size();
682  for ( Size k = 0; k < nb_iv; ++k )
683  {
684  tmp[ incident_vertices[ k ] ].insert( incident_vertices[ (k+1)%nb_iv ] );
685  tmp[ incident_vertices[ (k+1)%nb_iv ] ].insert( incident_vertices[ k ] );
686  }
687  }
688  // For each vertex, computes its neighboring vertices
689  for ( Index idx_v = 0; idx_v < nbVertices(); ++idx_v )
690  myNeighborVertices[ idx_v ] = Vertices( tmp[ idx_v ].cbegin(), tmp[ idx_v ].cend() );
691 
692  // For each face, computes its neighboring faces
693  Index idx_f = 0;
694  for ( auto incident_vertices : allIncidentVertices() )
695  {
696  std::set< Index > neighbor_faces_set;
697  std::sort( incident_vertices.begin(), incident_vertices.end() );
698  for ( auto idx_v : incident_vertices )
699  {
700  const auto & incident_faces = incidentFaces( idx_v );
701  for ( auto inc_f : incident_faces )
702  {
703  // Keep only faces incident to two vertices of f.
704  auto incident_vertices2 = incidentVertices( inc_f );
705  std::sort( incident_vertices2.begin(), incident_vertices2.end() );
706  Vertices common;
707  std::set_intersection( incident_vertices.cbegin(), incident_vertices.cend(),
708  incident_vertices2.cbegin(), incident_vertices2.cend(),
709  std::back_inserter( common ) );
710  if ( common.size() == 2 )
711  neighbor_faces_set.insert( inc_f );
712  }
713  }
714  neighbor_faces_set.erase( idx_f );
715  Faces neighbor_faces( neighbor_faces_set.begin(), neighbor_faces_set.end() );
716  myNeighborFaces[ idx_f++ ] = neighbor_faces;
717  }
718 }
719 
720 //-----------------------------------------------------------------------------
721 template <typename TRealPoint, typename TRealVector>
722 void
723 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
724 computeEdges()
725 {
726  std::map< VertexPair, std::vector<Face> > edge2face_right;
727  std::map< VertexPair, std::vector<Face> > edge2face_left;
728  std::set< VertexPair > edges;
729  Index idx_f = 0;
730  for ( auto incident_vertices : allIncidentVertices() )
731  {
732  const Size n = incident_vertices.size();
733  for ( Size i = 0; i < n; i++ )
734  {
735  VertexPair e = std::make_pair( incident_vertices[ i ],
736  incident_vertices[ (i+1) % n ] );
737  if ( e.first < e.second ) {
738  edge2face_left[ e ].push_back( idx_f );
739  } else {
740  std::swap( e.first, e.second );
741  edge2face_right[ e ].push_back( idx_f );
742  }
743  edges.insert( e );
744  }
745  idx_f++;
746  }
747  const Size nbe = edges.size();
748  myEdgeVertices.resize ( nbe );
749  myEdgeFaces.resize ( nbe );
750  myEdgeRightFaces.resize( nbe );
751  myEdgeLeftFaces.resize ( nbe );
752  Index idx_e = 0;
753  for ( auto e : edges )
754  {
755  myEdgeVertices [ idx_e ] = e;
756  myEdgeLeftFaces [ idx_e ] = edge2face_left [ e ];
757  myEdgeRightFaces[ idx_e ] = edge2face_right[ e ];
758  myEdgeFaces [ idx_e ] = myEdgeRightFaces[ idx_e ];
759  myEdgeFaces [ idx_e ].insert( myEdgeFaces[ idx_e ].end(),
760  myEdgeLeftFaces[ idx_e ].cbegin(),
761  myEdgeLeftFaces[ idx_e ].cend() );
762  idx_e++;
763  }
764 }
765 
766 //-----------------------------------------------------------------------------
767 template <typename TRealPoint, typename TRealVector>
768 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
769 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
770 computeManifoldBoundaryEdges() const
771 {
772  Edges edges;
773  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
774  {
775  if ( myEdgeFaces[ e ].size() == 1 )
776  edges.push_back( e );
777  }
778  return edges;
779 }
780 
781 //-----------------------------------------------------------------------------
782 template <typename TRealPoint, typename TRealVector>
783 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
784 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
785 computeManifoldInnerEdges() const
786 {
787  Edges edges;
788  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
789  {
790  if ( myEdgeFaces[ e ].size() == 2 )
791  edges.push_back( e );
792  }
793  return edges;
794 }
795 
796 //-----------------------------------------------------------------------------
797 template <typename TRealPoint, typename TRealVector>
798 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
799 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
800 computeManifoldInnerConsistentEdges() const
801 {
802  Edges edges;
803  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
804  {
805  if ( ( myEdgeRightFaces[ e ].size() == 1 )
806  && ( myEdgeLeftFaces[ e ].size() == 1 ) )
807  edges.push_back( e );
808  }
809  return edges;
810 }
811 
812 //-----------------------------------------------------------------------------
813 template <typename TRealPoint, typename TRealVector>
814 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
815 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
816 computeManifoldInnerUnconsistentEdges() const
817 {
818  Edges edges;
819  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
820  {
821  if ( ( myEdgeRightFaces[ e ].size() == 2 )
822  || ( myEdgeLeftFaces[ e ].size() == 2 ) )
823  edges.push_back( e );
824  }
825  return edges;
826 }
827 
828 //-----------------------------------------------------------------------------
829 template <typename TRealPoint, typename TRealVector>
830 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
831 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
832 computeNonManifoldEdges() const
833 {
834  Edges edges;
835  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
836  {
837  if ( myEdgeFaces[ e ].size() >= 3 )
838  edges.push_back( e );
839  }
840  return edges;
841 }
842 
843 ///////////////////////////////////////////////////////////////////////////////
844 // Implementation of inline functions //
845 
846 //-----------------------------------------------------------------------------
847 template <typename TRealPoint, typename TRealVector>
848 std::ostream&
849 DGtal::operator<< ( std::ostream & out,
850  const SurfaceMesh<TRealPoint, TRealVector> & object )
851 {
852  object.selfDisplay( out );
853  return out;
854 }
855 
856 
857 ///////////////////////////////////////////////////////////////////////////////
858 ///////////////////////////////////////////////////////////////////////////////