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.
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.
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/>.
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
24 * Implementation of inline methods defined in SurfaceMesh.h
26 * This file is part of the DGtal library.
30//////////////////////////////////////////////////////////////////////////////
33//////////////////////////////////////////////////////////////////////////////
35///////////////////////////////////////////////////////////////////////////////
36// IMPLEMENTATION of inline methods.
37///////////////////////////////////////////////////////////////////////////////
39///////////////////////////////////////////////////////////////////////////////
40// ----------------------- Standard services ------------------------------
42//-----------------------------------------------------------------------------
43template <typename TRealPoint, typename TRealVector>
44template <typename RealPointIterator, typename VerticesIterator>
45DGtal::SurfaceMesh<TRealPoint, TRealVector>::
46SurfaceMesh( RealPointIterator itPos, RealPointIterator itPosEnd,
47 VerticesIterator itVertices, VerticesIterator itVerticesEnd )
49 bool ok = init( itPos, itPosEnd, itVertices, itVerticesEnd );
53//-----------------------------------------------------------------------------
54template <typename TRealPoint, typename TRealVector>
55template <typename RealPointIterator, typename VerticesIterator>
57DGtal::SurfaceMesh<TRealPoint, TRealVector>::
58init( RealPointIterator itPos, RealPointIterator itPosEnd,
59 VerticesIterator itVertices, VerticesIterator itVerticesEnd )
62 myPositions = std::vector< RealPoint >( itPos, itPosEnd );
63 myIncidentFaces.resize( myPositions.size() );
64 Index f = 0; // current face index
66 for ( ; itVertices != itVerticesEnd; ++itVertices, ++f )
69 for ( auto it = itVertices->begin(), itE = itVertices->end(); it != itE; ++it )
72 if ( vtx >= nbVertices() )
74 trace.warning() << "[SurfaceMesh::init] Invalid vtx "
75 << vtx << " at face " << f
76 << " since #V=" << nbVertices()
77 << ". Ignoring vertex." << std::endl;
82 myIncidentFaces[ vtx ].push_back( f );
83 f_vtcs.push_back( vtx );
86 myIncidentVertices.push_back( f_vtcs );
93//-----------------------------------------------------------------------------
94template <typename TRealPoint, typename TRealVector>
96DGtal::SurfaceMesh<TRealPoint, TRealVector>::
99 myIncidentVertices.clear();
100 myIncidentFaces.clear();
102 myVertexNormals.clear();
103 myFaceNormals.clear();
104 myNeighborFaces.clear();
105 myNeighborVertices.clear();
106 myEdgeVertices.clear();
108 myEdgeRightFaces.clear();
109 myEdgeLeftFaces.clear();
112//-----------------------------------------------------------------------------
113template <typename TRealPoint, typename TRealVector>
114template <typename RealVectorIterator>
116DGtal::SurfaceMesh<TRealPoint, TRealVector>::
117setVertexNormals( RealVectorIterator itN, RealVectorIterator itNEnd )
119 myVertexNormals = std::vector< RealVector >( itN, itNEnd );
120 return myVertexNormals.size() == myPositions.size();
123//-----------------------------------------------------------------------------
124template <typename TRealPoint, typename TRealVector>
125template <typename RealVectorIterator>
127DGtal::SurfaceMesh<TRealPoint, TRealVector>::
128setFaceNormals( RealVectorIterator itN, RealVectorIterator itNEnd )
130 myFaceNormals = std::vector< RealVector >( itN, itNEnd );
131 return myFaceNormals.size() == myIncidentVertices.size();
134//-----------------------------------------------------------------------------
135template <typename TRealPoint, typename TRealVector>
137DGtal::SurfaceMesh<TRealPoint, TRealVector>::
138computeFaceNormalsFromPositions()
140 myFaceNormals.resize( myIncidentVertices.size() );
142 for ( auto face : myIncidentVertices )
144 RealPoint p; // barycenter
145 RealVector n; // normal
146 // compute barycenter
147 for ( auto idx : face ) p += myPositions[ idx ];
149 // compute normal as sum of triangle normal vectors.
150 for ( Index i = 0; i < face.size(); ++i )
152 const Index j = face[ i ];
153 const Index nj = face[ (i+1) % face.size() ];
154 n += ( myPositions[ j ] - p ).crossProduct( myPositions[ nj ] - p );
156 auto n_norm = n.norm();
157 myFaceNormals[ f ] = n_norm != 0.0 ? n / n_norm : n;
162//-----------------------------------------------------------------------------
163template <typename TRealPoint, typename TRealVector>
165DGtal::SurfaceMesh<TRealPoint, TRealVector>::
166computeFaceNormalsFromVertexNormals()
168 if ( myVertexNormals.empty() ) return;
169 myFaceNormals.resize( myIncidentVertices.size() );
171 for ( auto face : myIncidentVertices )
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;
180//-----------------------------------------------------------------------------
181template <typename TRealPoint, typename TRealVector>
183DGtal::SurfaceMesh<TRealPoint, TRealVector>::
184computeVertexNormalsFromFaceNormals()
186 if ( myFaceNormals.empty() ) return;
187 myVertexNormals.resize( myIncidentFaces.size() );
189 for ( auto vertex : myIncidentFaces )
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;
198//-----------------------------------------------------------------------------
199template <typename TRealPoint, typename TRealVector>
201DGtal::SurfaceMesh<TRealPoint, TRealVector>::
202computeVertexNormalsFromFaceNormalsWithMaxWeights()
204 if ( myFaceNormals.empty() ) return;
205 myVertexNormals.resize( myIncidentFaces.size() );
207 for ( auto incident_faces : myIncidentFaces )
209 RealVector n; // normal
210 const auto weights = getMaxWeights( v );
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;
219//-----------------------------------------------------------------------------
220template <typename TRealPoint, typename TRealVector>
221typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalars
222DGtal::SurfaceMesh<TRealPoint, TRealVector>::
223getMaxWeights( Index v ) const
226 const auto & neighbors = myNeighborVertices[ v ];
227 const RealPoint x = myPositions[ v ];
228 for ( auto idx_f : myIncidentFaces[ v ] )
230 // Find adjacent vertices to v
231 std::vector< Index > adj_vertices;
232 for ( auto idx_v : myIncidentVertices[ idx_f ] )
234 auto it = std::find( neighbors.cbegin(), neighbors.cend(), idx_v );
235 if ( it != neighbors.cend() ) adj_vertices.push_back( *it );
237 if ( adj_vertices.size() != 2 )
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;
245 if (adj_vertices.size() >= 2 )
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 );
253 else weights.push_back( 0.0 );
258//-----------------------------------------------------------------------------
259template <typename TRealPoint, typename TRealVector>
260template <typename AnyRing>
262DGtal::SurfaceMesh<TRealPoint, TRealVector>::
263computeFaceValuesFromVertexValues( const std::vector<AnyRing>& vvalues ) const
265 ASSERT( vvalues.size() == nbVertices() );
266 std::vector<AnyRing> fvalues( nbFaces() );
268 for ( auto face : myIncidentVertices )
270 AnyRing n = NumberTraits<AnyRing>::ZERO;
271 for ( auto idx : face ) n += vvalues[ idx ];
272 fvalues[ f++ ] = n / face.size();
277template <typename TRealPoint, typename TRealVector>
278template <typename AnyRing>
280DGtal::SurfaceMesh<TRealPoint, TRealVector>::
281computeVertexValuesFromFaceValues( const std::vector<AnyRing>& fvalues ) const
283 ASSERT( fvalues.size() == nbFaces() );
284 std::vector<AnyRing> vvalues( nbVertices() );
286 for ( auto vertex : myIncidentFaces )
288 AnyRing n = NumberTraits<AnyRing>::ZERO;
289 for ( auto idx : vertex ) n += fvalues[ idx ];
290 vvalues[ v++ ] = n / vertex.size();
295//-----------------------------------------------------------------------------
296template <typename TRealPoint, typename TRealVector>
297std::vector<TRealVector>
298DGtal::SurfaceMesh<TRealPoint, TRealVector>::
299computeFaceUnitVectorsFromVertexUnitVectors
300( const std::vector<RealVector>& vuvectors ) const
302 ASSERT( vuvectors.size() == nbVertices() );
303 std::vector<RealVector> fuvectors( nbFaces() );
305 for ( auto face : myIncidentVertices )
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;
315//-----------------------------------------------------------------------------
316template <typename TRealPoint, typename TRealVector>
317std::vector<TRealVector>
318DGtal::SurfaceMesh<TRealPoint, TRealVector>::
319computeVertexUnitVectorsFromFaceUnitVectors
320( const std::vector<RealVector>& fuvectors ) const
322 ASSERT( fuvectors.size() == nbFaces() );
323 std::vector<RealVector> vuvectors( nbVertices() );
325 for ( auto vertex : myIncidentFaces )
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;
336//-----------------------------------------------------------------------------
337template <typename TRealPoint, typename TRealVector>
338typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edge
339DGtal::SurfaceMesh<TRealPoint, TRealVector>::
340makeEdge( Vertex i, Vertex j ) const
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();
348//-----------------------------------------------------------------------------
349template <typename TRealPoint, typename TRealVector>
350typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
351DGtal::SurfaceMesh<TRealPoint, TRealVector>::
352averageEdgeLength() const
354 double lengths = 0.0;
355 for ( Edge e = 0; e < nbEdges(); ++e )
357 auto vtcs = edgeVertices( e );
358 const RealPoint p = myPositions[ vtcs.first ];
359 const RealVector pq = myPositions[ vtcs.second ] - p;
360 lengths += pq.norm();
362 lengths /= nbEdges();
366//-----------------------------------------------------------------------------
367template <typename TRealPoint, typename TRealVector>
368typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
369DGtal::SurfaceMesh<TRealPoint, TRealVector>::
370localWindow( Face f ) const
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();
379//-----------------------------------------------------------------------------
380template <typename TRealPoint, typename TRealVector>
382DGtal::SurfaceMesh<TRealPoint, TRealVector>::
383perturbateWithUniformRandomNoise( Scalar p )
385 for ( auto& x : myPositions )
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;
394//-----------------------------------------------------------------------------
395template <typename TRealPoint, typename TRealVector>
397DGtal::SurfaceMesh<TRealPoint, TRealVector>::
398perturbateWithAdaptiveUniformRandomNoise( Scalar p )
400 Scalars local_average_lengths( nbVertices() );
401 for ( Index v = 0; v < nbVertices(); ++v )
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();
409 for ( Index v = 0; v < nbVertices(); ++v )
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;
418//-----------------------------------------------------------------------------
419template <typename TRealPoint, typename TRealVector>
420typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::RealPoint
421DGtal::SurfaceMesh<TRealPoint, TRealVector>::
422faceCentroid( Index f ) const
425 for ( auto v : myIncidentVertices[ f ] )
426 c += myPositions[ v ];
427 return c / myIncidentVertices[ f ].size();
430//-----------------------------------------------------------------------------
431template <typename TRealPoint, typename TRealVector>
432typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::RealPoint
433DGtal::SurfaceMesh<TRealPoint, TRealVector>::
434edgeCentroid( Index e ) const
436 const auto& vtcs = myEdgeVertices[ e ];
437 RealPoint c = myPositions[ vtcs.first ] + myPositions[ vtcs.second ];
441//-----------------------------------------------------------------------------
442template <typename TRealPoint, typename TRealVector>
443typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
444DGtal::SurfaceMesh<TRealPoint, TRealVector>::
445faceArea( Index f ) const
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();
457//-----------------------------------------------------------------------------
458template <typename TRealPoint, typename TRealVector>
459typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces
460DGtal::SurfaceMesh<TRealPoint, TRealVector>::
461computeFacesInclusionsInBall( Scalar r, Index f ) const
463 return computeFacesInclusionsInBall( r, f, faceCentroid( f ) );
466//-----------------------------------------------------------------------------
467template <typename TRealPoint, typename TRealVector>
468typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces
469DGtal::SurfaceMesh<TRealPoint, TRealVector>::
470computeFacesInclusionsInBall( Scalar r, Index f, RealPoint p ) const
472 WeightedFaces result;
475 result.push_back( std::make_pair( f, 0.000001 ) );
478 std::unordered_set< Index > marked;
479 std::queue< Index > active;
482 while ( ! active.empty() )
484 Index current = active.front();
486 Scalar weight = faceInclusionRatio( p, r, current );
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() )
502//-----------------------------------------------------------------------------
503template <typename TRealPoint, typename TRealVector>
505< typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Vertices,
506 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedEdges,
507 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces >
508DGtal::SurfaceMesh<TRealPoint, TRealVector>::
509computeCellsInclusionsInBall( Scalar r, Index f ) const
511 return computeCellsInclusionsInBall( r, f, faceCentroid( f ) );
514//-----------------------------------------------------------------------------
515template <typename TRealPoint, typename TRealVector>
517< typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Vertices,
518 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedEdges,
519 typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces >
520DGtal::SurfaceMesh<TRealPoint, TRealVector>::
521computeCellsInclusionsInBall( Scalar r, Index f, RealPoint p ) const
524 std::set<Vertex> set_v;
525 WeightedEdges result_e;
526 WeightedFaces result_f;
529 result_f.push_back( std::make_pair( f, 0.000001 ) );
530 return std::make_tuple( result_v, result_e, result_f );
532 std::set< Index > marked;
533 std::queue< Index > active;
536 while ( ! active.empty() )
538 Index current = active.front();
540 Scalar fweight = faceInclusionRatio( p, r, current );
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() )
552 // Taking care of edges and vertices
553 const auto& inc_v = myIncidentVertices[ current ];
554 for ( Size i = 0; i < inc_v.size(); ++i )
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 )
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;
566 Scalar eweight = edgeInclusionRatio( p, r, e_ij );
568 result_e.push_back( std::make_pair( e_ij, eweight ) );
572 result_v = Vertices( set_v.cbegin(), set_v.cend() );
573 return std::make_tuple( result_v, result_e, result_f );
576//-----------------------------------------------------------------------------
577template <typename TRealPoint, typename TRealVector>
578typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
579DGtal::SurfaceMesh<TRealPoint, TRealVector>::
580faceInclusionRatio( RealPoint p, Scalar r, Index f ) const
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 )
588 Scalar d = ( myPositions[ v ] - p ).norm();
589 d_max = std::max( d_max, d );
590 d_min = std::min( d_min, d );
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 );
597//-----------------------------------------------------------------------------
598template <typename TRealPoint, typename TRealVector>
599typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
600DGtal::SurfaceMesh<TRealPoint, TRealVector>::
601edgeInclusionRatio( RealPoint p, Scalar r, Index e ) const
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 );
616//-----------------------------------------------------------------------------
617template <typename TRealPoint, typename TRealVector>
618typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
619DGtal::SurfaceMesh<TRealPoint, TRealVector>::
620vertexInclusionRatio( RealPoint p, Scalar r, Index v ) const
622 const RealPoint b = myPositions[ v ];
623 return ( ( b - p ).norm() <= r ) ? 1.0 : 0.0;
627///////////////////////////////////////////////////////////////////////////////
628// Interface - public :
630//-----------------------------------------------------------------------------
631template <typename TRealPoint, typename TRealVector>
633DGtal::SurfaceMesh<TRealPoint, TRealVector>::
634selfDisplay ( std::ostream & out ) const
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();
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();
649 nb_nv /= nbVertices();
651 out << " E[IF]=" << nb_nf << " E[IV]=" << nb_nv << " E[IFE]=" << nb_nfe;
655//-----------------------------------------------------------------------------
656template <typename TRealPoint, typename TRealVector>
658DGtal::SurfaceMesh<TRealPoint, TRealVector>::
661 return myPositions.size() == myIncidentFaces.size()
662 && ( myVertexNormals.size() == 0
663 || ( myVertexNormals.size() == myPositions.size() ) )
664 && ( myFaceNormals.size() == 0
665 || ( myFaceNormals.size() == myIncidentVertices.size() ) );
669//-----------------------------------------------------------------------------
670template <typename TRealPoint, typename TRealVector>
672DGtal::SurfaceMesh<TRealPoint, TRealVector>::
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() )
681 const Size nb_iv = incident_vertices.size();
682 for ( Size k = 0; k < nb_iv; ++k )
684 tmp[ incident_vertices[ k ] ].insert( incident_vertices[ (k+1)%nb_iv ] );
685 tmp[ incident_vertices[ (k+1)%nb_iv ] ].insert( incident_vertices[ k ] );
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() );
692 // For each face, computes its neighboring faces
694 for ( auto incident_vertices : allIncidentVertices() )
696 std::set< Index > neighbor_faces_set;
697 std::sort( incident_vertices.begin(), incident_vertices.end() );
698 for ( auto idx_v : incident_vertices )
700 const auto & incident_faces = incidentFaces( idx_v );
701 for ( auto inc_f : incident_faces )
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() );
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 );
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;
720//-----------------------------------------------------------------------------
721template <typename TRealPoint, typename TRealVector>
723DGtal::SurfaceMesh<TRealPoint, TRealVector>::
726 std::map< VertexPair, std::vector<Face> > edge2face_right;
727 std::map< VertexPair, std::vector<Face> > edge2face_left;
728 std::set< VertexPair > edges;
730 for ( auto incident_vertices : allIncidentVertices() )
732 const Size n = incident_vertices.size();
733 for ( Size i = 0; i < n; i++ )
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 );
740 std::swap( e.first, e.second );
741 edge2face_right[ e ].push_back( idx_f );
747 const Size nbe = edges.size();
748 myEdgeVertices.resize ( nbe );
749 myEdgeFaces.resize ( nbe );
750 myEdgeRightFaces.resize( nbe );
751 myEdgeLeftFaces.resize ( nbe );
753 for ( auto e : edges )
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() );
766//-----------------------------------------------------------------------------
767template <typename TRealPoint, typename TRealVector>
768typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
769DGtal::SurfaceMesh<TRealPoint, TRealVector>::
770computeManifoldBoundaryEdges() const
773 for ( Index e = 0; e < myEdgeFaces.size(); ++e )
775 if ( myEdgeFaces[ e ].size() == 1 )
776 edges.push_back( e );
781//-----------------------------------------------------------------------------
782template <typename TRealPoint, typename TRealVector>
783typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
784DGtal::SurfaceMesh<TRealPoint, TRealVector>::
785computeManifoldInnerEdges() const
788 for ( Index e = 0; e < myEdgeFaces.size(); ++e )
790 if ( myEdgeFaces[ e ].size() == 2 )
791 edges.push_back( e );
796//-----------------------------------------------------------------------------
797template <typename TRealPoint, typename TRealVector>
798typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
799DGtal::SurfaceMesh<TRealPoint, TRealVector>::
800computeManifoldInnerConsistentEdges() const
803 for ( Index e = 0; e < myEdgeFaces.size(); ++e )
805 if ( ( myEdgeRightFaces[ e ].size() == 1 )
806 && ( myEdgeLeftFaces[ e ].size() == 1 ) )
807 edges.push_back( e );
812//-----------------------------------------------------------------------------
813template <typename TRealPoint, typename TRealVector>
814typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
815DGtal::SurfaceMesh<TRealPoint, TRealVector>::
816computeManifoldInnerUnconsistentEdges() const
819 for ( Index e = 0; e < myEdgeFaces.size(); ++e )
821 if ( ( myEdgeRightFaces[ e ].size() == 2 )
822 || ( myEdgeLeftFaces[ e ].size() == 2 ) )
823 edges.push_back( e );
828//-----------------------------------------------------------------------------
829template <typename TRealPoint, typename TRealVector>
830typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
831DGtal::SurfaceMesh<TRealPoint, TRealVector>::
832computeNonManifoldEdges() const
835 for ( Index e = 0; e < myEdgeFaces.size(); ++e )
837 if ( myEdgeFaces[ e ].size() >= 3 )
838 edges.push_back( e );
843///////////////////////////////////////////////////////////////////////////////
844// Implementation of inline functions //
846//-----------------------------------------------------------------------------
847template <typename TRealPoint, typename TRealVector>
849DGtal::operator<< ( std::ostream & out,
850 const SurfaceMesh<TRealPoint, TRealVector> & object )
852 object.selfDisplay( out );
857///////////////////////////////////////////////////////////////////////////////
858///////////////////////////////////////////////////////////////////////////////