DGtal  1.0.0
testDigitalSurfaceBoostGraphInterface.cpp
Go to the documentation of this file.
1 
30 #include <iostream>
32 #include <boost/property_map/property_map.hpp>
33 #include <boost/pending/queue.hpp>
34 #include "DGtal/base/Common.h"
35 #include "DGtal/math/MPolynomial.h"
36 #include "DGtal/graph/DigitalSurfaceBoostGraphInterface.h"
37 #include "DGtal/helpers/StdDefs.h"
38 #include "DGtal/topology/helpers/Surfaces.h"
39 #include "DGtal/topology/DigitalSurface.h"
40 #include "DGtal/topology/SetOfSurfels.h"
41 #include "DGtal/shapes/GaussDigitizer.h"
42 #include "DGtal/shapes/implicit/ImplicitPolynomial3Shape.h"
43 #include "DGtal/io/readers/MPolynomialReader.h"
47 #include <boost/graph/graph_concepts.hpp>
48 #include <boost/graph/adjacency_list.hpp>
49 #include <boost/graph/copy.hpp>
50 #include <boost/graph/breadth_first_search.hpp>
51 #include <boost/graph/connected_components.hpp>
52 #include <boost/graph/stoer_wagner_min_cut.hpp>
53 #include <boost/graph/boykov_kolmogorov_max_flow.hpp>
55 
56 using namespace std;
57 using namespace DGtal;
58 
60 // Functions for testing class DigitalSurfaceBoostGraphInterface.
62 
63 
64 struct surfel_position_t {
65  typedef boost::vertex_property_tag kind;
66 };
67 
68 struct surfel_position {
69  Z3i::Point myP;
70  surfel_position() : myP()
71  {
72  }
73 };
74 
75 typedef boost::property< boost::vertex_index_t, std::size_t,
76  boost::property<surfel_position_t, surfel_position> > VertexProperties;
77 
78 
79 template <typename Graph1, typename Graph2, typename VertexIndexMap>
80 struct my_vertex_copier {
81  typedef typename boost::property_map< Graph2, boost::vertex_index_t>::type graph_vertex_index_map;
82  typedef typename boost::property_map< Graph2, surfel_position_t>::type graph_vertex_position_map;
83  typedef typename boost::graph_traits< Graph1 >::vertex_descriptor Vertex1;
84  typedef typename boost::graph_traits< Graph2 >::vertex_descriptor Vertex2;
85 
86  const Graph1 & myG1;
87  graph_vertex_index_map graph_vertex_index;
88  graph_vertex_position_map graph_vertex_position;
89  VertexIndexMap & myIndexMap;
90 
91  my_vertex_copier(const Graph1& g1, Graph2& g2, VertexIndexMap & indexMap )
92  : myG1( g1 ),
93  graph_vertex_index( boost::get( boost::vertex_index_t(), g2 ) ),
94  graph_vertex_position( boost::get( surfel_position_t(), g2) ),
95  myIndexMap( indexMap )
96  {}
97 
98  void operator()( const Vertex1& v1, const Vertex2& v2 ) const {
99  // std::size_t idx = myIndexMap[ v1 ];
100  // Does not work !
101  // put( graph_vertex_index, v2, idx);
102  surfel_position pos;
103  pos.myP = v1.myCoordinates;
104  //std::cout << "surfel " << idx << " at " << pos.myP << std::endl;
105  put( graph_vertex_position, v2, pos);
106  }
107 };
108 template <typename Graph1, typename Graph2>
109 struct my_edge_copier {
110  my_edge_copier(const Graph1& UNUSED(g1), Graph2& UNUSED(g2))
111  {}
112  template <typename Edge1, typename Edge2>
113  void operator()(const Edge1& /*v1*/, Edge2& /*v2*/) const {
114  // does nothing
115  }
116 };
117 
123 {
124  unsigned int nbok = 0;
125  unsigned int nb = 0;
126 
127 
128  using namespace Z3i;
129 
130  // Construct the implicit shape
131  typedef Space::RealPoint RealPoint;
132  typedef RealPoint::Coordinate Ring;
137 
138  // Implicit shape is an ellipse
139  trace.beginBlock( "Constructing implicit shape." );
140  double p1[] = {-2,-2,-2};
141  double p2[] = { 2, 2, 2};
142  std::string poly_str = "x*x+y*y+z*z-1";
143  double step = 0.4;
144  Polynomial3 P;
145  Polynomial3Reader reader;
146  std::string::const_iterator iter
147  = reader.read( P, poly_str.begin(), poly_str.end() );
148  if ( iter != poly_str.end() )
149  {
150  std::cerr << "ERROR: I read only <"
151  << poly_str.substr( 0, iter - poly_str.begin() )
152  << ">, and I built P=" << P << std::endl;
153  return 1;
154  }
155  trace.info() << "- P = " << P << std::endl;
156  ImplicitShape ishape( P );
157  DigitalShape dshape;
158  dshape.attach( ishape );
159  dshape.init( RealPoint( p1 ), RealPoint( p2 ), step );
160  Domain domain = dshape.getDomain();
161  trace.endBlock();
162 
163  // Construct the Khalimsky space from the image domain
164  trace.beginBlock ( "Construct the Khalimsky space from the image domain ..." );
165  KSpace K;
166  bool space_ok = K.init( domain.lowerBound(),
167  domain.upperBound(), true // necessary
168  );
169  if (!space_ok)
170  {
171  trace.error() << "Error in the Khamisky space construction."<<std::endl;
172  return false;
173  }
174  trace.endBlock();
175 
177  typedef SurfelAdjacency<KSpace::dimension> MySurfelAdjacency;
178  typedef KSpace::Surfel Surfel;
180  typedef SetOfSurfels< KSpace, SurfelSet > MySetOfSurfels;
182  trace.beginBlock ( "Extract surface ..." );
183 
184  MySurfelAdjacency surfAdj( true ); // interior in all directions.
185  MySetOfSurfels theSetOfSurfels( K, surfAdj );
186  Surfel bel = Surfaces<KSpace>::findABel( K, dshape, 100000 );
187  Surfaces<KSpace>::trackBoundary( theSetOfSurfels.surfelSet(),
188  K, surfAdj,
189  dshape, bel );
190  MyDigitalSurface digSurf( theSetOfSurfels );
191  trace.info() << "Digital surface has " << digSurf.size() << " surfels."
192  << std::endl;
193  trace.endBlock();
194 
195  trace.beginBlock ( "Testing Graph concepts for DigitalSurface ..." );
196  typedef MyDigitalSurface Graph;
197  typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor; // ie DigitalSurface::Vertex
198  typedef boost::graph_traits<Graph>::edge_descriptor edge_descriptor; // ie DigitalSurface::Arc
199  typedef boost::graph_traits<Graph>::vertices_size_type vertices_size_type; // ie DigitalSurface::Size
200  typedef boost::graph_traits<Graph>::vertex_iterator vertex_iterator;
201  typedef boost::graph_traits<Graph>::out_edge_iterator out_edge_iterator;
202  typedef boost::graph_traits<Graph>::edge_iterator edge_iterator;
203 
204  BOOST_CONCEPT_ASSERT(( boost::VertexListGraphConcept<Graph> ));
205  BOOST_CONCEPT_ASSERT(( boost::AdjacencyGraphConcept<Graph> ));
206  BOOST_CONCEPT_ASSERT(( boost::IncidenceGraphConcept<Graph> ));
207  BOOST_CONCEPT_ASSERT(( boost::EdgeListGraphConcept<Graph> ));
208  // BOOST_CONCEPT_ASSERT(( BidirectionalGraphConcept<Graph> ));
209  trace.endBlock();
210 
211  trace.beginBlock ( "Testing IncidenceGraph interface with breadth_first_visit ..." );
212  // get the property map for coloring vertices.
213  typedef std::map< vertex_descriptor, boost::default_color_type > StdColorMap;
214  StdColorMap colorMap;
215  boost::associative_property_map< StdColorMap > propColorMap( colorMap );
216  // get the property map for storing distances
217  typedef std::map< vertex_descriptor, uint64_t > StdDistanceMap;
218  StdDistanceMap distanceMap;
219  boost::associative_property_map< StdDistanceMap > propDistanceMap( distanceMap );
220  boost::queue< vertex_descriptor > Q; // std::queue does not have top().
221  vertex_descriptor start = *( digSurf.begin() );
222  boost::breadth_first_visit // boost graph breadth first visiting algorithm.
223  ( digSurf, // the graph
224  start, // the starting vertex
225  Q, // the buffer for breadth first queueing
226  boost::make_bfs_visitor( boost::record_distances( propDistanceMap, boost::on_tree_edge() ) ), // only record distances
227  propColorMap // necessary for the visiting vertices
228  );
229  trace.endBlock();
230 
231  // Display results
232  trace.beginBlock ( "Display breadth_first_visit result ..." );
233  uint64_t maxD = 0;
234  vertex_descriptor furthest = start;
235  uint64_t nbV = 0;
236  for ( std::pair<vertex_iterator, vertex_iterator>
237  vp = boost::vertices( digSurf ); vp.first != vp.second; ++vp.first, ++nbV )
238  {
239  uint64_t d = boost::get( propDistanceMap, *vp.first );
240  if ( d > maxD )
241  {
242  maxD = d;
243  furthest = *vp.first;
244  }
245  }
246  trace.info() << "- d[ " << start << " ] = " << boost::get( propDistanceMap, start ) << std::endl;
247  trace.info() << "- d[ " << furthest << " ] = " << maxD << std::endl;
248  ++nb; nbok += ( nbV == digSurf.size() ) ? 1 : 0;
249  trace.info() << "(" << nbok << "/" << nb << ") "
250  << "nb vertices is ok" << std::endl;
251  ++nb; nbok += ( maxD == 12 ) ? 1 : 0;
252  trace.info() << "(" << nbok << "/" << nb << ") "
253  << "maxD == 12" << std::endl;
254  trace.endBlock();
255 
256  trace.beginBlock ( "Testing VertexListGraph interface with connected_components ..." );
257  // get the property map for labelling vertices.
258  typedef std::map< vertex_descriptor, vertices_size_type > StdComponentMap;
259  StdComponentMap componentMap;
260  boost::associative_property_map< StdComponentMap > propComponentMap( componentMap );
261  vertices_size_type nbComp =
262  boost::connected_components // boost graph connected components algorithm.
263  ( digSurf, // the graph
264  propComponentMap, // the mapping vertex -> label
265  boost::color_map( propColorMap ) // this map is used internally when computing connected components.
266  );
267  trace.info() << "- nbComponents = " << nbComp << std::endl;
268  ++nb; nbok += ( nbComp == 1 ) ? 1 : 0;
269  trace.info() << "(" << nbok << "/" << nb << ") "
270  << "nbComp == 1" << std::endl;
271  trace.endBlock();
272 
273  trace.beginBlock ( "Testing UndirectedGraph interface with Wagner-Stoer min cut ..." );
274  // get the property map for weighting edges.
275  typedef double weight_type;
276  typedef std::map< edge_descriptor, weight_type > StdWeightMap;
277  StdWeightMap weightMap;
278  boost::associative_property_map< StdWeightMap > propWeightMap( weightMap );
279  typedef std::map< vertex_descriptor, vertices_size_type > StdVertexIndexMap;
280  StdVertexIndexMap vertexIndexMap;
281  boost::associative_property_map< StdVertexIndexMap > propVertexIndexMap( vertexIndexMap );
282  vertices_size_type idxV = 0;
283  // The weight is smaller for edges traversing plane z=0 than anywhere else.
284  // The min cut thus cuts the sphere in two approximate halves.
285  for ( std::pair<vertex_iterator, vertex_iterator>
286  vp = boost::vertices( digSurf ); vp.first != vp.second; ++vp.first, ++idxV )
287  {
288  vertex_descriptor v1 = *vp.first;
289  vertexIndexMap[ v1 ] = idxV;
290  for ( std::pair<out_edge_iterator, out_edge_iterator>
291  ve = boost::out_edges( v1, digSurf ); ve.first != ve.second; ++ve.first )
292  {
293  vertex_descriptor v2 = boost::target( *ve.first, digSurf );
294  if ( v1 < v2 )
295  {
296  KSpace::SCell sep = digSurf.separator( *ve.first );
297  weight_type weight = ( K.sKCoord( sep, 2 ) == 0 ) ? 0.01 : 1.0;
298  weightMap[ *ve.first ] = weight;
299  weightMap[ digSurf.opposite( *ve.first ) ] = weight;
300  }
301  }
302  }
303  // get the parity map for assigning a set to each vertex.
304  typedef std::map< vertex_descriptor, bool > StdParityMap;
305  StdParityMap parityMap;
306  boost::associative_property_map< StdParityMap > propParityMap( parityMap );
307 
308  weight_type total_weight =
309  boost::stoer_wagner_min_cut // boost wagner stoer min cut algorithm.
310  ( digSurf, // the graph
311  propWeightMap, // the mapping edge -> weight
312  boost::parity_map( propParityMap ) // this map stores the vertex assignation
313  .vertex_index_map( propVertexIndexMap )
314  );
315  trace.info() << "- total weight = " << total_weight << std::endl;
316  uint64_t nb0 = 0;
317  uint64_t nb1 = 0;
318  for ( std::pair<vertex_iterator, vertex_iterator>
319  vp = boost::vertices( digSurf ); vp.first != vp.second; ++vp.first, ++idxV )
320  {
321  vertex_descriptor v1 = *vp.first;
322  //trace.info() << "- " << v1 << " in " << parityMap[ v1 ] << std::endl;
323  if ( parityMap[ v1 ] ) ++nb1;
324  else ++nb0;
325  }
326  ++nb; nbok += ( total_weight < 1.0 ) ? 1 : 0;
327  trace.info() << "(" << nbok << "/" << nb << ") "
328  << "total_weight < 1.0"
329  << ", nb0=" << nb0 << " nb1=" << nb1 << std::endl;
330  trace.info() << "- parity[ " << start << " ] = " << parityMap[ start ] << std::endl;
331  trace.info() << "- parity[ " << furthest << " ] = " << parityMap[ furthest ] << std::endl;
332  ++nb; nbok += ( parityMap[ start ] != parityMap[ furthest ] ) ? 1 : 0;
333  trace.info() << "(" << nbok << "/" << nb << ") "
334  << "parityMap[ start ] != parityMap[ furthest ]" << std::endl;
335  trace.endBlock();
336 
337  trace.beginBlock ( "Testing EdgeListGraph and IncidenceGraph interfaces with Boykov-Kolmogorov max flow ..." );
338  typedef double capacity_type;
339  // get the property map for edge capacity.
340  typedef std::map< edge_descriptor, weight_type > StdEdgeCapacityMap;
341  StdEdgeCapacityMap edgeCapacityMap;
342  boost::associative_property_map< StdEdgeCapacityMap > propEdgeCapacityMap( edgeCapacityMap );
343  // get the property map for edge residual capacity.
344  typedef std::map< edge_descriptor, weight_type > StdEdgeResidualCapacityMap;
345  StdEdgeResidualCapacityMap edgeResidualCapacityMap;
346  boost::associative_property_map< StdEdgeResidualCapacityMap > propEdgeResidualCapacityMap( edgeResidualCapacityMap );
347  // get the property map for reverse edge.
348  typedef std::map< edge_descriptor, edge_descriptor > StdReversedEdgeMap;
349  StdReversedEdgeMap reversedEdgeMap;
350  boost::associative_property_map< StdReversedEdgeMap > propReversedEdgeMap( reversedEdgeMap );
351  // get the property map for vertex predecessor.
352  typedef std::map< vertex_descriptor, edge_descriptor > StdPredecessorMap;
353  StdPredecessorMap predecessorMap;
354  boost::associative_property_map< StdPredecessorMap > propPredecessorMap( predecessorMap );
355  // We already have vertex color map, vertex distance map and vertex index map.
356  uint64_t nbEdges = 0;
357  // The weight is smaller for edges traversing plane z=0 than anywhere else.
358  // The min cut thus cuts the sphere in two approximate halves.
359  for ( std::pair<edge_iterator, edge_iterator>
360  ve = boost::edges( digSurf ); ve.first != ve.second; ++ve.first, ++nbEdges )
361  {
362  edge_descriptor e = *ve.first;
363  edge_descriptor rev_e = digSurf.opposite( e );
364  vertex_descriptor v1 = boost::source( e, digSurf );
365  vertex_descriptor v2 = boost::target( e, digSurf );
366  ASSERT( boost::source( rev_e, digSurf ) == v2 );
367  ASSERT( boost::target( rev_e, digSurf ) == v1 );
368  if ( v1 < v2 )
369  {
370  KSpace::SCell sep = digSurf.separator( *ve.first );
371  capacity_type capacity = ( K.sKCoord( sep, 2 ) == 0 ) ? 0.01 : 1.0;
372  edgeCapacityMap[ e ] = capacity;
373  edgeCapacityMap[ digSurf.opposite( e ) ] = capacity;
374  reversedEdgeMap[ e ] = digSurf.opposite( e );
375  reversedEdgeMap[ digSurf.opposite( e ) ] = e;
376  }
377  }
378  trace.info() << "- nb edges = " << nbEdges << std::endl;
379  distanceMap.clear();
380  colorMap.clear();
381  capacity_type max_flow =
382  boost::boykov_kolmogorov_max_flow // boykov kolmogorov max flow algorithm.
383  ( digSurf, // the graph
384  propEdgeCapacityMap, propEdgeResidualCapacityMap,
385  propReversedEdgeMap, propPredecessorMap, propColorMap, propDistanceMap, propVertexIndexMap,
386  start, furthest );
387  trace.info() << "- max flow = " << max_flow << std::endl;
388  ++nb; nbok += ( abs( max_flow - total_weight ) < 0.0000001 ) ? 1 : 0;
389  trace.info() << "(" << nbok << "/" << nb << ") "
390  << "max_flow == min_cut, Duality max-flow/min-cut." << std::endl;
391  trace.endBlock();
392 
393  trace.beginBlock ( "Testing AdjacencyListGraph with copy_graph ..." );
394  typedef boost::adjacency_list< boost::vecS, boost::vecS, boost::undirectedS, VertexProperties > BGraph;
395  BGraph bG;
396  boost::copy_graph( digSurf, bG,
397  boost::vertex_copy( my_vertex_copier<Graph,BGraph,StdVertexIndexMap>( digSurf, bG, vertexIndexMap ) )
398  .edge_copy( my_edge_copier<Graph,BGraph>( digSurf, bG ) )
399  .vertex_index_map( propVertexIndexMap )
400  );
401  typedef boost::graph_traits< BGraph >::vertex_descriptor vertex_descriptor_2;
402  typedef boost::graph_traits< BGraph >::vertex_iterator vertex_iterator_2;
403  typedef boost::property_map< BGraph, surfel_position_t>::type GraphSurfelPositionMap;
404  GraphSurfelPositionMap surfelPos( boost::get( surfel_position_t(), bG) );
405  for ( std::pair<vertex_iterator_2, vertex_iterator_2>
406  vp = boost::vertices( bG ); vp.first != vp.second; ++vp.first )
407  {
408  vertex_descriptor_2 v1 = *vp.first;
409  surfel_position pos = boost::get( surfelPos, v1 );
410  trace.info() << "- " << v1 << " was at " << pos.myP << std::endl;
411  }
412  ++nb; nbok += boost::num_vertices( bG ) == boost::num_vertices( digSurf ) ? 1 : 0;
413  trace.info() << "(" << nbok << "/" << nb << ") "
414  << "after copy: Boost graph has " << num_vertices( bG )
415  << " vertices." << std::endl;
416  trace.endBlock();
417 
418  return nbok == nb;
419 }
420 
421 
423 // Standard services - public :
424 
425 int main( int /* argc */, char** /* argv */ )
426 {
427  trace.beginBlock ( "Testing class DigitalSurfaceBoostGraphInterface" );
428 
429  bool res = testDigitalSurfaceBoostGraphInterface(); // && ... other tests
430  trace.emphase() << ( res ? "Passed." : "Error." ) << endl;
431  trace.endBlock();
432  return res ? 0 : 1;
433 }
434 // //
void beginBlock(const std::string &keyword="")
Iterator read(Polynomial &p, Iterator begin, Iterator end)
graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_descriptor target(typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::edge_descriptor edge, const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
Go to http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/EdgeListGraph.html.
Definition: Boost.dox:171
MyDigitalSurface::SurfelSet SurfelSet
Go to http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/IncidenceGraph.html.
Definition: Boost.dox:168
Integer sKCoord(const SCell &c, Dimension k) const
Return its Khalimsky coordinate along [k].
Trace trace
Definition: Common.h:144
Component Coordinate
Type for Point elements.
Definition: PointVector.h:617
Aim: A utility class for constructing surfaces (i.e. set of (n-1)-cells).
Definition: Surfaces.h:78
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::edge_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::edge_iterator > edges(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
void attach(ConstAlias< EuclideanShape > shape)
Go to http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/AdjacencyGraph.html.
Definition: Boost.dox:165
Aim: This class converts a string polynomial expression in a multivariate polynomial.
const Point & upperBound() const
double endBlock()
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Domain domain
PointVector< 3, double > RealPoint
graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertices_size_type num_vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
bool testDigitalSurfaceBoostGraphInterface()
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::out_edge_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::out_edge_iterator > out_edges(typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_descriptor u, const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
std::ostream & emphase()
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
Aim: model of CEuclideanOrientedShape concepts to create a shape from a polynomial.
int main(int, char **)
DGtal is the top-level namespace which contains all DGtal functions and types.
boost::property< boost::vertex_index_t, std::size_t, boost::property< surfel_position_t, surfel_position > > VertexProperties
graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_descriptor source(typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::edge_descriptor edge, const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
std::ostream & info()
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
Go to http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/VertexListGraph.html.
Definition: Boost.dox:162
boost::uint64_t uint64_t
unsigned 64-bit integer.
Definition: BasicTypes.h:65
const Point & lowerBound() const
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
Definition: SetOfSurfels.h:73
KSpace K
Domain getDomain() const
std::ostream & error()
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface