DGtal 1.3.0
Loading...
Searching...
No Matches
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//-----------------------------------------------------------------------------
43template <typename TRealPoint, typename TRealVector>
44template <typename RealPointIterator, typename VerticesIterator>
45DGtal::SurfaceMesh<TRealPoint, TRealVector>::
46SurfaceMesh( 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//-----------------------------------------------------------------------------
54template <typename TRealPoint, typename TRealVector>
55template <typename RealPointIterator, typename VerticesIterator>
56bool
57DGtal::SurfaceMesh<TRealPoint, TRealVector>::
58init( 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//-----------------------------------------------------------------------------
94template <typename TRealPoint, typename TRealVector>
95void
96DGtal::SurfaceMesh<TRealPoint, TRealVector>::
97clear()
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//-----------------------------------------------------------------------------
113template <typename TRealPoint, typename TRealVector>
114template <typename RealVectorIterator>
115bool
116DGtal::SurfaceMesh<TRealPoint, TRealVector>::
117setVertexNormals( RealVectorIterator itN, RealVectorIterator itNEnd )
118{
119 myVertexNormals = std::vector< RealVector >( itN, itNEnd );
120 return myVertexNormals.size() == myPositions.size();
121}
122
123//-----------------------------------------------------------------------------
124template <typename TRealPoint, typename TRealVector>
125template <typename RealVectorIterator>
126bool
127DGtal::SurfaceMesh<TRealPoint, TRealVector>::
128setFaceNormals( RealVectorIterator itN, RealVectorIterator itNEnd )
129{
130 myFaceNormals = std::vector< RealVector >( itN, itNEnd );
131 return myFaceNormals.size() == myIncidentVertices.size();
132}
133
134//-----------------------------------------------------------------------------
135template <typename TRealPoint, typename TRealVector>
136void
137DGtal::SurfaceMesh<TRealPoint, TRealVector>::
138computeFaceNormalsFromPositions()
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//-----------------------------------------------------------------------------
163template <typename TRealPoint, typename TRealVector>
164void
165DGtal::SurfaceMesh<TRealPoint, TRealVector>::
166computeFaceNormalsFromVertexNormals()
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//-----------------------------------------------------------------------------
181template <typename TRealPoint, typename TRealVector>
182void
183DGtal::SurfaceMesh<TRealPoint, TRealVector>::
184computeVertexNormalsFromFaceNormals()
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//-----------------------------------------------------------------------------
199template <typename TRealPoint, typename TRealVector>
200void
201DGtal::SurfaceMesh<TRealPoint, TRealVector>::
202computeVertexNormalsFromFaceNormalsWithMaxWeights()
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//-----------------------------------------------------------------------------
220template <typename TRealPoint, typename TRealVector>
221typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalars
222DGtal::SurfaceMesh<TRealPoint, TRealVector>::
223getMaxWeights( 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//-----------------------------------------------------------------------------
259template <typename TRealPoint, typename TRealVector>
260template <typename AnyRing>
261std::vector<AnyRing>
262DGtal::SurfaceMesh<TRealPoint, TRealVector>::
263computeFaceValuesFromVertexValues( 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
277template <typename TRealPoint, typename TRealVector>
278template <typename AnyRing>
279std::vector<AnyRing>
280DGtal::SurfaceMesh<TRealPoint, TRealVector>::
281computeVertexValuesFromFaceValues( 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//-----------------------------------------------------------------------------
296template <typename TRealPoint, typename TRealVector>
297std::vector<TRealVector>
298DGtal::SurfaceMesh<TRealPoint, TRealVector>::
299computeFaceUnitVectorsFromVertexUnitVectors
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//-----------------------------------------------------------------------------
316template <typename TRealPoint, typename TRealVector>
317std::vector<TRealVector>
318DGtal::SurfaceMesh<TRealPoint, TRealVector>::
319computeVertexUnitVectorsFromFaceUnitVectors
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//-----------------------------------------------------------------------------
337template <typename TRealPoint, typename TRealVector>
338typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edge
339DGtal::SurfaceMesh<TRealPoint, TRealVector>::
340makeEdge( 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//-----------------------------------------------------------------------------
349template <typename TRealPoint, typename TRealVector>
350typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
351DGtal::SurfaceMesh<TRealPoint, TRealVector>::
352averageEdgeLength() 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//-----------------------------------------------------------------------------
367template <typename TRealPoint, typename TRealVector>
368typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
369DGtal::SurfaceMesh<TRealPoint, TRealVector>::
370localWindow( 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//-----------------------------------------------------------------------------
380template <typename TRealPoint, typename TRealVector>
381void
382DGtal::SurfaceMesh<TRealPoint, TRealVector>::
383perturbateWithUniformRandomNoise( 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//-----------------------------------------------------------------------------
395template <typename TRealPoint, typename TRealVector>
396void
397DGtal::SurfaceMesh<TRealPoint, TRealVector>::
398perturbateWithAdaptiveUniformRandomNoise( 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//-----------------------------------------------------------------------------
419template <typename TRealPoint, typename TRealVector>
420typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::RealPoint
421DGtal::SurfaceMesh<TRealPoint, TRealVector>::
422faceCentroid( 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//-----------------------------------------------------------------------------
431template <typename TRealPoint, typename TRealVector>
432typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::RealPoint
433DGtal::SurfaceMesh<TRealPoint, TRealVector>::
434edgeCentroid( 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//-----------------------------------------------------------------------------
442template <typename TRealPoint, typename TRealVector>
443typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
444DGtal::SurfaceMesh<TRealPoint, TRealVector>::
445faceArea( 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//-----------------------------------------------------------------------------
458template <typename TRealPoint, typename TRealVector>
459typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces
460DGtal::SurfaceMesh<TRealPoint, TRealVector>::
461computeFacesInclusionsInBall( Scalar r, Index f ) const
462{
463 return computeFacesInclusionsInBall( r, f, faceCentroid( f ) );
464}
465
466//-----------------------------------------------------------------------------
467template <typename TRealPoint, typename TRealVector>
468typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::WeightedFaces
469DGtal::SurfaceMesh<TRealPoint, TRealVector>::
470computeFacesInclusionsInBall( 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//-----------------------------------------------------------------------------
503template <typename TRealPoint, typename TRealVector>
504std::tuple
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
510{
511 return computeCellsInclusionsInBall( r, f, faceCentroid( f ) );
512}
513
514//-----------------------------------------------------------------------------
515template <typename TRealPoint, typename TRealVector>
516std::tuple
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
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//-----------------------------------------------------------------------------
577template <typename TRealPoint, typename TRealVector>
578typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
579DGtal::SurfaceMesh<TRealPoint, TRealVector>::
580faceInclusionRatio( 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//-----------------------------------------------------------------------------
598template <typename TRealPoint, typename TRealVector>
599typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
600DGtal::SurfaceMesh<TRealPoint, TRealVector>::
601edgeInclusionRatio( 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//-----------------------------------------------------------------------------
617template <typename TRealPoint, typename TRealVector>
618typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Scalar
619DGtal::SurfaceMesh<TRealPoint, TRealVector>::
620vertexInclusionRatio( 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//-----------------------------------------------------------------------------
631template <typename TRealPoint, typename TRealVector>
632void
633DGtal::SurfaceMesh<TRealPoint, TRealVector>::
634selfDisplay ( 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//-----------------------------------------------------------------------------
656template <typename TRealPoint, typename TRealVector>
657bool
658DGtal::SurfaceMesh<TRealPoint, TRealVector>::
659isValid() 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//-----------------------------------------------------------------------------
670template <typename TRealPoint, typename TRealVector>
671void
672DGtal::SurfaceMesh<TRealPoint, TRealVector>::
673computeNeighbors()
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//-----------------------------------------------------------------------------
721template <typename TRealPoint, typename TRealVector>
722void
723DGtal::SurfaceMesh<TRealPoint, TRealVector>::
724computeEdges()
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//-----------------------------------------------------------------------------
767template <typename TRealPoint, typename TRealVector>
768typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
769DGtal::SurfaceMesh<TRealPoint, TRealVector>::
770computeManifoldBoundaryEdges() 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//-----------------------------------------------------------------------------
782template <typename TRealPoint, typename TRealVector>
783typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
784DGtal::SurfaceMesh<TRealPoint, TRealVector>::
785computeManifoldInnerEdges() 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//-----------------------------------------------------------------------------
797template <typename TRealPoint, typename TRealVector>
798typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
799DGtal::SurfaceMesh<TRealPoint, TRealVector>::
800computeManifoldInnerConsistentEdges() 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//-----------------------------------------------------------------------------
813template <typename TRealPoint, typename TRealVector>
814typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
815DGtal::SurfaceMesh<TRealPoint, TRealVector>::
816computeManifoldInnerUnconsistentEdges() 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//-----------------------------------------------------------------------------
829template <typename TRealPoint, typename TRealVector>
830typename DGtal::SurfaceMesh<TRealPoint, TRealVector>::Edges
831DGtal::SurfaceMesh<TRealPoint, TRealVector>::
832computeNonManifoldEdges() 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//-----------------------------------------------------------------------------
847template <typename TRealPoint, typename TRealVector>
848std::ostream&
849DGtal::operator<< ( std::ostream & out,
850 const SurfaceMesh<TRealPoint, TRealVector> & object )
851{
852 object.selfDisplay( out );
853 return out;
854}
855
856
857///////////////////////////////////////////////////////////////////////////////
858///////////////////////////////////////////////////////////////////////////////