DGtal  1.5.beta
LatticeSetByIntervals.h
1 
17 #pragma once
27 #if defined(LatticeSetByIntervals_RECURSES)
28 #error Recursive header files inclusion detected in LatticeSetByIntervals.h
29 #else // defined(LatticeSetByIntervals_RECURSES)
31 #define LatticeSetByIntervals_RECURSES
32 
33 #if !defined LatticeSetByIntervals_h
35 #define LatticeSetByIntervals_h
36 
37 #include <unordered_map>
38 #include <boost/iterator/iterator_facade.hpp>
39 #include <climits>
40 #include "DGtal/base/Common.h"
41 #include "DGtal/kernel/CUnsignedNumber.h"
42 #include "DGtal/kernel/CBoundedNumber.h"
43 #include "DGtal/kernel/CSpace.h"
44 #include "DGtal/kernel/PointVector.h"
45 #include "DGtal/kernel/IntegralIntervals.h"
46 
47 namespace DGtal
48 {
49 
51  // template class LatticeSetByIntervals
60  template < typename TSpace >
62  {
63  public:
65 
66  typedef TSpace Space;
68  using Point = typename Space::Point;
69  using Vector = typename Space::Vector;
70  using Integer = typename Space::Integer;
71  using PointRange= std::vector< Point >;
73  using Interval = typename Intervals::Interval;
74  using Container = std::map< Point, Intervals >;
75  using RowIterator = typename Container::iterator;
76  using RowConstIterator = typename Container::const_iterator;
77  using LatticeSetByInterval = std::map< Point, Interval >;
78  using Size = std::size_t;
79  using size_type = Size;
80  static const Dimension dimension = Space::dimension;
81 
82  //------------------- standard services (construction, move) -------------------
83  public:
86 
90  : myAxis( axis ), myData() {}
91 
94  LatticeSetByIntervals( const Self & other ) = default;
95 
98  LatticeSetByIntervals( Self&& other ) = default;
99 
103  Self& operator=( const Self & other ) = default;
104 
108  Self& operator=( Self&& other ) = default;
109 
114  template <typename PointIterator>
115  LatticeSetByIntervals( PointIterator it, PointIterator itE, Dimension axis = 0 )
116  : myAxis( axis )
117  {
118  for ( ; it != itE; ++it ) {
119  Point q = *it;
120  Integer x = q[ axis ];
121  q[ axis ] = 0;
122  myData[ q ].insert( x );
123  }
124  }
125 
133  : myAxis( axis )
134  {
135  for ( const auto& aRow : aSet )
136  myData[ aRow.first ].data().push_back( aRow.second );
137  }
138 
140  void clear() { myData.clear(); }
141 
148  {
149  myData.clear();
150  myAxis = axis;
151  }
152 
154  Dimension axis() const
155  { return myAxis; }
156 
158  Container& data() { return myData; }
159 
161 
162  //------------------- conversion services -----------------------------
163  public:
166 
169  {
170  PointRange X;
171  X.reserve( size() );
172  for ( const auto& pV : myData )
173  {
174  auto p = pV.first;
175  auto iVec = pV.second.integerVector();
176  for ( auto x : iVec )
177  {
178  p[ myAxis ] = x;
179  X.push_back( p );
180  }
181  }
182  return X;
183  }
184 
186 
187  //------------------- capacity services -----------------------------
188  public:
191 
194  bool empty() const
195  {
196  return myData.empty();
197  }
198 
202  Size size() const
203  {
204  Size nb = 0;
205  for ( const auto& pV : myData )
206  nb += pV.second.size();
207  return nb;
208  }
209 
213  Size max_size() const
214  {
215  return myData.max_size();
216  }
217 
219 
220  //------------------- modifier services -----------------------------
221  public:
224 
227  void insert( Point p )
228  {
229  Integer x = p[ myAxis ];
230  p[ myAxis ] = 0;
231  myData[ p ].insert( x );
232  }
233 
236  void erase( Point p )
237  {
238  Integer x = p[ myAxis ];
239  p[ myAxis ] = 0;
240  auto it = myData.find( p );
241  if ( it != myData.end() )
242  {
243  it->second.erase( x );
244  if ( it->second.empty() )
245  myData.erase( it ); // purge element if it was the last in the row.
246  }
247  }
248 
249 
251  void purge()
252  {
253  for ( auto it = myData.begin(), itE = myData.end(); it != itE; )
254  if ( it->second.empty() )
255  it = myData.erase( it );
256  else ++it;
257  }
258 
260 
261  //------------------- set operations --------------------------------
262  public:
265 
272  Self& add( const Self& other )
273  {
274  if ( other.axis() != axis() )
275  {
276  trace.error() << "[LatticeSetByInterval::add] "
277  << "Both lattice sets should share the same axis: "
278  << axis() << " != " << other.axis() << std::endl;
279  return *this;
280  }
281  for ( const auto& pV : other.myData )
282  {
283  const Point& p = pV.first;
284  auto it = myData.find( p );
285  if ( it != myData.end() )
286  it->second.add( pV.second );
287  else
288  myData[ p ] = pV.second;
289  }
290  return *this;
291  }
292 
299  Self& subtract( const Self& other )
300  {
301  if ( other.axis() != axis() )
302  {
303  trace.error() << "[LatticeSetByInterval::subtract] "
304  << "Both lattice sets should share the same axis: "
305  << axis() << " != " << other.axis() << std::endl;
306  return *this;
307  }
308  for ( const auto& pV : other.myData )
309  {
310  const Point& p = pV.first;
311  auto it = myData.find( p );
312  if ( it != myData.end() )
313  it->second.subtract( pV.second );
314  }
315  return *this;
316  }
317 
324  Self set_union( const Self& other ) const
325  {
326  Self U = *this;
327  U.add( other );
328  return U;
329  }
330 
337  Self set_difference( const Self& other ) const
338  {
339  Self U = *this;
340  U.subtract( other );
341  return U;
342  }
343 
350  Self set_intersection( const Self& other ) const
351  {
352  Self A_plus_B = set_union( other );
353  Self A_delta_B = set_symmetric_difference( other );
354  return A_plus_B.subtract( A_delta_B );
355  }
356 
363  Self set_symmetric_difference( const Self& other ) const
364  {
365  Self A_minus_B = *this;
366  A_minus_B.subtract( other );
367  Self B_minus_A = other;
368  B_minus_A.subtract( *this );
369  return A_minus_B.add( B_minus_A );
370  }
371 
377  bool includes( const Self& other ) const
378  {
379  if ( other.axis() != axis() )
380  {
381  trace.error() << "[LatticeSetByInterval::subtract] "
382  << "Both lattice sets should share the same axis: "
383  << axis() << " != " << other.axis() << std::endl;
384  return false;
385  }
386  for ( const auto& pV : other.myData )
387  {
388  const Point& p = pV.first;
389  const auto it = myData.find( p );
390  if ( it == myData.cend() ) return false;
391  if ( ! it->second.includes( pV.second ) ) return false;
392  }
393  return true;
394  }
395 
398  bool equals( const Self& other ) const
399  {
400  if ( other.axis() != axis() )
401  {
402  trace.error() << "[LatticeSetByInterval::subtract] "
403  << "Both lattice sets should share the same axis: "
404  << axis() << " != " << other.axis() << std::endl;
405  return false;
406  }
407  if ( myData.size() != other.myData.size() ) return false;
408  auto it = myData.cbegin();
409  for ( const auto& I : other.myData )
410  {
411  if ( it->first != I.first ) return false;
412  if ( ! it->second.equals( I.second ) ) return false;
413  ++it;
414  }
415  return true;
416  }
417 
419 
420  //------------------- topology operations --------------------------------
421  public:
424 
436  {
437  Self C( myAxis );
438  // First step, place points as pointels and insert their star along
439  // dimension a.
440  for ( auto& pV : myData )
441  {
442  const Point q = 2 * pV.first;
443  C.myData[ q ] = pV.second.starOfPoints();
444  }
445  // Second step, dilate along remaining directions
446  for ( Dimension k = 0; k < dimension; k++ )
447  {
448  if ( k == myAxis ) continue;
449  for ( const auto& value : C.myData )
450  {
451  Point q = value.first;
452  if ( q[ k ] & 0x1 ) continue;
453  q[ k ] -= 1;
454  C.myData[ q ].add( value.second );
455  q[ k ] += 2;
456  C.myData[ q ].add( value.second );
457  }
458  }
459  return C;
460  }
461 
468  {
469  Self C( *this );
470  // First step, compute star along dimension a.
471  for ( auto& pV : C.myData )
472  pV.second = pV.second.starOfCells();
473  // Second step, dilate along remaining directions
474  for ( Dimension k = 0; k < dimension; k++ )
475  {
476  if ( k == myAxis ) continue;
477  for ( const auto& value : C.myData )
478  {
479  Point q = value.first;
480  if ( q[ k ] & 0x1 ) continue;
481  q[ k ] -= 1;
482  C.myData[ q ].add( value.second );
483  q[ k ] += 2;
484  C.myData[ q ].add( value.second );
485  }
486  }
487  return C;
488  }
489 
496  {
497  Self S( *this );
498  // Now extracting implicitly its Skel
499  for ( const auto& pV : myData )
500  {
501  const Point& p = pV.first;
502  const auto& V = pV.second;
503  for ( Dimension k = 0; k < dimension; k++ )
504  {
505  if ( k == myAxis ) continue;
506  if ( ( p[ k ] & 0x1 ) != 0 ) continue; // if open along axis continue
507  // if closed, check upper incident cells along direction k
508  Point q = p; q[ k ] -= 1;
509  Point r = p; r[ k ] += 1;
510  auto itq = S.myData.find( q );
511  if ( itq != S.myData.end() )
512  {
513  auto& W = itq->second;
514  W.subtract( V );
515  }
516  auto itr = S.myData.find( r );
517  if ( itr != S.myData.end() )
518  {
519  auto& W = itr->second;
520  W.subtract( V );
521  }
522  }
523  }
524  // Extract skel along main axis
525  for ( auto& value : S.myData )
526  {
527  auto & V = value.second;
528  Intervals sub_to_V;
529  for ( auto I : V.data() )
530  {
531  if ( ( I.first & 0x1 ) != 0 )
532  {
533  if ( I.first != I.second )
534  sub_to_V.data().push_back( Interval{ I.first, I.first } );
535  I.first += 1;
536  }
537  if ( ( I.second & 0x1 ) != 0 ) I.second -= 1;
538  for ( auto x = I.first; x <= I.second; x += 2 )
539  sub_to_V.data().push_back( Interval{ x+1, x+1 } );
540  }
541  V.subtract( sub_to_V );
542  }
543  // Erase empty stacks
544  S.purge();
545  return S;
546  }
547 
551  {
552  typedef std::vector< Integer > Coordinates;
553  std::map< Point, Coordinates > E;
554  // Now extracting vertices along axis direction.
555  for ( const auto& pV : myData )
556  {
557  const Point& p = pV.first;
558  const auto& V = pV.second;
559  E[ p ] = V.extremaOfCells();
560  }
561  std::map< Point, Coordinates > next_E;
562  for ( Dimension k = 0; k != dimension; k++ )
563  {
564  if ( k == myAxis ) continue;
565  for ( const auto& pC : E )
566  {
567  const auto& p = pC.first;
568  Point q = p;
569  q[ k ] = q[ k ] >> 1;
570  bool odd = ( p[ k ] & 0x1 ) != 0;
571  { // odd/even always copy
572  auto it = next_E.find( q );
573  if ( it == next_E.end() ) next_E[ q ] = pC.second;
574  else
575  {
576  Coordinates F;
577  std::set_union( pC.second.cbegin(), pC.second.cend(),
578  it->second.cbegin(), it->second.cend(),
579  std::back_inserter( F ) );
580  it->second = F;
581  }
582  }
583  if ( odd )
584  { // odd: must copy forward also
585  q[ k ] += 1;
586  auto it = next_E.find( q );
587  if ( it == next_E.end() ) next_E[ q ] = pC.second;
588  else
589  {
590  Coordinates F;
591  std::set_union( pC.second.cbegin(), pC.second.cend(),
592  it->second.cbegin(), it->second.cend(),
593  std::back_inserter( F ) );
594  it->second = F;
595  }
596  }
597  } // for ( const auto& pC : E )
598  E.swap( next_E );
599  next_E.clear();
600  }
601  // Build point range.
602  PointRange R;
603  for ( const auto& pC : E )
604  {
605  Point p = pC.first;
606  for ( auto&& x : pC.second )
607  {
608  p[ myAxis ] = x;
609  R.push_back( p );
610  }
611  }
612  std::sort( R.begin(), R.end() );
613  return R;
614  }
615 
617 
618  //------------------- specific services (interval insertion, removal) -------------
619  public:
622 
625  size_type memory_usage() const noexcept
626  {
627  size_type nb = 0;
628  for ( const auto& pV : myData )
629  {
630  nb += sizeof( Interval ) * pV.second.capacity() + sizeof( void* );
631  nb += sizeof( pV ) + sizeof( void* );
632  }
633  nb += sizeof( Self );
634  return nb;
635  }
636 
643  {
644  q[ myAxis ] = 0;
645  return myData[ q ];
646  }
647 
653  void reset( Point p )
654  {
655  p[ myAxis ] = 0;
656  auto it = myData.find( p );
657  if ( it != myData.end() ) myData.erase( it );
658  }
659 
660 
665  };
666 
667 } // namespace DGtal
668 
669 
671 // Includes inline functions.
672 
673 // //
675 
676 #endif // !defined LatticeSetByIntervals_h
677 
678 #undef LatticeSetByIntervals_RECURSES
679 #endif // else defined(LatticeSetByIntervals_RECURSES)
std::pair< Integer, Integer > Interval
void clear()
Clears the data structure.
size_type memory_usage() const noexcept
LatticeSetByIntervals(const LatticeSetByInterval &aSet, Dimension axis)
typename Space::Integer Integer
typename Container::iterator RowIterator
typename Intervals::Interval Interval
typename Container::const_iterator RowConstIterator
std::map< Point, Intervals > Container
LatticeSetByIntervals(const Self &other)=default
Self set_difference(const Self &other) const
BOOST_CONCEPT_ASSERT((concepts::CSpace< TSpace >))
Self set_union(const Self &other) const
Container myData
Associate to each point its sequences of intervals.
Self set_symmetric_difference(const Self &other) const
Self & operator=(Self &&other)=default
LatticeSetByIntervals(Self &&other)=default
std::map< Point, Interval > LatticeSetByInterval
bool includes(const Self &other) const
Self & add(const Self &other)
Self set_intersection(const Self &other) const
LatticeSetByIntervals(PointIterator it, PointIterator itE, Dimension axis=0)
Self & operator=(const Self &other)=default
bool equals(const Self &other) const
static const Dimension dimension
Dimension myAxis
The axis along which data is stacked in intervals.
void purge()
Eliminates rows that contains no element.
Self & subtract(const Self &other)
std::ostream & error()
std::vector< Point > PointRange
DigitalPlane::Point Vector
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
Trace trace
Definition: Common.h:153
Aim: Defines the concept describing a digital space, ie a cartesian product of integer lines.
Definition: CSpace.h:106
MyPointD Point
Definition: testClone2.cpp:383