DGtal 1.3.0
Loading...
Searching...
No Matches
Surfaces.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 Surfaces.ih
19 * @author Jacques-Olivier Lachaud (\c jacques-olivier.lachaud@univ-savoie.fr )
20 * Laboratory of Mathematics (CNRS, UMR 5807), University of Savoie, France
21 *
22 * @date 2011/03/19
23 *
24 * Implementation of inline methods defined in Surfaces.h
25 *
26 * This file is part of the DGtal library.
27 */
28
29
30//////////////////////////////////////////////////////////////////////////////
31#include <cstdlib>
32#include <vector>
33#include <queue>
34#include <algorithm>
35#include "DGtal/kernel/CPointPredicate.h"
36#include "DGtal/images/imagesSetsUtils/ImageFromSet.h"
37#include "DGtal/images/ImageSelector.h"
38#include "DGtal/topology/CSurfelPredicate.h"
39#include "DGtal/helpers/StdDefs.h"
40
41
42//////////////////////////////////////////////////////////////////////////////
43
44///////////////////////////////////////////////////////////////////////////////
45// IMPLEMENTATION of inline methods.
46///////////////////////////////////////////////////////////////////////////////
47
48///////////////////////////////////////////////////////////////////////////////
49// ----------------------- Standard services ------------------------------
50
51/**
52 * Destructor.
53 */
54template <typename TKSpace>
55inline
56DGtal::Surfaces<TKSpace>::~Surfaces()
57{
58}
59//-----------------------------------------------------------------------------
60template <typename TKSpace>
61template <typename PointPredicate>
62typename DGtal::Surfaces<TKSpace>::SCell
63DGtal::Surfaces<TKSpace>::
64findABel
65( const KSpace & K,
66 const PointPredicate & pp,
67 unsigned int nbtries)
68{
69 BOOST_CONCEPT_ASSERT(( concepts::CPointPredicate<PointPredicate> ));
70
71 DGtal::InputException dgtalerror;
72 Point sizes = K.upperBound() - K.lowerBound();
73 Point x1 = K.lowerBound();
74 Point x2;
75 // (1) Find two candidates in the space.
76 bool val_v1 = pp( x1 ); // dset.find( x1 ) != dset.end();
77 bool found = false;
78 Integer r;
79 for ( unsigned int j = 0; j < nbtries; ++j )
80 {
81 for ( Dimension i = 0; i < K.dimension; ++i )
82 {
83 r = rand();
84 x2[ i ] = ( r % sizes[ i ] ) + K.min( i );
85 }
86 bool val_v2 = pp( x2 ); // dset.find( x2 ) != dset.end();
87 if ( val_v2 != val_v1 )
88 {
89 found = true;
90 break;
91 }
92 }
93 if ( ! found ) throw dgtalerror;
94 return findABel( K, pp, x1, x2 );
95}
96//-----------------------------------------------------------------------------
97template <typename TKSpace>
98template <typename PointPredicate>
99typename DGtal::Surfaces<TKSpace>::SCell
100DGtal::Surfaces<TKSpace>::
101findABel
102( const KSpace & K,
103 const PointPredicate & pp,
104 Point x1, Point x2 )
105{
106 BOOST_CONCEPT_ASSERT(( concepts::CPointPredicate<PointPredicate> ));
107 // (1) Checks the two candidates in the space.
108 bool val_v1 = pp( x1 ); // dset.find( x1 ) != dset.end();
109 ASSERT( val_v1 != pp( x2 ) );
110 // (2) Find two candidates on the same axis.
111 Dimension d = 0;
112 bool alreadyOnSameAxis = true;
113 for ( Dimension i = 0; i < K.dimension; ++i )
114 {
115 if ( x1[ i ] != x2[ i ] )
116 {
117 for ( Dimension j = i + 1; j < K.dimension; ++j )
118 {
119 if ( x1[ j ] != x2[ j ] )
120 {
121 alreadyOnSameAxis = false;
122 Integer c = x2[ j ];
123 x2[ j ] = x1[ j ];
124 bool val_v2 = pp( x2 ); // dset.find( x2 ) != dset.end();
125 if ( val_v2 != val_v1 )
126 { // v2 is updated.
127 d = i;
128 }
129 else
130 { // v1 is updated.
131 x1 = x2;
132 x2[ j ] = c;
133 d = j;
134 }
135 } // if ( x1[ j ] != x2[ j ] )
136 } // for ( Dimension j = i + 1; j < K.dimension; ++j )
137 if ( alreadyOnSameAxis )
138 d = i;
139 } // if ( x1[ i ] != x2[ i ] )
140 } // for ( Dimension i = 0; i < K.dimension; ++i )
141
142 // (3) Check result.
143 for ( Dimension i = 0; i < K.dimension; ++i )
144 {
145 ASSERT( ! ( ( i == d ) && ( x1[ i ] == x2[ i ] ) ) && "[DGtal::Surfaces::findABel] Same position on the main axis." );
146 ASSERT( ! ( ( i != d ) && ( x1[ i ] != x2[ i ] ) ) && "[DGtal::Surfaces::findABel] Different position on other axis." );
147 }
148
149 // (4) Dichotomy
150 Point xmid = x1;
151 while ( true )
152 {
153 xmid[ d ] = ( x1[ d ] + x2[ d ] ) / 2;
154 if ( ( xmid[ d ] == x1[ d ] ) || ( xmid[ d ] == x2[ d ] ) )
155 break;
156 bool val_mid = pp( xmid ); // dset.find( xmid ) != dset.end();
157 if ( val_mid != val_v1 )
158 x2[ d ] = xmid[ d ];
159 else
160 x1[ d ] = xmid[ d ];
161 }
162
163 // (5) Check result.
164 for ( Dimension i = 0; i < K.dimension; ++i )
165 {
166 ASSERT( ! ( ( i == d ) && ( x1[ i ] != x2[ i ] - 1 ) && ( x1[ i ] != x2[ i ] + 1 ) )
167 && "[DGtal::Surfaces::findABel] Points should be adjacent on main axis." );
168 ASSERT( ! ( ( i != d ) && ( x1[ i ] != x2[ i ] ) )
169 && "[DGtal::Surfaces::findABel] Points should have the same coordinates on other axes." );
170 }
171
172 // (6) Computes bel.
173 SCell spel1 = K.sSpel( x1, val_v1 ? K.POS : K.NEG );
174 return K.sIncident( spel1, d, ( x1[ d ] == x2[ d ] - 1 ) );
175}
176//-----------------------------------------------------------------------------
177template <typename TKSpace>
178template <typename SCellSet, typename PointPredicate >
179void
180DGtal::Surfaces<TKSpace>::
181trackBoundary( SCellSet & surface,
182 const KSpace & K,
183 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
184 const PointPredicate & pp,
185 const SCell & start_surfel )
186{
187 BOOST_CONCEPT_ASSERT(( concepts::CPointPredicate<PointPredicate> ));
188
189 SCell b; // current surfel
190 SCell bn; // neighboring surfel
191 ASSERT( K.sIsSurfel( start_surfel ) );
192 surface.clear(); // boundary being extracted.
193
194 SurfelNeighborhood<KSpace> SN;
195 SN.init( &K, &surfel_adj, start_surfel );
196 std::queue<SCell> qbels;
197 qbels.push( start_surfel );
198 surface.insert( start_surfel );
199 // For all pending bels
200 while ( ! qbels.empty() )
201 {
202 b = qbels.front();
203 qbels.pop();
204 SN.setSurfel( b );
205 for ( DirIterator q = K.sDirs( b ); q != 0; ++q )
206 {
207 Dimension track_dir = *q;
208 // ----- 1st pass with positive orientation ------
209 if ( SN.getAdjacentOnPointPredicate( bn, pp, track_dir, true ) )
210 {
211 if ( surface.find( bn ) == surface.end() )
212 {
213 surface.insert( bn );
214 qbels.push( bn );
215 }
216 }
217 // ----- 2nd pass with negative orientation ------
218 if ( SN.getAdjacentOnPointPredicate( bn, pp, track_dir, false ) )
219 {
220 if ( surface.find( bn ) == surface.end() )
221 {
222 surface.insert( bn );
223 qbels.push( bn );
224 }
225 }
226 } // for ( DirIterator q = K.sDirs( b ); q != 0; ++q )
227 } // while ( ! qbels.empty() )
228}
229//-----------------------------------------------------------------------------
230template <typename TKSpace>
231template <typename SCellSet, typename SurfelPredicate >
232void
233DGtal::Surfaces<TKSpace>::
234trackSurface( SCellSet & surface,
235 const KSpace & K,
236 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
237 const SurfelPredicate & sp,
238 const SCell & start_surfel )
239{
240 BOOST_CONCEPT_ASSERT(( concepts::CSurfelPredicate<SurfelPredicate> ));
241
242 SCell b; // current surfel
243 SCell bn; // neighboring surfel
244 ASSERT( K.sIsSurfel( start_surfel ) );
245 surface.clear(); // boundary being extracted.
246
247 SurfelNeighborhood<KSpace> SN;
248 SN.init( &K, &surfel_adj, start_surfel );
249 std::queue<SCell> qbels;
250 qbels.push( start_surfel );
251 surface.insert( start_surfel );
252 // For all pending bels
253 while ( ! qbels.empty() )
254 {
255 b = qbels.front();
256 qbels.pop();
257 SN.setSurfel( b );
258 for ( DirIterator q = K.sDirs( b ); q != 0; ++q )
259 {
260 Dimension track_dir = *q;
261 // ----- 1st pass with positive orientation ------
262 if ( SN.getAdjacentOnSurfelPredicate( bn, sp, track_dir, true ) )
263 {
264 if ( surface.find( bn ) == surface.end() )
265 {
266 surface.insert( bn );
267 qbels.push( bn );
268 }
269 }
270 // ----- 2nd pass with negative orientation ------
271 if ( SN.getAdjacentOnSurfelPredicate( bn, sp, track_dir, false ) )
272 {
273 if ( surface.find( bn ) == surface.end() )
274 {
275 surface.insert( bn );
276 qbels.push( bn );
277 }
278 }
279 } // for ( DirIterator q = K.sDirs( b ); q != 0; ++q )
280 } // while ( ! qbels.empty() )
281}
282
283//-----------------------------------------------------------------------------
284template <typename TKSpace>
285template <typename SCellSet, typename SurfelPredicate >
286void
287DGtal::Surfaces<TKSpace>::
288trackClosedSurface( SCellSet & surface,
289 const KSpace & K,
290 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
291 const SurfelPredicate & sp,
292 const SCell & start_surfel )
293{
294 BOOST_CONCEPT_ASSERT(( concepts::CSurfelPredicate<SurfelPredicate> ));
295
296 SCell b; // current surfel
297 SCell bn; // neighboring surfel
298 ASSERT( K.sIsSurfel( start_surfel ) );
299 surface.clear(); // boundary being extracted.
300
301 SurfelNeighborhood<KSpace> SN;
302 SN.init( &K, &surfel_adj, start_surfel );
303 std::queue<SCell> qbels;
304 qbels.push( start_surfel );
305 surface.insert( start_surfel );
306 // For all pending bels
307 while ( ! qbels.empty() )
308 {
309 b = qbels.front();
310 qbels.pop();
311 SN.setSurfel( b );
312 for ( DirIterator q = K.sDirs( b ); q != 0; ++q )
313 {
314 Dimension track_dir = *q;
315 // ----- direct orientation ------
316 if ( SN.getAdjacentOnSurfelPredicate( bn, sp, track_dir,
317 K.sDirect( b, track_dir ) ) )
318 {
319 if ( surface.find( bn ) == surface.end() )
320 {
321 surface.insert( bn );
322 qbels.push( bn );
323 }
324 }
325 } // for ( DirIterator q = K.sDirs( b ); q != 0; ++q )
326 } // while ( ! qbels.empty() )
327}
328
329
330//-----------------------------------------------------------------------------
331template <typename TKSpace>
332template <typename PointPredicate>
333void
334DGtal::Surfaces<TKSpace>::track2DBoundary
335( std::vector<SCell> & aSCellContour2D,
336 const KSpace & K,
337 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
338 const PointPredicate & pp,
339 const SCell & start_surfel )
340{
341 BOOST_CONCEPT_ASSERT(( concepts::CPointPredicate<PointPredicate> ));
342 ASSERT( K.dimension == 2 );
343
344 SCell b= start_surfel; // current surfel
345 SCell bn; // neighboring surfel
346 ASSERT( K.sIsSurfel( start_surfel ) );
347 // std::set<SCell> setSurface;
348 // setSurface.insert(start_surfel);
349 aSCellContour2D.clear(); // boundary being extracted.
350 aSCellContour2D.push_back(start_surfel);
351 SurfelNeighborhood<KSpace> SN;
352 SN.init( &K, &surfel_adj, start_surfel );
353 SN.setSurfel( b );
354 // search along indirect orientation.
355 bool hasPrevNeighbor = true;
356 while ( hasPrevNeighbor )
357 {
358 hasPrevNeighbor=false;
359 Dimension track_dir = *(K.sDirs( b ));
360 SN.setSurfel( b );
361 if ( SN.getAdjacentOnPointPredicate( bn, pp, track_dir,
362 ! K.sDirect( b, track_dir ) ) )
363 {
364 if ( bn != start_surfel )
365 // if ( setSurface.find( bn ) == setSurface.end() )
366 {
367 hasPrevNeighbor=true;
368 aSCellContour2D.push_back( bn );
369 // setSurface.insert(bn);
370 }
371 }
372 b = bn;
373 }
374 // since the contour is not necessary closed we search in the other direction.
375 reverse(aSCellContour2D.begin(), aSCellContour2D.end());
376 if ( b != start_surfel )
377 { // the contour is necessarily open.
378 b = start_surfel;
379 bool hasNextNeighbor = true;
380 while ( hasNextNeighbor )
381 {
382 hasNextNeighbor=false;
383 Dimension track_dir = *(K.sDirs( b ));
384 SN.setSurfel( b );
385 if ( SN.getAdjacentOnPointPredicate( bn, pp, track_dir,
386 K.sDirect( b, track_dir ) ) )
387 {
388 // if ( setSurface.find( bn ) == setSurface.end() )
389 // {
390 aSCellContour2D.push_back( bn );
391 hasNextNeighbor=true;
392 // setSurface.insert(bn);
393 // }
394 }
395 b=bn;
396 }
397 }
398}
399
400
401
402//-----------------------------------------------------------------------------
403template <typename TKSpace>
404template <typename PointPredicate>
405void
406DGtal::Surfaces<TKSpace>::track2DSliceBoundary
407( std::vector<SCell> & aSCellContour2D,
408 const KSpace & K,
409 const Dimension & trackDir,
410 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
411 const PointPredicate & pp,
412 const SCell & start_surfel )
413{
414 BOOST_CONCEPT_ASSERT(( concepts::CPointPredicate<PointPredicate> ));
415 SCell b= start_surfel; // current surfel
416 SCell bn; // neighboring surfel
417 ASSERT( K.sIsSurfel( start_surfel ) );
418 // std::set<SCell> setSurface;
419 // setSurface.insert(start_surfel);
420 aSCellContour2D.clear(); // boundary being extracted.
421 aSCellContour2D.push_back(start_surfel);
422 SurfelNeighborhood<KSpace> SN;
423 SN.init( &K, &surfel_adj, start_surfel );
424 SN.setSurfel( b );
425 Dimension orthDir = K.sOrthDir( start_surfel );
426 bool hasPrevNeighbor = true;
427 while ( hasPrevNeighbor )
428 {
429 hasPrevNeighbor=false;
430 // search a tracking direction compatible with track/orth direction
431 Dimension track_dir = K.sOrthDir( b ) == orthDir ? trackDir : orthDir;
432 SN.setSurfel( b );
433 if ( SN.getAdjacentOnPointPredicate( bn, pp, track_dir,
434 !K.sDirect( b, track_dir ) ) )
435 {
436 if ( bn != start_surfel )
437 // if ( setSurface.find( bn ) == setSurface.end() )
438 {
439 hasPrevNeighbor=true;
440 aSCellContour2D.push_back( bn );
441 // setSurface.insert(bn);
442 }
443 }
444 b = bn;
445 }
446 // since the contour is not necessary closed we search in the other direction.
447 reverse(aSCellContour2D.begin(), aSCellContour2D.end());
448 if ( b != start_surfel )
449 { // the contour is necessarily open.
450 b = start_surfel;
451 bool hasNextNeighbor = true;
452 while ( hasNextNeighbor )
453 {
454 hasNextNeighbor=false;
455 // search a tracking direction compatible with constant direction
456 Dimension track_dir = K.sOrthDir( b ) == orthDir ? trackDir : orthDir;
457 SN.setSurfel( b );
458 if ( SN.getAdjacentOnPointPredicate( bn, pp, track_dir,
459 K.sDirect( b, track_dir ) ) )
460 {
461 // if ( setSurface.find( bn ) == setSurface.end() )
462 // {
463 aSCellContour2D.push_back( bn );
464 // setSurface.insert(bn);
465 hasNextNeighbor=true;
466 // }
467 }
468 b=bn;
469 }
470 }
471}
472
473//-----------------------------------------------------------------------------
474template <typename TKSpace>
475template <typename SurfelPredicate >
476inline
477void
478DGtal::Surfaces<TKSpace>::
479track2DSurface( std::vector<SCell> & aSCellContour,
480 const KSpace & K,
481 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
482 const SurfelPredicate & sp,
483 const SCell & start_surfel )
484{
485 BOOST_CONCEPT_ASSERT(( concepts::CSurfelPredicate<SurfelPredicate> ));
486 ASSERT( KSpace::dimension == 2 );
487
488 SCell b= start_surfel; // current surfel
489 SCell bn; // neighboring surfel
490 ASSERT( K.sIsSurfel( start_surfel ) );
491 aSCellContour.clear(); // boundary being extracted.
492 aSCellContour.push_back(start_surfel);
493 SurfelNeighborhood<KSpace> SN;
494 SN.init( &K, &surfel_adj, start_surfel );
495 SN.setSurfel( b );
496 // search along indirect orientation.
497 bool hasPrevNeighbor = true;
498 while ( hasPrevNeighbor )
499 {
500 hasPrevNeighbor=false;
501 Dimension track_dir = *(K.sDirs( b ));
502 SN.setSurfel( b );
503 if ( SN.getAdjacentOnSurfelPredicate( bn, sp, track_dir,
504 ! K.sDirect( b, track_dir ) ) )
505 {
506 if ( bn != start_surfel )
507 {
508 hasPrevNeighbor=true;
509 aSCellContour.push_back( bn );
510 }
511 }
512 b = bn;
513 }
514 // since the contour is not necessary closed we search in the other direction.
515 reverse( aSCellContour.begin(), aSCellContour.end() );
516 if ( b != start_surfel )
517 { // the contour is necessarily open.
518 b = start_surfel;
519 bool hasNextNeighbor = true;
520 while ( hasNextNeighbor )
521 {
522 hasNextNeighbor=false;
523 Dimension track_dir = *(K.sDirs( b ));
524 SN.setSurfel( b );
525 if ( SN.getAdjacentOnSurfelPredicate( bn, sp, track_dir,
526 K.sDirect( b, track_dir ) ) )
527 {
528 aSCellContour.push_back( bn );
529 hasNextNeighbor=true;
530 }
531 b=bn;
532 }
533 }
534}
535
536//-----------------------------------------------------------------------------
537template <typename TKSpace>
538template <typename SurfelPredicate >
539inline
540void
541DGtal::Surfaces<TKSpace>::
542track2DSliceSurface( std::vector<SCell> & aSCellContour,
543 const KSpace & K,
544 const Dimension & trackDir,
545 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
546 const SurfelPredicate & sp,
547 const SCell & start_surfel )
548{
549 BOOST_CONCEPT_ASSERT(( concepts::CSurfelPredicate<SurfelPredicate> ));
550 SCell b= start_surfel; // current surfel
551 SCell bn; // neighboring surfel
552 ASSERT( K.sIsSurfel( start_surfel ) );
553 aSCellContour.clear(); // boundary being extracted.
554 aSCellContour.push_back(start_surfel);
555 SurfelNeighborhood<KSpace> SN;
556 SN.init( &K, &surfel_adj, start_surfel );
557 SN.setSurfel( b );
558 Dimension orthDir = K.sOrthDir( start_surfel );
559 bool hasPrevNeighbor = true;
560 while ( hasPrevNeighbor )
561 {
562 hasPrevNeighbor=false;
563 // search a tracking direction compatible with track/orth direction
564 Dimension track_dir = K.sOrthDir( b ) == orthDir ? trackDir : orthDir;
565 SN.setSurfel( b );
566 if ( SN.getAdjacentOnSurfelPredicate( bn, sp, track_dir,
567 !K.sDirect( b, track_dir ) ) )
568 {
569 if ( bn != start_surfel )
570 {
571 hasPrevNeighbor=true;
572 aSCellContour.push_back( bn );
573 }
574 }
575 b = bn;
576 }
577 // since the contour is not necessary closed we search in the other direction.
578 reverse( aSCellContour.begin(), aSCellContour.end() );
579 if ( b != start_surfel )
580 { // the contour is necessarily open.
581 b = start_surfel;
582 bool hasNextNeighbor = true;
583 while ( hasNextNeighbor )
584 {
585 hasNextNeighbor=false;
586 // search a tracking direction compatible with constant direction
587 Dimension track_dir = K.sOrthDir( b ) == orthDir ? trackDir : orthDir;
588 SN.setSurfel( b );
589 if ( SN.getAdjacentOnSurfelPredicate( bn, sp, track_dir,
590 K.sDirect( b, track_dir ) ) )
591 {
592 aSCellContour.push_back( bn );
593 hasNextNeighbor=true;
594 }
595 b=bn;
596 }
597 }
598}
599
600
601
602//-----------------------------------------------------------------------------
603template <typename TKSpace>
604template <typename PointPredicate>
605inline
606void
607DGtal::Surfaces<TKSpace>::track2DBoundaryPoints
608( std::vector<Point> & aVectorOfPoints,
609 const KSpace & K,
610 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
611 const PointPredicate & pp,
612 const SCell & start_surfel )
613{
614 aVectorOfPoints.clear();
615
616 // Getting the consecutive surfels of the 2D boundary
617 std::vector<SCell> vectBdrySCell;
618 Surfaces<KSpace>::track2DBoundary( vectBdrySCell,
619 K, surfel_adj, pp, start_surfel );
620 typedef typename std::vector<SCell>::const_iterator SCellConstIterator;
621 for ( SCellConstIterator it = vectBdrySCell.begin(),
622 it_end = vectBdrySCell.end();
623 it != it_end; ++it )
624 {
625 Dimension track = *( K.sDirs( *it ) );
626 SCell pointel = K.sIndirectIncident( *it, track );
627 aVectorOfPoints.push_back( K.sCoords( pointel ) );
628 }
629}
630
631
632
633
634//-----------------------------------------------------------------------------
635template <typename TKSpace>
636template <typename PointPredicate>
637void
638DGtal::Surfaces<TKSpace>::
639extractAll2DSCellContours( std::vector< std::vector<SCell> > & aVectSCellContour2D,
640 const KSpace & aKSpace,
641 const SurfelAdjacency<KSpace::dimension> & aSurfelAdj,
642 const PointPredicate & pp )
643{
644 std::set<SCell> bdry;
645 sMakeBoundary( bdry, aKSpace, pp,
646 aKSpace.lowerBound(), aKSpace.upperBound() );
647 aVectSCellContour2D.clear();
648 while( ! bdry.empty() )
649 {
650 std::vector<SCell> aContour;
651 SCell aCell = *(bdry.begin());
652 track2DBoundary( aContour, aKSpace, aSurfelAdj, pp, aCell );
653 aVectSCellContour2D.push_back( aContour );
654 // removing cells from boundary;
655 for( unsigned int i = 0; i < aContour.size(); i++ )
656 {
657 SCell sc = aContour.at(i);
658 bdry.erase(sc);
659 }
660 }
661}
662
663
664
665//-----------------------------------------------------------------------------
666template <typename TKSpace>
667template <typename PointPredicate>
668void
669DGtal::Surfaces<TKSpace>::orientSCellExterior(std::vector<SCell> & aVectOfSCell,
670 const KSpace & aKSpace, const PointPredicate & pp ){
671 for( typename std::vector<SCell>::iterator it = aVectOfSCell.begin();
672 it!=aVectOfSCell.end(); it++){
673 SCell incidUpperDim = aKSpace.sDirectIncident(*it, aKSpace.sOrthDir(*it));
674 if( pp( aKSpace.sCoords(incidUpperDim) )){
675 aKSpace.sSetSign (*it, !aKSpace.sDirect(*it, aKSpace.sOrthDir(*it)) );
676 }else{
677 aKSpace.sSetSign (*it, aKSpace.sDirect(*it, !aKSpace.sOrthDir(*it)) );
678 }
679 }
680}
681
682
683
684
685//-----------------------------------------------------------------------------
686template <typename TKSpace>
687template <typename PointPredicate>
688void
689DGtal::Surfaces<TKSpace>::
690extractAllConnectedSCell
691( std::vector< std::vector<SCell> > & aVectConnectedSCell,
692 const KSpace & aKSpace,
693 const SurfelAdjacency<KSpace::dimension> & aSurfelAdj,
694 const PointPredicate & pp,
695 bool forceOrientCellExterior )
696{
697 std::set<SCell> bdry;
698 sMakeBoundary( bdry, aKSpace, pp,
699 aKSpace.lowerBound(),
700 aKSpace.upperBound() );
701 aVectConnectedSCell.clear();
702 while(!bdry.empty()){
703 std::set<SCell> aConnectedSCellSet;
704 SCell aCell = *(bdry.begin());
705 trackBoundary(aConnectedSCellSet, aKSpace, aSurfelAdj, pp, aCell );
706 //transform into vector<SCell>
707 std::vector<SCell> vCS;
708 for(typename std::set<SCell>::iterator it = aConnectedSCellSet.begin(); it!= aConnectedSCellSet.end(); ++it){
709 vCS.push_back(*it);
710 // removing cells from boundary;
711 bdry.erase(*it);
712 }
713 if(forceOrientCellExterior){
714 orientSCellExterior(vCS, aKSpace, pp);
715 }
716 aVectConnectedSCell.push_back(vCS);
717 }
718}
719
720
721
722
723//-----------------------------------------------------------------------------
724template <typename TKSpace>
725template <typename PointPredicate>
726void
727DGtal::Surfaces<TKSpace>::
728extractAllPointContours4C( std::vector< std::vector< Point > > & aVectPointContour2D,
729 const KSpace & aKSpace,
730 const PointPredicate & pp,
731 const SurfelAdjacency<2> & aSAdj)
732{
733 aVectPointContour2D.clear();
734
735 std::vector< std::vector<SCell> > vectContoursBdrySCell;
736 extractAll2DSCellContours( vectContoursBdrySCell,
737 aKSpace, aSAdj, pp );
738
739 for(unsigned int i=0; i< vectContoursBdrySCell.size(); i++){
740 std::vector< Point > aContour;
741 for(unsigned int j=0; j< vectContoursBdrySCell.at(i).size(); j++){
742 SCell sc = vectContoursBdrySCell.at(i).at(j);
743 float x = (float)
744 ( NumberTraits<typename TKSpace::Integer>::castToInt64_t( aKSpace.sKCoord(sc, 0) ) >> 1 );
745 float y = (float)
746 ( NumberTraits<typename TKSpace::Integer>::castToInt64_t( aKSpace.sKCoord(sc, 1) ) >> 1 );
747 bool xodd = ( aKSpace.sKCoord(sc, 0) & 1 );
748 bool yodd = ( aKSpace.sKCoord(sc, 1) & 1 );
749 double x0 = !xodd ? x - 0.5 : (!aKSpace.sSign(sc)? x - 0.5: x + 0.5) ;
750 double y0 = !yodd ? y - 0.5 : (!aKSpace.sSign(sc)? y - 0.5: y + 0.5);
751 double x1 = !xodd ? x - 0.5 : (aKSpace.sSign(sc)? x - 0.5: x + 0.5) ;
752 double y1 = !yodd ? y - 0.5 : (aKSpace.sSign(sc)? y - 0.5: y + 0.5);
753
754 Point ptA((const int)(x0+0.5), (const int)(y0-0.5));
755 Point ptB((const int)(x1+0.5), (const int)(y1-0.5)) ;
756 aContour.push_back(ptA);
757 if(sc== vectContoursBdrySCell.at(i).at(vectContoursBdrySCell.at(i).size()-1)){
758 aContour.push_back(ptB);
759 }
760 }
761 aVectPointContour2D.push_back(aContour);
762 }
763}
764
765
766
767//-----------------------------------------------------------------------------
768template <typename TKSpace>
769template <typename SCellSet, typename PointPredicate >
770void
771DGtal::Surfaces<TKSpace>::
772trackClosedBoundary( SCellSet & surface,
773 const KSpace & K,
774 const SurfelAdjacency<KSpace::dimension> & surfel_adj,
775 const PointPredicate & pp,
776 const SCell & start_surfel )
777{
778 SCell b; // current surfel
779 SCell bn; // neighboring surfel
780 ASSERT( K.sIsSurfel( start_surfel ) );
781 surface.clear(); // boundary being extracted.
782
783 SurfelNeighborhood<KSpace> SN;
784 SN.init( &K, &surfel_adj, start_surfel );
785 std::queue<SCell> qbels;
786 qbels.push( start_surfel );
787 surface.insert( start_surfel );
788 // For all pending bels
789 while ( ! qbels.empty() )
790 {
791 b = qbels.front();
792 qbels.pop();
793 SN.setSurfel( b );
794 for ( DirIterator q = K.sDirs( b ); q != 0; ++q )
795 {
796 Dimension track_dir = *q;
797 // ----- One pass, look for direct orientation ------
798 if ( SN.getAdjacentOnPointPredicate( bn, pp, track_dir,
799 K.sDirect( b, track_dir ) ) )
800 {
801 if ( surface.find( bn ) == surface.end() )
802 {
803 surface.insert( bn );
804 qbels.push( bn );
805 }
806 }
807 } // for ( DirIterator q = K.sDirs( b ); q != 0; ++q )
808 } // while ( ! qbels.empty() )
809}
810
811//-----------------------------------------------------------------------------
812template <typename TKSpace>
813template <typename CellSet, typename PointPredicate >
814void
815DGtal::Surfaces<TKSpace>::
816uMakeBoundary( CellSet & aBoundary,
817 const KSpace & aKSpace,
818 const PointPredicate & pp,
819 const Point & aLowerBound,
820 const Point & aUpperBound )
821{
822 unsigned int k;
823 bool in_here, in_further;
824 for ( k = 0; k < aKSpace.dimension; ++k )
825 {
826 Cell dir_low_uid = aKSpace.uSpel( aLowerBound );
827 Cell dir_up_uid = aKSpace.uGetDecr( aKSpace.uSpel( aUpperBound ), k);
828 Cell p = dir_low_uid;
829 do
830 {
831 in_here = pp( aKSpace.uCoords(p) );
832 in_further = pp( aKSpace.uCoords(aKSpace.uGetIncr( p, k )) );
833 if ( in_here != in_further ) // boundary element
834 { // add it to the set.
835 aBoundary.insert( aKSpace.uIncident( p, k, true ));
836 }
837 }
838 while ( aKSpace.uNext( p, dir_low_uid, dir_up_uid ) );
839 }
840}
841
842
843
844//-----------------------------------------------------------------------------
845template <typename TKSpace>
846template <typename SCellSet, typename PointPredicate >
847void
848DGtal::Surfaces<TKSpace>::
849sMakeBoundary( SCellSet & aBoundary,
850 const KSpace & aKSpace,
851 const PointPredicate & pp,
852 const Point & aLowerBound,
853 const Point & aUpperBound )
854{
855 unsigned int k;
856 bool in_here, in_further;
857
858 for ( k = 0; k < aKSpace.dimension; ++k )
859 {
860 Cell dir_low_uid = aKSpace.uSpel( aLowerBound );
861 Cell dir_up_uid = aKSpace.uGetDecr( aKSpace.uSpel( aUpperBound ), k);
862 Cell p = dir_low_uid;
863 do
864 {
865 in_here = pp( aKSpace.uCoords(p) );
866 in_further = pp( aKSpace.uCoords(aKSpace.uGetIncr( p, k )) );
867 if ( in_here != in_further ) // boundary element
868 { // add it to the set.
869 aBoundary.insert( aKSpace.sIncident( aKSpace.signs( p, in_here ),
870 k, true ));
871 }
872 }
873 while ( aKSpace.uNext( p, dir_low_uid, dir_up_uid ) );
874 }
875}
876
877
878//-----------------------------------------------------------------------------
879template <typename TKSpace>
880template <typename OutputIterator, typename PointPredicate >
881void
882DGtal::Surfaces<TKSpace>::
883uWriteBoundary( OutputIterator & out_it,
884 const KSpace & aKSpace,
885 const PointPredicate & pp,
886 const Point & aLowerBound, const Point & aUpperBound )
887{
888 unsigned int k;
889 bool in_here, in_further;
890 for ( k = 0; k < aKSpace.dimension; ++k )
891 {
892 Cell dir_low_uid = aKSpace.uSpel( aLowerBound );
893 Cell dir_up_uid = aKSpace.uGetDecr( aKSpace.uSpel( aUpperBound ), k);
894 Cell p = dir_low_uid;
895 do
896 {
897 in_here = pp( aKSpace.uCoords(p) );
898 in_further = pp( aKSpace.uCoords(aKSpace.uGetIncr( p, k )) );
899 if ( in_here != in_further ) // boundary element
900 { // writes it into the output iterator.
901 *out_it++ = aKSpace.uIncident( p, k, true );
902 }
903 }
904 while ( aKSpace.uNext( p, dir_low_uid, dir_up_uid ) );
905 }
906}
907
908
909
910//-----------------------------------------------------------------------------
911template <typename TKSpace>
912template <typename OutputIterator, typename PointPredicate >
913void
914DGtal::Surfaces<TKSpace>::
915sWriteBoundary( OutputIterator & out_it,
916 const KSpace & aKSpace,
917 const PointPredicate & pp,
918 const Point & aLowerBound, const Point & aUpperBound )
919{
920 // bool in_here, in_further;
921 bool in_here = false, in_before = false;
922
923 typedef typename KSpace::Space Space;
924 typedef HyperRectDomain<Space> Domain;
925 std::vector< Dimension > axes( aKSpace.dimension );
926 for ( Dimension k = 0; k < aKSpace.dimension; ++k )
927 axes[ k ] = k;
928
929 // We look for surfels in every direction.
930 for ( Dimension k = 0; k < aKSpace.dimension; ++k )
931 {
932 // When looking for surfels, the visiting must follow the k-th
933 // axis first so as to reuse the predicate "pp( p )".
934 std::swap( axes[ 0 ], axes[ k ] );
935
936 // We keep direct coordinates manipulation (instead of using KSpace methods)
937 // to allow correct domain span even with periodic Khalimsky space.
938 Point low = aLowerBound; ++low[ k ];
939 Point up = aUpperBound;
940 const Domain domain( low, up );
941 const Integer x = low[ k ];
942
943 for ( auto const& p : domain.subRange( axes ) )
944 {
945 auto cell = aKSpace.sSpel( p, true );
946
947 if ( p[ k ] == x)
948 {
949 in_here = pp( aKSpace.sCoords( cell ) );
950 in_before = pp( aKSpace.sCoords( aKSpace.sGetDecr( cell, k ) ) );
951 }
952 else
953 {
954 in_before = in_here;
955 in_here = pp( aKSpace.sCoords( cell ) );
956 }
957 if ( in_here != in_before ) // boundary element
958 { // writes it into the output iterator.
959 aKSpace.sSetSign( cell, in_here );
960 *out_it++ = aKSpace.sIncident( cell, k, false );
961 }
962 }
963 }
964}
965
966template <typename TKSpace>
967template <typename SurfelPredicate, typename TImageContainer>
968unsigned int
969DGtal::Surfaces<TKSpace>::uFillInterior( const KSpace & aKSpace,
970 const SurfelPredicate & aSurfPred,
971 TImageContainer & anImage,
972 const typename TImageContainer::Value & aValue,
973 bool empty_is_inside,
974 bool incrementMode)
975{
976 unsigned int numberCellFilled = 0;
977 std::deque<Cell> pixels;
978 Cell p = aKSpace.uFirst( typename KSpace::PreCell( Point::diagonal(1) ) );
979 const Cell lowerCell = aKSpace.uFirst( p );
980 const Cell upperCell = aKSpace.uLast( p );
981 bool in_mode = empty_is_inside;
982 do
983 {
984 if ( p != aKSpace.uGetMax( p, 0 ) )
985 {
986 pixels.push_front( p );
987 SCell b = aKSpace.sIncident( aKSpace.signs( p, KSpace::POS ), 0, true );
988 if ( aSurfPred( b ) )
989 {
990 while ( ! pixels.empty() )
991 {
992 Point aPoint = aKSpace.uCoords(pixels.front());
993 anImage.setValue(aPoint, aValue + (incrementMode ? anImage(aPoint) :0));
994 pixels.pop_front();
995 numberCellFilled++;
996 }
997 in_mode = false;
998 }
999 else if ( aSurfPred( aKSpace.sOpp( b ) ) )
1000 {
1001 pixels.clear();
1002 in_mode = true;
1003
1004 }
1005 }
1006 else
1007 {
1008 pixels.push_front( p );
1009 if ( in_mode )
1010 {
1011 while ( ! pixels.empty() )
1012 {
1013 Point aPoint = aKSpace.uCoords(pixels.front());
1014 anImage.setValue(aPoint, aValue + (incrementMode ? anImage(aPoint) :0));
1015 pixels.pop_front();
1016 numberCellFilled++;
1017 }
1018 }
1019 else
1020 pixels.clear();
1021 in_mode = empty_is_inside;
1022 }
1023 } while ( aKSpace.uNext(p, lowerCell, upperCell) );
1024
1025 return numberCellFilled;
1026}
1027
1028
1029
1030
1031template <typename TKSpace>
1032template <typename SurfelPredicate, typename TImageContainer>
1033unsigned int
1034DGtal::Surfaces<TKSpace>::uFillExterior( const KSpace & aKSpace,
1035 const SurfelPredicate & aSurfPred,
1036 TImageContainer & anImage,
1037 const typename TImageContainer::Value & aValue,
1038 bool empty_is_outside,
1039 bool incrementMode)
1040{
1041 unsigned int numberCellFilled=0;
1042 std::deque<Cell> pixels;
1043 Cell p = aKSpace.uFirst( typename KSpace::PreCell( Point::diagonal(1) ) );
1044 const Cell lowerCell = aKSpace.uFirst( p );
1045 const Cell upperCell = aKSpace.uLast( p );
1046 bool in_mode = empty_is_outside;
1047 do
1048 {
1049 if ( p != aKSpace.uGetMax( p, 0 ) )
1050 {
1051 pixels.push_front( p );
1052 SCell b = aKSpace.sIncident( aKSpace.signs( p, KSpace::POS ), 0, true );
1053 if ( aSurfPred( b ) )
1054 {
1055 pixels.clear();
1056 in_mode = true;
1057 }
1058 else if ( aSurfPred( aKSpace.sOpp( b ) ) )
1059 {
1060 while ( ! pixels.empty() )
1061 {
1062 Point aPoint = aKSpace.uCoords(pixels.front());
1063 anImage.setValue(aPoint, aValue + (incrementMode ? anImage(aPoint) :0));
1064 pixels.pop_front();
1065 numberCellFilled++;
1066 }
1067 in_mode = false;
1068 }
1069 }
1070 else
1071 {
1072 pixels.push_front( p );
1073 if ( in_mode )
1074 {
1075 while ( ! pixels.empty() )
1076 {
1077 Point aPoint = aKSpace.uCoords(pixels.front());
1078 anImage.setValue(aPoint, aValue + (incrementMode ? anImage(aPoint) :0));
1079 pixels.pop_front();
1080 numberCellFilled++;
1081 }
1082 }
1083 else
1084 pixels.clear();
1085 in_mode = empty_is_outside;
1086 }
1087 } while ( aKSpace.uNext(p, lowerCell, upperCell) );
1088 return numberCellFilled;
1089}
1090
1091template <typename TKSpace>
1092typename DGtal::Surfaces<TKSpace>::CellRange
1093DGtal::Surfaces<TKSpace>::getPrimalVertices( const KSpace& K, const SCell& s )
1094{
1095 auto faces = K.uFaces( K.unsigns( s ) );
1096 CellRange primal_vtcs;
1097 for ( auto&& f : faces ) {
1098 if ( K.uDim( f ) == 0 ) primal_vtcs.push_back( f );
1099 }
1100 return primal_vtcs;
1101}
1102
1103template <typename TKSpace>
1104typename DGtal::Surfaces<TKSpace>::CellRange
1105DGtal::Surfaces<TKSpace>::getPrimalVertices( const KSpace& K, const Surfel& s, bool ccw )
1106{
1107 BOOST_STATIC_ASSERT(( KSpace::dimension == 3 ));
1108 CellRange vtcs = getPrimalVertices( K, s );
1109 std::swap( vtcs[ 2 ], vtcs[ 3 ] );
1110 auto orth_dir = K.sOrthDir( s );
1111 auto direct = K.sDirect( s, orth_dir ) ? ccw : ! ccw;
1112 Vector s0s1 = K.uCoords( vtcs[ 1 ] ) - K.uCoords( vtcs[ 0 ] );
1113 Vector s0s2 = K.uCoords( vtcs[ 2 ] ) - K.uCoords( vtcs[ 0 ] );
1114 Vector t = s0s1.crossProduct( s0s2 );
1115 if ( ( ( t[ orth_dir ] > 0.0 ) && direct )
1116 || ( ( t[ orth_dir ] < 0.0 ) && ! direct ) )
1117 std::reverse( vtcs.begin(), vtcs.end() );
1118 return vtcs;
1119}
1120
1121
1122
1123
1124///////////////////////////////////////////////////////////////////////////////
1125// Interface - public :
1126
1127/**
1128 * Writes/Displays the object on an output stream.
1129 * @param out the output stream where the object is written.
1130 */
1131template <typename TKSpace>
1132inline
1133void
1134DGtal::Surfaces<TKSpace>::selfDisplay ( std::ostream & out ) const
1135{
1136 out << "[Surfaces]";
1137}
1138
1139/**
1140 * Checks the validity/consistency of the object.
1141 * @return 'true' if the object is valid, 'false' otherwise.
1142 */
1143template <typename TKSpace>
1144inline
1145bool
1146DGtal::Surfaces<TKSpace>::isValid() const
1147{
1148 return true;
1149}
1150
1151
1152
1153///////////////////////////////////////////////////////////////////////////////
1154// Implementation of inline functions //
1155
1156template <typename TKSpace>
1157inline
1158std::ostream&
1159DGtal::operator<< ( std::ostream & out,
1160 const Surfaces<TKSpace> & object )
1161{
1162 object.selfDisplay( out );
1163 return out;
1164}
1165
1166// //
1167///////////////////////////////////////////////////////////////////////////////
1168
1169
Generated on Sun Nov 27 2022 15:15:52 for DGtal by doxygen 1.9.5