2 * Everything that needs the Boost Graph Library and C++ templates etc.
3 * (and what a crazy set of stuff that all is)
15 * edge descriptor f = 00 | e | y | x
18 * e is 0..5. The edge is edge e out of vertex (x,y).
20 * BGL expects an undirected graph's edges to have two descriptors
21 * each, one in each direction.
25 * We use BGL's implementation of Dijkstra's single source shortest
26 * paths. We really want all pairs shortest paths, so Johnson All
27 * Pairs Shortest Paths would seem sensible. But actually Johnson's
28 * algorithm is just a wrapper around Dijkstra's; the extra
29 * functionality is just to deal with -ve edge weights, which we don't
30 * have. So we can use Dijkstra directly and save some cpu (and some
31 * code: we don't have to supply all of the machinery needed for
32 * Johnson's invocation of Bellman-Ford). The overall time cost is
33 * O(VE log V); I think the space used is O(E).
36 #define VMASK (YMASK|XMASK)
37 #define ESHIFT (YBITS|XBITS)
39 class Graph { }; // this is a dummy as our graph has no actual representation
42 // We make Graph a model of various BGL Graph concepts.
43 // This mainly means that graph_traits<Graph> has lots of stuff.
45 // First, some definitions used later:
47 struct layout_graph_traversal_category :
48 public virtual incidence_graph_tag,
49 public virtual vertex_list_graph_tag,
50 public virtual edge_list_graph_tag { };
52 struct OutEdgeIncrable {
54 OutEdgeIncrable& operator++() { f += 1<<ESHIFT; return self; }
55 OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
58 struct graph_traits<Graph> {
62 typedef int vertex_descriptor; /* vertex number, -1 => none */
63 typedef int edge_descriptor; /* see above */
64 typedef undirected_tag directed_category;
65 typedef disallow_parallel_ege edge_parallel_category;
66 typedef layout_graph_traversal_category traversal_category;
67 inline int null_vertex() { return -1; }
69 // Concept IncidenceGraph:
71 typedef counting_iterator<OutEdgeIncrable,
72 forward_iterator_tag> out_edge_iterator;
73 typedef int degree_size_type;
75 inline int source(int f, const Graph&) { return f&VMASK; }
76 inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
77 inline std::pair<out_edge_iterator,out_edge_iterator>
78 out_edges(int v, const Graph&) {
79 return std::make_pair(out_edge_iterator(OutEdgeIncrable(v, VE_MIN(v))),
80 out_edge_iterator(OutEdgeIncrable(v, VE_MAX(v))));
82 inline out_degree(int v, const Graph&) { return VE_MAX(v) - VE_MIN(v); }
84 // Concept VertexListGraph:
85 typedef counting_iterator<int> vertex_iterator;
86 typedef unsigned vertices_size_type;
87 inline std::pair<vertex_iterator,vertex_iterator>
88 vertices(const Graph&) {
89 return std::make_pair(vertex_iterator(0), vertex_iterator(N));
91 inline unsigned num_vertices(const Graph&) { return N; }
95 static void single_source_shortest_paths(int v1,
96 const double edge_weights[/*f*/],
97 double vertex_distances[/*v*/]) {
100 boost::dijkstra_shortest_paths(g, v1,
101 weight_map(edge_weights).
102 vertex_index_map(identity_property_map()).
103 distance_map(vertex_distances));
106 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
107 /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
108 * along edges, and actual distance d_ij = |vi-vj|.
110 * We will also use the `vertex areas': for each vertex vi the
111 * vertex area a_vi is the mean area of the incident triangles.
112 * This is computed elsewhere.
114 * Energy contribution is proportional to
117 * a a . d . [ (s/d) - 1 ]
120 * (In practice we compute d^2+epsilon and use it for the
121 * divisions, to avoid division by zero.)
123 static const d2_epsilon= 1e-6;
125 double edge_weights[N*V6], vertex_distances[N], total_cost;
130 edge_weights[f]= NaN)
131 edge_weights[f]= hypotD(v[v1], v[v2]);
134 double a1= vertex_areas[v1];
135 single_source_shortest_paths(v1, edge_weights, vertex_distances);
137 double a2= vertex_areas[v2];
138 double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
139 double sd= vertex_distances[v2] / d2;
141 total_cost += a1*a2 * (sd2 - 1) / (d2*d2);