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