2 * Everything that needs the Boost Graph Library and C++ templates etc.
3 * (and what a crazy set of stuff that all is)
10 #include <boost/config.hpp>
11 #include <boost/iterator/iterator_facade.hpp>
12 #include <boost/graph/graph_traits.hpp>
13 #include <boost/graph/graph_concepts.hpp>
14 #include <boost/graph/dijkstra_shortest_paths.hpp>
15 #include <boost/graph/properties.hpp>
16 #include <boost/iterator/counting_iterator.hpp>
17 #include <boost/iterator/iterator_categories.hpp>
25 * edge descriptor f = 00 | e | y | x
28 * e is 0..5. The edge is edge e out of vertex (x,y).
30 * BGL expects an undirected graph's edges to have two descriptors
31 * each, one in each direction.
35 * We use BGL's implementation of Dijkstra's single source shortest
36 * paths. We really want all pairs shortest paths, so Johnson All
37 * Pairs Shortest Paths would seem sensible. But actually Johnson's
38 * algorithm is just a wrapper around Dijkstra's; the extra
39 * functionality is just to deal with -ve edge weights, which we don't
40 * have. So we can use Dijkstra directly and save some cpu (and some
41 * code: we don't have to supply all of the machinery needed for
42 * Johnson's invocation of Bellman-Ford). The overall time cost is
43 * O(VE log V); I think the space used is O(E).
46 #define VMASK (YMASK|XMASK)
47 #define ESHIFT (YBITS|XBITS)
49 class Graph { }; // this is a dummy as our graph has no actual representation
51 using namespace boost;
53 struct OutEdgeIterator :
54 public iterator_facade<
60 void increment() { f += 1<<ESHIFT; }
61 bool equal(OutEdgeIterator const& other) const { return f == other.f; }
62 int const& dereference() const { return f; }
64 OutEdgeIterator(int _f) : f(_f) { }
65 OutEdgeIterator(int v, int e) : f(e << ESHIFT | v) { }
68 typedef counting_iterator<int> VertexIterator;
71 // We make Graph a model of various BGL Graph concepts.
72 // This mainly means that graph_traits<Graph> has lots of stuff.
74 // First, some definitions used later:
76 struct layout_graph_traversal_category :
77 public virtual incidence_graph_tag,
78 public virtual vertex_list_graph_tag,
79 public virtual edge_list_graph_tag { };
81 struct graph_traits<Graph> {
83 typedef int vertex_descriptor; /* vertex number, -1 => none */
84 typedef int edge_descriptor; /* see above */
85 typedef undirected_tag directed_category;
86 typedef disallow_parallel_edge_tag edge_parallel_category;
87 typedef layout_graph_traversal_category traversal_category;
89 // Concept IncidenceGraph:
90 typedef OutEdgeIterator out_edge_iterator;
91 typedef unsigned degree_size_type;
93 // Concept VertexListGraph:
94 typedef VertexIterator vertex_iterator;
95 typedef unsigned vertices_size_type;
99 inline int null_vertex() { return -1; }
101 // Concept IncidenceGraph:
102 inline int source(int f, const Graph&) { return f&VMASK; }
103 inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
104 inline std::pair<OutEdgeIterator,OutEdgeIterator>
105 out_edges(int v, const Graph&) {
106 return std::make_pair(OutEdgeIterator(v, VE_MIN(v)),
107 OutEdgeIterator(v, VE_MAX(v)));
109 inline unsigned out_degree(int v, const Graph&) {
110 return VE_MAX(v) - VE_MIN(v);
113 // Concept VertexListGraph:
114 inline std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
115 return std::make_pair(VertexIterator(0), VertexIterator(N));
117 inline unsigned num_vertices(const Graph&) { return N; }
120 static void single_source_shortest_paths(int v1,
121 const double edge_weights[/*f*/],
122 double vertex_distances[/*v*/]) {
125 dijkstra_shortest_paths(g, v1,
126 weight_map(edge_weights).
127 vertex_index_map(identity_property_map()).
128 distance_map(vertex_distances));
131 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
132 /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
133 * along edges, and actual distance d_ij = |vi-vj|.
135 * We will also use the `vertex areas': for each vertex vi the
136 * vertex area a_vi is the mean area of the incident triangles.
137 * This is computed elsewhere.
139 * Energy contribution is proportional to
142 * a a . d . [ (s/d) - 1 ]
145 * (In practice we compute d^2+epsilon and use it for the
146 * divisions, to avoid division by zero.)
148 static const double d2_epsilon= 1e-6;
150 double edge_weights[N*V6], vertex_distances[N], total_cost=0;
156 edge_weights[f]= NAN)
157 edge_weights[f]= hypotD(v[v1], v[v2]);
160 double a1= vertex_areas[v1];
161 single_source_shortest_paths(v1, edge_weights, vertex_distances);
163 double a2= vertex_areas[v2];
164 double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
165 double sd= vertex_distances[v2] / d2;
167 total_cost += a1*a2 * (sd2 - 1) / (d2*d2);