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 = 0000 | e | y | x
28 * e is 0..6. The edge is edge e out of vertex (x,y), or if
29 * e==6 it's the `at end' value for the out edge iterator.
31 * BGL expects an undirected graph's edges to have two descriptors
32 * each, one in each direction (otherwise e would be just 0..2).
36 * We use BGL's implementation of Dijkstra's single source shortest
37 * paths. We really want all pairs shortest paths, so Johnson All
38 * Pairs Shortest Paths would seem sensible. But actually Johnson's
39 * algorithm is just a wrapper around Dijkstra's; the extra
40 * functionality is just to deal with -ve edge weights, which we don't
41 * have. So we can use Dijkstra directly and save some cpu (and some
42 * code: we don't have to supply all of the machinery needed for
43 * Johnson's invocation of Bellman-Ford). The overall time cost is
44 * O(VE log V); I think the space used is O(E).
47 #define VMASK (YMASK|XMASK)
48 #define ESHIFT (YBITS+XBITS)
50 using namespace boost;
53 * We iterate over edges in the following order:
60 * #4/ #5\ and finally #6 is V6
63 * This ordering permits the order-4 nodes at the strip's edge
64 * to have a contiguous edge iterator values. The iterator
65 * starts at #0 which is edge 2 (see mgraph.h), or #2 (edge 3).
67 static const int oei_edge_delta[V6]=
68 /* 0 1 2 3 4 5 initial e
69 * #3 #1 #0 #2 #4 #5 initial ix
70 * #4 #2 #1 #3 #5 #6 next ix
73 4<<ESHIFT, 2<<ESHIFT, -1<<ESHIFT,
74 -3<<ESHIFT, 1<<ESHIFT, (V6-5)<<ESHIFT
77 class OutEdgeIterator :
78 public iterator_facade<
86 //printf("incrementing f=%03x..",f);
87 f += oei_edge_delta[f>>ESHIFT];
90 bool equal(OutEdgeIterator const& other) const { return f == other.f; }
91 int const& dereference() const { return f; }
93 OutEdgeIterator(int _f) : f(_f) { }
94 OutEdgeIterator(int v, int e) : f(e<<ESHIFT | v) {
95 //printf("constructed v=%02x e=%x f=%03x\n",v,e,f);
98 static int voe_min(int _v) { return (_v & YMASK) ? 2 : 3; }
99 static int voe_max(int _v) { return (_v & YMASK)==(Y-1) ? V6 : 4; }
100 static int voe_degree(int _v) { return RIM_VERTEX_P(_v) ? 4 : V6; }
103 typedef counting_iterator<int> VertexIterator;
106 class Graph { }; // this is a dummy as our graph has no actual representation
108 // We make Graph a model of various BGL Graph concepts.
109 // This mainly means that graph_traits<Graph> has lots of stuff.
111 // First, some definitions used later:
113 struct layout_graph_traversal_category :
114 public virtual incidence_graph_tag,
115 public virtual vertex_list_graph_tag,
116 public virtual edge_list_graph_tag { };
119 struct graph_traits<Graph> {
121 typedef int vertex_descriptor; /* vertex number, -1 => none */
122 typedef int edge_descriptor; /* see above */
123 typedef undirected_tag directed_category;
124 typedef disallow_parallel_edge_tag edge_parallel_category;
125 typedef layout_graph_traversal_category traversal_category;
127 // Concept IncidenceGraph:
128 typedef OutEdgeIterator out_edge_iterator;
129 typedef unsigned degree_size_type;
131 // Concept VertexListGraph:
132 typedef VertexIterator vertex_iterator;
133 typedef unsigned vertices_size_type;
137 inline int null_vertex() { return -1; }
139 // Concept IncidenceGraph:
140 inline int source(int f, const Graph&) { return f&VMASK; }
141 inline int target(int f, const Graph&) {
142 int v2= EDGE_END2(f&VMASK, f>>ESHIFT);
143 //printf("traversed %03x..%02x\n",f,v2);
146 inline std::pair<OutEdgeIterator,OutEdgeIterator>
147 out_edges(int v, const Graph&) {
148 return std::make_pair(OutEdgeIterator(v, OutEdgeIterator::voe_min(v)),
149 OutEdgeIterator(v, OutEdgeIterator::voe_max(v)));
151 inline unsigned out_degree(int v, const Graph&) {
152 return OutEdgeIterator::voe_degree(v);
155 // Concept VertexListGraph:
157 std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
158 return std::make_pair(VertexIterator(0), VertexIterator(N));
160 inline unsigned num_vertices(const Graph&) { return N; }
163 static void single_source_shortest_paths(int v1,
164 const double edge_weights[/*f*/],
165 double vertex_distances[/*v*/]) {
168 dijkstra_shortest_paths(g, v1,
169 weight_map(edge_weights).
170 vertex_index_map(identity_property_map()).
171 distance_map(vertex_distances));
174 static int distances[N][N];
176 void graph_layout_prepare() {
181 int *d= distances[v1];
182 FOR_VERTEX(v2) d[v2]= -1;
186 vertex_index_map(identity_property_map()).
187 visitor(make_bfs_visitor(record_distances(d,on_tree_edge()))));
189 FOR_VERTEX(v2) printf(" %02x:%d",v2,d[v2]);
194 int *d= distances[v1];
196 FOR_VERTEX(v2) printf(" %02x:%d",v2,d[v2]);
201 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
202 /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
203 * along edges, and actual distance d_ij = |vi-vj|.
205 * We will also use the `vertex areas': for each vertex vi the
206 * vertex area a_vi is the mean area of the incident triangles.
207 * This is computed elsewhere.
209 * Energy contribution is proportional to
212 * a a . d . [ (s/d) - 1 ]
215 * (In practice we compute d^2+epsilon and use it for the
216 * divisions, to avoid division by zero.)
218 static const double d2_epsilon= 1e-6;
220 // double edge_weights[V6<<ESHIFT], vertex_distances[N],
222 int v1,v2,e, nedges=0;
223 double totaledgelength=0, meanedgelength;
226 totaledgelength += hypotD(v[v1], v[v2]);
230 meanedgelength= totaledgelength / nedges;
234 if (v1 == v2) continue;
236 double d= hypotD(v[v1],v[v2]);
238 int dist= distances[v1][v2];
241 double s= dist * meanedgelength * 0.03;
243 double enoughdistance= d - s;
244 if (enoughdistance > 1e-6) continue;
246 /* energy = 1/2 stiffness deviation^2
247 * where stiffness = 1/d
250 double cost= pow(enoughdistance,4);
252 //double s2= s*s + d2_epsilon;
253 //double sd2= s2 / d2;
254 //double cost_contrib= a1*a2 * (sd2 - 1) / (d2*d2);
255 //double cost_contrib= sd2;
257 printf("layout %03x..%03x dist=%d mean=%g s=%g d=%g enough=%g"
258 " cost+=%g\n", v1,v2, dist, meanedgelength,
259 s,d, enoughdistance, cost);