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>
64struct surfel_position_t {
65 typedef boost::vertex_property_tag kind;
68struct surfel_position {
70 surfel_position() : myP()
75typedef boost::property< boost::vertex_index_t, std::size_t,
79template <
typename Graph1,
typename Graph2,
typename VertexIndexMap>
80struct 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;
87 graph_vertex_index_map graph_vertex_index;
88 graph_vertex_position_map graph_vertex_position;
89 VertexIndexMap & myIndexMap;
91 my_vertex_copier(
const Graph1& g1, Graph2& g2, VertexIndexMap & indexMap )
93 graph_vertex_index( boost::get( boost::vertex_index_t(), g2 ) ),
94 graph_vertex_position( boost::get( surfel_position_t(), g2) ),
95 myIndexMap( indexMap )
98 void operator()(
const Vertex1& v1,
const Vertex2& v2 )
const {
103 pos.myP = v1.myCoordinates;
105 put( graph_vertex_position, v2, pos);
108template <
typename Graph1,
typename Graph2>
109struct my_edge_copier {
110 my_edge_copier(
const Graph1& UNUSED(g1), Graph2& UNUSED(g2))
112 template <
typename Edge1,
typename Edge2>
113 void operator()(
const Edge1& , Edge2& )
const {
124 unsigned int nbok = 0;
140 double p1[] = {-2,-2,-2};
141 double p2[] = { 2, 2, 2};
142 std::string poly_str =
"x*x+y*y+z*z-1";
146 std::string::const_iterator iter
147 = reader.
read( P, poly_str.begin(), poly_str.end() );
148 if ( iter != poly_str.end() )
150 std::cerr <<
"ERROR: I read only <"
151 << poly_str.substr( 0, iter - poly_str.begin() )
152 <<
">, and I built P=" << P << std::endl;
155 trace.
info() <<
"- P = " << P << std::endl;
164 trace.
beginBlock (
"Construct the Khalimsky space from the image domain ..." );
171 trace.
error() <<
"Error in the Khamisky space construction."<<std::endl;
184 MySurfelAdjacency surfAdj(
true );
185 MySetOfSurfels theSetOfSurfels(
K, surfAdj );
191 trace.
info() <<
"Digital surface has " << digSurf.
size() <<
" surfels."
197 typedef boost::graph_traits<Graph>::vertex_descriptor vertex_descriptor;
198 typedef boost::graph_traits<Graph>::edge_descriptor edge_descriptor;
199 typedef boost::graph_traits<Graph>::vertices_size_type vertices_size_type;
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;
211 trace.
beginBlock (
"Testing IncidenceGraph interface with breadth_first_visit ..." );
213 typedef std::map< vertex_descriptor, boost::default_color_type > StdColorMap;
214 StdColorMap colorMap;
215 boost::associative_property_map< StdColorMap > propColorMap( colorMap );
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;
221 vertex_descriptor start = *( digSurf.
begin() );
222 boost::breadth_first_visit
226 boost::make_bfs_visitor( boost::record_distances( propDistanceMap, boost::on_tree_edge() ) ),
234 vertex_descriptor furthest = start;
236 for ( std::pair<vertex_iterator, vertex_iterator>
237 vp =
boost::vertices( digSurf ); vp.first != vp.second; ++vp.first, ++nbV )
239 uint64_t d = boost::get( propDistanceMap, *vp.first );
243 furthest = *vp.first;
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;
256 trace.
beginBlock (
"Testing VertexListGraph interface with connected_components ..." );
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
265 boost::color_map( propColorMap )
267 trace.
info() <<
"- nbComponents = " << nbComp << std::endl;
268 ++nb; nbok += ( nbComp == 1 ) ? 1 : 0;
269 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") "
270 <<
"nbComp == 1" << std::endl;
273 trace.
beginBlock (
"Testing UndirectedGraph interface with Wagner-Stoer min cut ..." );
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;
285 for ( std::pair<vertex_iterator, vertex_iterator>
286 vp =
boost::vertices( digSurf ); vp.first != vp.second; ++vp.first, ++idxV )
288 vertex_descriptor v1 = *vp.first;
289 vertexIndexMap[ v1 ] = idxV;
290 for ( std::pair<out_edge_iterator, out_edge_iterator>
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;
304 typedef std::map< vertex_descriptor, bool > StdParityMap;
305 StdParityMap parityMap;
306 boost::associative_property_map< StdParityMap > propParityMap( parityMap );
308 weight_type total_weight =
309 boost::stoer_wagner_min_cut
312 boost::parity_map( propParityMap )
313 .vertex_index_map( propVertexIndexMap )
315 trace.
info() <<
"- total weight = " << total_weight << std::endl;
318 for ( std::pair<vertex_iterator, vertex_iterator>
319 vp =
boost::vertices( digSurf ); vp.first != vp.second; ++vp.first, ++idxV )
321 vertex_descriptor v1 = *vp.first;
323 if ( parityMap[ v1 ] ) ++nb1;
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;
337 trace.
beginBlock (
"Testing EdgeListGraph and IncidenceGraph interfaces with Boykov-Kolmogorov max flow ..." );
338 typedef double capacity_type;
340 typedef std::map< edge_descriptor, weight_type > StdEdgeCapacityMap;
341 StdEdgeCapacityMap edgeCapacityMap;
342 boost::associative_property_map< StdEdgeCapacityMap > propEdgeCapacityMap( edgeCapacityMap );
344 typedef std::map< edge_descriptor, weight_type > StdEdgeResidualCapacityMap;
345 StdEdgeResidualCapacityMap edgeResidualCapacityMap;
346 boost::associative_property_map< StdEdgeResidualCapacityMap > propEdgeResidualCapacityMap( edgeResidualCapacityMap );
348 typedef std::map< edge_descriptor, edge_descriptor > StdReversedEdgeMap;
349 StdReversedEdgeMap reversedEdgeMap;
350 boost::associative_property_map< StdReversedEdgeMap > propReversedEdgeMap( reversedEdgeMap );
352 typedef std::map< vertex_descriptor, edge_descriptor > StdPredecessorMap;
353 StdPredecessorMap predecessorMap;
354 boost::associative_property_map< StdPredecessorMap > propPredecessorMap( predecessorMap );
359 for ( std::pair<edge_iterator, edge_iterator>
360 ve =
boost::edges( digSurf ); ve.first != ve.second; ++ve.first, ++nbEdges )
362 edge_descriptor e = *ve.first;
363 edge_descriptor rev_e = digSurf.
opposite( e );
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;
378 trace.
info() <<
"- nb edges = " << nbEdges << std::endl;
381 capacity_type max_flow =
382 boost::boykov_kolmogorov_max_flow
384 propEdgeCapacityMap, propEdgeResidualCapacityMap,
385 propReversedEdgeMap, propPredecessorMap, propColorMap, propDistanceMap, propVertexIndexMap,
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;
394 typedef boost::adjacency_list< boost::vecS, boost::vecS, boost::undirectedS, VertexProperties > BGraph;
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 )
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>
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;
413 trace.
info() <<
"(" << nbok <<
"/" << nb <<
") "
414 <<
"after copy: Boost graph has " << num_vertices( bG )
415 <<
" vertices." << std::endl;
427 trace.
beginBlock (
"Testing class DigitalSurfaceBoostGraphInterface" );
430 trace.
emphase() << ( res ?
"Passed." :
"Error." ) << endl;
Aim: Represents a set of n-1-cells in a nD space, together with adjacency relation between these cell...
Arc opposite(const Arc &a) const
ConstIterator begin() const
SCell separator(const Arc &a) const
Aim: A class for computing the Gauss digitization of some Euclidean shape, i.e. its intersection with...
void attach(ConstAlias< EuclideanShape > shape)
void init(const RealPoint &xLow, const RealPoint &xUp, typename RealVector::Component gridStep)
const Point & lowerBound() const
const Point & upperBound() const
Aim: model of CEuclideanOrientedShape concepts to create a shape from a polynomial.
Aim: This class is a model of CCellularGridSpaceND. It represents the cubical grid as a cell complex,...
std::set< SCell > SurfelSet
Preferred type for defining a set of surfels (always signed cells).
bool init(const Point &lower, const Point &upper, bool isClosed)
Specifies the upper and lower bounds for the maximal cells in this space.
Integer sKCoord(const SCell &c, Dimension k) const
Return its Khalimsky coordinate along [k].
Aim: This class converts a string polynomial expression in a multivariate polynomial.
Iterator read(Polynomial &p, Iterator begin, Iterator end)
Aim: Represents a multivariate polynomial, i.e. an element of , where K is some ring or field.
Component Coordinate
Type for Point elements.
Aim: A model of CDigitalSurfaceContainer which defines the digital surface as connected surfels....
static void trackBoundary(SCellSet &surface, const KSpace &K, const SurfelAdjacency< KSpace::dimension > &surfel_adj, const PointPredicate &pp, const SCell &start_surfel)
static SCell findABel(const KSpace &K, const PointPredicate &pp, unsigned int nbtries=1000)
Aim: Represent adjacencies between surfel elements, telling if it follows an interior to exterior ord...
void beginBlock(const std::string &keyword="")
DigitalSurface< MyDigitalSurfaceContainer > MyDigitalSurface
MyDigitalSurface::SurfelSet SurfelSet
DGtal is the top-level namespace which contains all DGtal functions and types.
boost::uint64_t uint64_t
unsigned 64-bit integer.
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::edge_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::edge_iterator > edges(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertices_size_type num_vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_descriptor source(typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::edge_descriptor edge, const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
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)
std::pair< typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator, typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_iterator > vertices(const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::vertex_descriptor target(typename graph_traits< DGtal::DigitalSurface< TDigitalSurfaceContainer > >::edge_descriptor edge, const DGtal::DigitalSurface< TDigitalSurfaceContainer > &digSurf)
Represents a signed cell in a cellular grid space by its Khalimsky coordinates and a boolean value.
Go to http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/AdjacencyGraph.html.
Go to http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/EdgeListGraph.html.
Go to http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/IncidenceGraph.html.
Go to http://www.boost.org/doc/libs/1_52_0/libs/graph/doc/VertexListGraph.html.
bool testDigitalSurfaceBoostGraphInterface()
boost::property< boost::vertex_index_t, std::size_t, boost::property< surfel_position_t, surfel_position > > VertexProperties