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