DGtal  1.5.beta
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  */
53 template <typename TKSpace>
54 inline
55 DGtal::Surfaces<TKSpace>::~Surfaces()
56 {
57 }
58 //-----------------------------------------------------------------------------
59 template <typename TKSpace>
60 template <typename PointPredicate>
61 typename DGtal::Surfaces<TKSpace>::SCell
62 DGtal::Surfaces<TKSpace>::
63 findABel
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 //-----------------------------------------------------------------------------
96 template <typename TKSpace>
97 template <typename PointPredicate>
98 typename DGtal::Surfaces<TKSpace>::SCell
99 DGtal::Surfaces<TKSpace>::
100 findABel
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 //-----------------------------------------------------------------------------
176 template <typename TKSpace>
177 template <typename SCellSet, typename PointPredicate >
178 void
179 DGtal::Surfaces<TKSpace>::
180 trackBoundary( 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 //-----------------------------------------------------------------------------
229 template <typename TKSpace>
230 template <typename SCellSet, typename SurfelPredicate >
231 void
232 DGtal::Surfaces<TKSpace>::
233 trackSurface( 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 //-----------------------------------------------------------------------------
283 template <typename TKSpace>
284 template <typename SCellSet, typename SurfelPredicate >
285 void
286 DGtal::Surfaces<TKSpace>::
287 trackClosedSurface( 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 //-----------------------------------------------------------------------------
330 template <typename TKSpace>
331 template <typename PointPredicate>
332 void
333 DGtal::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 //-----------------------------------------------------------------------------
402 template <typename TKSpace>
403 template <typename PointPredicate>
404 void
405 DGtal::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 //-----------------------------------------------------------------------------
473 template <typename TKSpace>
474 template <typename SurfelPredicate >
475 inline
476 void
477 DGtal::Surfaces<TKSpace>::
478 track2DSurface( 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 //-----------------------------------------------------------------------------
536 template <typename TKSpace>
537 template <typename SurfelPredicate >
538 inline
539 void
540 DGtal::Surfaces<TKSpace>::
541 track2DSliceSurface( 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 //-----------------------------------------------------------------------------
602 template <typename TKSpace>
603 template <typename PointPredicate>
604 inline
605 void
606 DGtal::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 //-----------------------------------------------------------------------------
634 template <typename TKSpace>
635 template <typename PointPredicate>
636 void
637 DGtal::Surfaces<TKSpace>::
638 extractAll2DSCellContours( 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 //-----------------------------------------------------------------------------
665 template <typename TKSpace>
666 template <typename PointPredicate>
667 void
668 DGtal::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 //-----------------------------------------------------------------------------
685 template <typename TKSpace>
686 template <typename PointPredicate>
687 void
688 DGtal::Surfaces<TKSpace>::
689 extractAllConnectedSCell
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 //-----------------------------------------------------------------------------
723 template <typename TKSpace>
724 template <typename PointPredicate>
725 void
726 DGtal::Surfaces<TKSpace>::
727 extractAllPointContours4C( 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 //-----------------------------------------------------------------------------
767 template <typename TKSpace>
768 template <typename SCellSet, typename PointPredicate >
769 void
770 DGtal::Surfaces<TKSpace>::
771 trackClosedBoundary( 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 //-----------------------------------------------------------------------------
811 template <typename TKSpace>
812 template <typename CellSet, typename PointPredicate >
813 void
814 DGtal::Surfaces<TKSpace>::
815 uMakeBoundary( 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 //-----------------------------------------------------------------------------
844 template <typename TKSpace>
845 template <typename SCellSet, typename PointPredicate >
846 void
847 DGtal::Surfaces<TKSpace>::
848 sMakeBoundary( 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 //-----------------------------------------------------------------------------
878 template <typename TKSpace>
879 template <typename OutputIterator, typename PointPredicate >
880 void
881 DGtal::Surfaces<TKSpace>::
882 uWriteBoundary( 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 //-----------------------------------------------------------------------------
910 template <typename TKSpace>
911 template <typename OutputIterator, typename PointPredicate >
912 void
913 DGtal::Surfaces<TKSpace>::
914 sWriteBoundary( 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 
965 template <typename TKSpace>
966 template <typename SurfelPredicate, typename TImageContainer>
967 unsigned int
968 DGtal::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 
1030 template <typename TKSpace>
1031 template <typename SurfelPredicate, typename TImageContainer>
1032 unsigned int
1033 DGtal::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 
1090 template <typename TKSpace>
1091 typename DGtal::Surfaces<TKSpace>::CellRange
1092 DGtal::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 
1102 template <typename TKSpace>
1103 typename DGtal::Surfaces<TKSpace>::CellRange
1104 DGtal::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  */
1130 template <typename TKSpace>
1131 inline
1132 void
1133 DGtal::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  */
1142 template <typename TKSpace>
1143 inline
1144 bool
1145 DGtal::Surfaces<TKSpace>::isValid() const
1146 {
1147  return true;
1148 }
1149 
1150 
1151 
1152 ///////////////////////////////////////////////////////////////////////////////
1153 // Implementation of inline functions //
1154 
1155 template <typename TKSpace>
1156 inline
1157 std::ostream&
1158 DGtal::operator<< ( std::ostream & out,
1159  const Surfaces<TKSpace> & object )
1160 {
1161  object.selfDisplay( out );
1162  return out;
1163 }
1164 
1165 // //
1166 ///////////////////////////////////////////////////////////////////////////////
1167 
1168