DGtal  1.2.0
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  WeightedFaces result;
464  if ( r < 0.000001 )
465  {
466  result.push_back( std::make_pair( f, 0.000001 ) );
467  return result;
468  }
469  RealPoint p = faceCentroid( f );
470  std::unordered_set< Index > marked;
471  std::queue< Index > active;
472  active.push( f );
473  marked.insert( f );
474  while ( ! active.empty() )
475  {
476  Index current = active.front();
477  active.pop();
478  Scalar weight = faceInclusionRatio( p, r, current );
479  if ( weight > 0.0 )
480  {
481  result.push_back( std::make_pair( current, weight ) );
482  auto neighbors = myNeighborFaces[ current ];
483  for ( auto n : neighbors )
484  if ( marked.find( n ) == marked.end() )
485  {
486  active.push( n );
487  marked.insert( n );
488  }
489  }
490  }
491  return result;
492 }
493 
494 //-----------------------------------------------------------------------------
495 template <typename TRealPoint, typename TRealVector>
496 std::tuple
497 < typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Vertices,
498  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedEdges,
499  typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces >
500 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
501 computeCellsInclusionsInBall( Scalar r, Index f ) const
502 {
503  Vertices result_v;
504  std::set<Vertex> set_v;
505  WeightedEdges result_e;
506  WeightedFaces result_f;
507  if ( r < 0.000001 )
508  {
509  result_f.push_back( std::make_pair( f, 0.000001 ) );
510  return std::make_tuple( result_v, result_e, result_f );
511  }
512  RealPoint p = faceCentroid( f );
513  std::set< Index > marked;
514  std::queue< Index > active;
515  active.push( f );
516  marked.insert( f );
517  while ( ! active.empty() )
518  {
519  Index current = active.front();
520  active.pop();
521  Scalar fweight = faceInclusionRatio( p, r, current );
522  if ( fweight > 0.0 )
523  {
524  result_f.push_back( std::make_pair( current, fweight ) );
525  // Taking care of faces, and the breadth-first traversal
526  const auto& neighbors = myNeighborFaces[ current ];
527  for ( auto n : neighbors )
528  if ( marked.find( n ) == marked.end() )
529  {
530  active.push( n );
531  marked.insert( n );
532  }
533  // Taking care of edges and vertices
534  const auto& inc_v = myIncidentVertices[ current ];
535  for ( Size i = 0; i < inc_v.size(); ++i )
536  {
537  const Vertex vi = inc_v[ i ];
538  const Vertex vn = inc_v[ (i+1) % inc_v.size() ];
539  if ( vertexInclusionRatio( p, r, vi ) > 0.0 )
540  set_v.insert( vi );
541  const Edge e_ij = makeEdge( vi, vn );
542  if ( e_ij >= nbEdges() ) {
543  trace.error() << "bad edge " << vi << " " << vn << std::endl;
544  continue;
545  }
546  Scalar eweight = edgeInclusionRatio( p, r, e_ij );
547  if ( eweight > 0.0 )
548  result_e.push_back( std::make_pair( e_ij, eweight ) );
549  }
550  }
551  }
552  result_v = Vertices( set_v.cbegin(), set_v.cend() );
553  return std::make_tuple( result_v, result_e, result_f );
554 }
555 
556 //-----------------------------------------------------------------------------
557 template <typename TRealPoint, typename TRealVector>
558 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
559 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
560 faceInclusionRatio( RealPoint p, Scalar r, Index f ) const
561 {
562  const auto vertices = myIncidentVertices[ f ];
563  const RealPoint b = faceCentroid( f );
564  Scalar d_min = ( b - p ).norm();
565  Scalar d_max = d_min;
566  for ( auto v : vertices )
567  {
568  Scalar d = ( myPositions[ v ] - p ).norm();
569  d_max = std::max( d_max, d );
570  d_min = std::min( d_min, d );
571  }
572  if ( d_max <= r ) return 1.0;
573  else if ( r <= d_min ) return 0.0;
574  return ( r - d_min ) / ( d_max - d_min );
575 }
576 
577 //-----------------------------------------------------------------------------
578 template <typename TRealPoint, typename TRealVector>
579 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
580 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
581 edgeInclusionRatio( RealPoint p, Scalar r, Index e ) const
582 {
583  const auto vertices = myEdgeVertices[ e ];
584  const RealPoint b = edgeCentroid( e );
585  const Scalar d0 = ( myPositions[ vertices.first ] - p ).norm();
586  const Scalar d1 = ( myPositions[ vertices.second ] - p ).norm();
587  Scalar d_min = ( b - p ).norm();
588  Scalar d_max = d_min;
589  d_max = std::max( d_max, std::max( d0, d1 ) );
590  d_min = std::min( d_min, std::min( d0, d1 ) );
591  if ( d_max <= r ) return 1.0;
592  else if ( r <= d_min ) return 0.0;
593  return ( r - d_min ) / ( d_max - d_min );
594 }
595 
596 //-----------------------------------------------------------------------------
597 template <typename TRealPoint, typename TRealVector>
598 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
599 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
600 vertexInclusionRatio( RealPoint p, Scalar r, Index v ) const
601 {
602  const RealPoint b = myPositions[ v ];
603  return ( ( b - p ).norm() <= r ) ? 1.0 : 0.0;
604 }
605 
606 
607 ///////////////////////////////////////////////////////////////////////////////
608 // Interface - public :
609 
610 //-----------------------------------------------------------------------------
611 template <typename TRealPoint, typename TRealVector>
612 void
613 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
614 selfDisplay ( std::ostream & out ) const
615 {
616  out << "[SurfaceMesh" << ( isValid() ? " (OK)" : " (KO)" )
617  << " #V=" << myPositions.size()
618  << " #VN=" << myVertexNormals.size()
619  << " #E=" << myEdgeVertices.size()
620  << " #F=" << myIncidentVertices.size()
621  << " #FN=" << myFaceNormals.size();
622  double nb_nf = 0.0;
623  double nb_nv = 0.0;
624  double nb_nfe = 0.0;
625  for ( auto nf : myNeighborFaces ) nb_nf += nf.size();
626  for ( auto nv : myNeighborVertices ) nb_nv += nv.size();
627  for ( auto nfe : myEdgeFaces ) nb_nfe += nfe.size();
628  nb_nf /= nbFaces();
629  nb_nv /= nbVertices();
630  nb_nfe /= nbEdges();
631  out << " E[IF]=" << nb_nf << " E[IV]=" << nb_nv << " E[IFE]=" << nb_nfe;
632  out << "]";
633 }
634 
635 //-----------------------------------------------------------------------------
636 template <typename TRealPoint, typename TRealVector>
637 bool
638 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
639 isValid() const
640 {
641  return myPositions.size() == myIncidentFaces.size()
642  && ( myVertexNormals.size() == 0
643  || ( myVertexNormals.size() == myPositions.size() ) )
644  && ( myFaceNormals.size() == 0
645  || ( myFaceNormals.size() == myIncidentVertices.size() ) );
646 }
647 
648 
649 //-----------------------------------------------------------------------------
650 template <typename TRealPoint, typename TRealVector>
651 void
652 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
653 computeNeighbors()
654 {
655  myNeighborFaces .resize( nbFaces() );
656  myNeighborVertices.resize( nbVertices() );
657  std::vector< std::set< Index > > tmp( nbVertices() );
658  // For each vertex, computes its neighboring vertices
659  for ( auto incident_vertices : allIncidentVertices() )
660  {
661  const Size nb_iv = incident_vertices.size();
662  for ( Size k = 0; k < nb_iv; ++k )
663  {
664  tmp[ incident_vertices[ k ] ].insert( incident_vertices[ (k+1)%nb_iv ] );
665  tmp[ incident_vertices[ (k+1)%nb_iv ] ].insert( incident_vertices[ k ] );
666  }
667  }
668  // For each vertex, computes its neighboring vertices
669  for ( Index idx_v = 0; idx_v < nbVertices(); ++idx_v )
670  myNeighborVertices[ idx_v ] = Vertices( tmp[ idx_v ].cbegin(), tmp[ idx_v ].cend() );
671 
672  // For each face, computes its neighboring faces
673  Index idx_f = 0;
674  for ( auto incident_vertices : allIncidentVertices() )
675  {
676  std::set< Index > neighbor_faces_set;
677  std::sort( incident_vertices.begin(), incident_vertices.end() );
678  for ( auto idx_v : incident_vertices )
679  {
680  const auto & incident_faces = incidentFaces( idx_v );
681  for ( auto inc_f : incident_faces )
682  {
683  // Keep only faces incident to two vertices of f.
684  auto incident_vertices2 = incidentVertices( inc_f );
685  std::sort( incident_vertices2.begin(), incident_vertices2.end() );
686  Vertices common;
687  std::set_intersection( incident_vertices.cbegin(), incident_vertices.cend(),
688  incident_vertices2.cbegin(), incident_vertices2.cend(),
689  std::back_inserter( common ) );
690  if ( common.size() == 2 )
691  neighbor_faces_set.insert( inc_f );
692  }
693  }
694  neighbor_faces_set.erase( idx_f );
695  Faces neighbor_faces( neighbor_faces_set.begin(), neighbor_faces_set.end() );
696  myNeighborFaces[ idx_f++ ] = neighbor_faces;
697  }
698 }
699 
700 //-----------------------------------------------------------------------------
701 template <typename TRealPoint, typename TRealVector>
702 void
703 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
704 computeEdges()
705 {
706  std::map< VertexPair, std::vector<Face> > edge2face_right;
707  std::map< VertexPair, std::vector<Face> > edge2face_left;
708  std::set< VertexPair > edges;
709  Index idx_f = 0;
710  for ( auto incident_vertices : allIncidentVertices() )
711  {
712  const Size n = incident_vertices.size();
713  for ( Size i = 0; i < n; i++ )
714  {
715  VertexPair e = std::make_pair( incident_vertices[ i ],
716  incident_vertices[ (i+1) % n ] );
717  if ( e.first < e.second ) {
718  edge2face_left[ e ].push_back( idx_f );
719  } else {
720  std::swap( e.first, e.second );
721  edge2face_right[ e ].push_back( idx_f );
722  }
723  edges.insert( e );
724  }
725  idx_f++;
726  }
727  const Size nbe = edges.size();
728  myEdgeVertices.resize ( nbe );
729  myEdgeFaces.resize ( nbe );
730  myEdgeRightFaces.resize( nbe );
731  myEdgeLeftFaces.resize ( nbe );
732  Index idx_e = 0;
733  for ( auto e : edges )
734  {
735  myEdgeVertices [ idx_e ] = e;
736  myEdgeLeftFaces [ idx_e ] = edge2face_left [ e ];
737  myEdgeRightFaces[ idx_e ] = edge2face_right[ e ];
738  myEdgeFaces [ idx_e ] = myEdgeRightFaces[ idx_e ];
739  myEdgeFaces [ idx_e ].insert( myEdgeFaces[ idx_e ].end(),
740  myEdgeLeftFaces[ idx_e ].cbegin(),
741  myEdgeLeftFaces[ idx_e ].cend() );
742  idx_e++;
743  }
744 }
745 
746 //-----------------------------------------------------------------------------
747 template <typename TRealPoint, typename TRealVector>
748 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
749 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
750 computeManifoldBoundaryEdges() const
751 {
752  Edges edges;
753  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
754  {
755  if ( myEdgeFaces[ e ].size() == 1 )
756  edges.push_back( e );
757  }
758  return edges;
759 }
760 
761 //-----------------------------------------------------------------------------
762 template <typename TRealPoint, typename TRealVector>
763 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
764 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
765 computeManifoldInnerEdges() const
766 {
767  Edges edges;
768  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
769  {
770  if ( myEdgeFaces[ e ].size() == 2 )
771  edges.push_back( e );
772  }
773  return edges;
774 }
775 
776 //-----------------------------------------------------------------------------
777 template <typename TRealPoint, typename TRealVector>
778 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
779 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
780 computeManifoldInnerConsistentEdges() const
781 {
782  Edges edges;
783  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
784  {
785  if ( ( myEdgeRightFaces[ e ].size() == 1 )
786  && ( myEdgeLeftFaces[ e ].size() == 1 ) )
787  edges.push_back( e );
788  }
789  return edges;
790 }
791 
792 //-----------------------------------------------------------------------------
793 template <typename TRealPoint, typename TRealVector>
794 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
795 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
796 computeManifoldInnerUnconsistentEdges() const
797 {
798  Edges edges;
799  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
800  {
801  if ( ( myEdgeRightFaces[ e ].size() == 2 )
802  || ( myEdgeLeftFaces[ e ].size() == 2 ) )
803  edges.push_back( e );
804  }
805  return edges;
806 }
807 
808 //-----------------------------------------------------------------------------
809 template <typename TRealPoint, typename TRealVector>
810 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
811 DGtal::SurfaceMesh<TRealPoint, TRealVector>::
812 computeNonManifoldEdges() const
813 {
814  Edges edges;
815  for ( Index e = 0; e < myEdgeFaces.size(); ++e )
816  {
817  if ( myEdgeFaces[ e ].size() >= 3 )
818  edges.push_back( e );
819  }
820  return edges;
821 }
822 
823 ///////////////////////////////////////////////////////////////////////////////
824 // Implementation of inline functions //
825 
826 //-----------------------------------------------------------------------------
827 template <typename TRealPoint, typename TRealVector>
828 std::ostream&
829 DGtal::operator<< ( std::ostream & out,
830  const SurfaceMesh<TRealPoint, TRealVector> & object )
831 {
832  object.selfDisplay( out );
833  return out;
834 }
835 
836 
837 ///////////////////////////////////////////////////////////////////////////////
838 ///////////////////////////////////////////////////////////////////////////////
Generated on Tue Jun 1 2021 17:05:44 for DGtal by doxygen 1.9.1