DGtal  0.9.2
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  */
54 template <typename TKSpace>
55 inline
56 DGtal::Surfaces<TKSpace>::~Surfaces()
57 {
58 }
59 //-----------------------------------------------------------------------------
60 template <typename TKSpace>
61 template <typename PointPredicate>
62 typename DGtal::Surfaces<TKSpace>::SCell
63 DGtal::Surfaces<TKSpace>::
64 findABel
65 ( const KSpace & K,
66  const PointPredicate & pp,
67  unsigned int nbtries) throw (DGtal::InputException)
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 //-----------------------------------------------------------------------------
97 template <typename TKSpace>
98 template <typename PointPredicate>
99 typename DGtal::Surfaces<TKSpace>::SCell
100 DGtal::Surfaces<TKSpace>::
101 findABel
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 //-----------------------------------------------------------------------------
177 template <typename TKSpace>
178 template <typename SCellSet, typename PointPredicate >
179 void
180 DGtal::Surfaces<TKSpace>::
181 trackBoundary( 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 //-----------------------------------------------------------------------------
230 template <typename TKSpace>
231 template <typename SCellSet, typename SurfelPredicate >
232 void
233 DGtal::Surfaces<TKSpace>::
234 trackSurface( 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 //-----------------------------------------------------------------------------
284 template <typename TKSpace>
285 template <typename SCellSet, typename SurfelPredicate >
286 void
287 DGtal::Surfaces<TKSpace>::
288 trackClosedSurface( 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 //-----------------------------------------------------------------------------
331 template <typename TKSpace>
332 template <typename PointPredicate>
333 void
334 DGtal::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 //-----------------------------------------------------------------------------
403 template <typename TKSpace>
404 template <typename PointPredicate>
405 void
406 DGtal::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 //-----------------------------------------------------------------------------
474 template <typename TKSpace>
475 template <typename SurfelPredicate >
476 inline
477 void
478 DGtal::Surfaces<TKSpace>::
479 track2DSurface( 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 //-----------------------------------------------------------------------------
537 template <typename TKSpace>
538 template <typename SurfelPredicate >
539 inline
540 void
541 DGtal::Surfaces<TKSpace>::
542 track2DSliceSurface( 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 //-----------------------------------------------------------------------------
603 template <typename TKSpace>
604 template <typename PointPredicate>
605 inline
606 void
607 DGtal::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 //-----------------------------------------------------------------------------
635 template <typename TKSpace>
636 template <typename PointPredicate>
637 void
638 DGtal::Surfaces<TKSpace>::
639 extractAll2DSCellContours( 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 //-----------------------------------------------------------------------------
666 template <typename TKSpace>
667 template <typename PointPredicate>
668 void
669 DGtal::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 //-----------------------------------------------------------------------------
686 template <typename TKSpace>
687 template <typename PointPredicate>
688 void
689 DGtal::Surfaces<TKSpace>::
690 extractAllConnectedSCell
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 //-----------------------------------------------------------------------------
724 template <typename TKSpace>
725 template <typename PointPredicate>
726 void
727 DGtal::Surfaces<TKSpace>::
728 extractAllPointContours4C( 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 //-----------------------------------------------------------------------------
768 template <typename TKSpace>
769 template <typename SCellSet, typename PointPredicate >
770 void
771 DGtal::Surfaces<TKSpace>::
772 trackClosedBoundary( 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 //-----------------------------------------------------------------------------
812 template <typename TKSpace>
813 template <typename CellSet, typename PointPredicate >
814 void
815 DGtal::Surfaces<TKSpace>::
816 uMakeBoundary( 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 //-----------------------------------------------------------------------------
845 template <typename TKSpace>
846 template <typename SCellSet, typename PointPredicate >
847 void
848 DGtal::Surfaces<TKSpace>::
849 sMakeBoundary( 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 //-----------------------------------------------------------------------------
879 template <typename TKSpace>
880 template <typename OutputIterator, typename PointPredicate >
881 void
882 DGtal::Surfaces<TKSpace>::
883 uWriteBoundary( 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 //-----------------------------------------------------------------------------
911 template <typename TKSpace>
912 template <typename OutputIterator, typename PointPredicate >
913 void
914 DGtal::Surfaces<TKSpace>::
915 sWriteBoundary( 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 
966 template <typename TKSpace>
967 template <typename SurfelPredicate, typename TImageContainer>
968 unsigned int
969 DGtal::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 
1031 template <typename TKSpace>
1032 template <typename SurfelPredicate, typename TImageContainer>
1033 unsigned int
1034 DGtal::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 
1091 
1092 
1093 
1094 ///////////////////////////////////////////////////////////////////////////////
1095 // Interface - public :
1096 
1097 /**
1098  * Writes/Displays the object on an output stream.
1099  * @param out the output stream where the object is written.
1100  */
1101 template <typename TKSpace>
1102 inline
1103 void
1104 DGtal::Surfaces<TKSpace>::selfDisplay ( std::ostream & out ) const
1105 {
1106  out << "[Surfaces]";
1107 }
1108 
1109 /**
1110  * Checks the validity/consistency of the object.
1111  * @return 'true' if the object is valid, 'false' otherwise.
1112  */
1113 template <typename TKSpace>
1114 inline
1115 bool
1116 DGtal::Surfaces<TKSpace>::isValid() const
1117 {
1118  return true;
1119 }
1120 
1121 
1122 
1123 ///////////////////////////////////////////////////////////////////////////////
1124 // Implementation of inline functions //
1125 
1126 template <typename TKSpace>
1127 inline
1128 std::ostream&
1129 DGtal::operator<< ( std::ostream & out,
1130  const Surfaces<TKSpace> & object )
1131 {
1132  object.selfDisplay( out );
1133  return out;
1134 }
1135 
1136 // //
1137 ///////////////////////////////////////////////////////////////////////////////
1138 
1139