DGtal  1.5.beta
SymmetricConvexExpander.h
1 
17 #pragma once
18 
29 #if defined(SymmetricConvexExpander_RECURSES)
30 #error Recursive header files inclusion detected in SymmetricConvexExpander.h
31 #else // defined(SymmetricConvexExpander_RECURSES)
33 #define SymmetricConvexExpander_RECURSES
34 
35 #if !defined SymmetricConvexExpander_h
37 #define SymmetricConvexExpander_h
38 
40 // Inclusions
41 #include <iostream>
42 #include "DGtal/base/Common.h"
43 #include "DGtal/topology/CCellularGridSpaceND.h"
44 #include "DGtal/geometry/volumes/DigitalConvexity.h"
46 
47 namespace DGtal
48 {
49 
51  // class SymmetricConvexExpander
62  template < typename TKSpace,
63  typename TPointPredicate >
65  {
67 
68  public:
70  typedef TKSpace KSpace;
71  typedef TPointPredicate PointPredicate;
72  typedef typename KSpace::Integer Integer;
73  typedef typename KSpace::Point Point;
74  typedef typename KSpace::Vector Vector;
75  typedef typename KSpace::Space Space;
76  typedef std::size_t Size;
79  typedef std::vector<Point> PointRange;
80  typedef std::vector<Vector> VectorRange;
81  typedef std::unordered_set<Point> PointSet;
82 
83  static const Dimension dimension = KSpace::dimension;
84 
85  typedef std::pair< Point, Integer > Node;
86 
87  // Inversion order since priority queue output max element.
88  struct NodeComparator {
90  NodeComparator() = default;
91  // p < q iff p.second < q.second
92  bool operator()( const Node& p, const Node& q ) const
93  {
94  return p.second > q.second;
95  }
96  };
97 
98  typedef std::priority_queue< Node, std::vector<Node>, NodeComparator > NodeQueue;
99 
102  const Point& kcenter,
103  const Point& lo,
104  const Point& hi )
105  : myPredicate( &predicate ), myConvexity( lo, hi )
106  {
107  init( kcenter );
108  }
109 
110  bool predicate( const Point& p ) const
111  {
112  return (*myPredicate)( p );
113  }
114 
115  void init( const Point& kcenter )
116  {
117  myKCenter = kcenter;
118  myPoints.clear();
119  myQ = NodeQueue();
120  myM.clear();
121  myPerfectSymmetry = true;
123  // The starting points depend on the parity of the coordinates of the center.
124  // There are from 1 to 2^d starting points.
125  PointRange points;
126  const auto x = myKCenter[ 0 ];
127  if ( x % 2 == 0 )
128  points.push_back( Point::base( 0, x / 2 ) );
129  else
130  {
131  points.push_back( Point::base( 0, (x-1) / 2 ) );
132  points.push_back( Point::base( 0, (x+1) / 2 ) );
133  }
134  for ( Dimension k = 1; k < dimension; k++ )
135  {
136  const auto n = points.size();
137  const auto y = myKCenter[ k ];
138  if ( y % 2 == 0 )
139  {
140  for ( auto i = 0; i < n; i++ )
141  points[ i ][ k ] = y / 2;
142  }
143  else
144  {
145  points.resize( 2*n );
146  const auto z = (y-1)/2;
147  const auto z1 = z + 1;
148  for ( auto i = 0; i < n; i++ )
149  {
150  points[ i ][ k ] = z;
151  Point q = points[ i ];
152  q[ k ] = z1;
153  points[ i+n ] = q;
154  }
155  }
156  }
157  // Keep only the points that satisfy the predicate.
158  for ( auto&& p : points )
159  {
160  const Point sp = symmetric( p );
161  if ( ! myM.count( p )
162  && predicate( p ) && predicate( sp ) )
163  {
164  Node n( p, (2*p - myKCenter).squaredNorm() );
165  myQ.push( n );
166  myM.insert( p );
167  myM.insert( sp );
168  }
169  }
170  }
171 
173  bool advance( bool enforce_full_convexity )
174  {
175  while ( ! finished() )
176  {
177  const auto p = current().first;
178  //NOT USED const auto d = current().second; // current ring distance
179  const auto sp = symmetric( p );
180  myPoints.insert( p );
181  myPoints.insert( sp );
182  PointRange X( myPoints.cbegin(), myPoints.cend() );
183  if ( enforce_full_convexity && ! myConvexity.isFullyConvex( X ) )
184  {
185  myPoints.erase( p );
186  myPoints.erase( sp );
187  ignore();
188  }
189  else
190  {
191  expand();
192  return true;
193  }
194  }
195  return false;
196  }
197 
205  const Node& current() const
206  {
207  return myQ.top();
208  }
209 
217  void ignore()
218  {
219  myQ.pop();
220  }
221 
227  void expand()
228  {
229  const Point p = current().first;
230  myQ.pop();
231  myPoints.insert( p );
232  myPoints.insert( symmetric( p ) );
233  const auto next_points = next( p );
234  for ( auto&& q : next_points )
235  {
236  if ( ! myM.count( q ) )
237  {
238  const auto sq = symmetric( q );
239  const bool q_in = predicate( q );
240  const bool sq_in = predicate( sq );
241  if ( q_in && sq_in )
242  {
243  Node n( q, (2*q - myKCenter).squaredNorm() );
244  myQ.push( n );
245  myM.insert( q );
246  myM.insert( sq );
247  if ( myPerfectSymmetry )
249  n.second );
250  }
251  else if ( ( q_in && ! sq_in ) || ( ! q_in && sq_in ) )
252  {
253  myPerfectSymmetry = false;
254  }
255  }
256  }
257  }
258 
262  bool finished() const
263  {
264  return myQ.empty();
265  }
266 
267 
268  Point symmetric( const Point& p ) const
269  {
270  return myKCenter - p;
271  }
272 
273  PointRange next( const Point& p ) const
274  {
275  PointRange N;
276  Point d = 2*p - myKCenter;
277  for ( Dimension i = 0; i < dimension; i++ )
278  {
279  if ( d[ i ] >= 0 ) N.push_back( p + Point::base( i, 1 ) );
280  if ( d[ i ] <= 0 ) N.push_back( p - Point::base( i, 1 ) );
281  }
282  return N;
283  }
284 
287 
290 
293 
296 
302 
307 
308  };
309 
310 
311 } // namespace DGtal
312 
313 
314 // //
316 
317 #endif // !defined SymmetricConvexExpander_h
318 
319 #endif // !defined SymmetricConvexExpander_RECURSES
Aim: Represents an nD lattice polytope, i.e. a convex polyhedron bounded with vertices with integer c...
Aim: Represents an nD rational polytope, i.e. a convex polyhedron bounded by vertices with rational c...
Aim: A helper class to build polytopes from digital sets and to check digital k-convexity and full co...
bool isFullyConvex(const PointRange &X, bool convex0=false) const
TInteger Integer
Arithmetic ring induced by (+,-,*) and Integer numbers.
static Self base(Dimension k, Component val=1)
static Dimension size()
Aim: SymmetricConvexExpander computes symmetric fully convex subsets of a given digital set.
std::unordered_set< Point > PointSet
SymmetricConvexExpander(const PointPredicate &predicate, const Point &kcenter, const Point &lo, const Point &hi)
Constructor from predicate and symmetry center point.
BOOST_CONCEPT_ASSERT((concepts::CCellularGridSpaceND< TKSpace >))
DGtal::BoundedRationalPolytope< Space > RationalPolytope
DGtal::BoundedLatticePolytope< Space > LatticePolytope
PointRange next(const Point &p) const
std::pair< Point, Integer > Node
const PointPredicate * myPredicate
The predicate that every point must satisfy.
Point symmetric(const Point &p) const
PointSet myM
Marked points, i.e. points already in the queue or in the object.
bool predicate(const Point &p) const
Point myKCenter
Symmetry center (with doubled coordinates to represent half-integers).
bool advance(bool enforce_full_convexity)
Advance of one symmetric point.
std::priority_queue< Node, std::vector< Node >, NodeComparator > NodeQueue
DigitalConvexity< KSpace > myConvexity
The digital convexity object.
bool myPerfectSymmetry
True iff the set and its local complement are symmetric.
Integer myPerfectSymmetryRadius
Upper bound on the max distance of perfect symmetry.
DigitalConvexity< TKSpace > Self
PointSet myPoints
Symmetric range of lattice points, sorted.
DGtal is the top-level namespace which contains all DGtal functions and types.
DGtal::uint32_t Dimension
Definition: Common.h:136
NodeComparator()=default
Default constructor.
bool operator()(const Node &p, const Node &q) const
Aim: This concept describes a cellular grid space in nD. In these spaces obtained by cartesian produc...
int max(int a, int b)