chiark / gitweb /
325617eca59b114fff476833f347eeb3cc718e90
1 /*
2  * Everything that needs the Boost Graph Library and C++ templates etc.
3  * (and what a crazy set of stuff that all is)
4  */
6 #include <math.h>
8 #include <iterator>
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>
19 extern "C" {
20 #include "bgl.h"
21 #include "mgraph.h"
22 }
24 /*
25  * edge descriptor f =   0000 | e | y     | x
26  *                              3  YBITS   XBITS
27  *
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.
30  *
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).
33  */
35 /*
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).
45  */
47 #define VMASK (YMASK|XMASK)
48 #define ESHIFT (YBITS+XBITS)
50 class Graph { }; // this is a dummy as our graph has no actual representation
52 using namespace boost;
54 /*
55  * We iterate over edges in the following order:
56  *
57  *          \#0  /1#
58  *           \  /
59  *        ___ 0   __
60  *       #2    1   #3
61  *           /  \
62  *        #4/  #5\     and finally #6 is V6
63  *
64  *
65  * This ordering permits the order-4 nodes at the strip's edge
66  * to have a contiguous edge iterator values.  The iterator
67  * starts at #0 which is edge 2 (see mgraph.h), or #2 (edge 3).
68  */
69 static const int oei_edge_delta[V6]=
70   /*     0       1       2       3       4       5      initial e
71    *    #3      #1      #0      #2      #4      #5      initial ix
72    *    #4      #2      #1      #3      #5      #6      next ix
73    *     4       3       1       0       5      V6      next e
74    */ {
75      4<<ESHIFT, 2<<ESHIFT, -1<<ESHIFT,
76                               -3<<ESHIFT, 1<<ESHIFT, (V6-5)<<ESHIFT
77 };
79 class OutEdgeIterator :
80   public iterator_facade<
81     OutEdgeIterator,
82     int const,
83     forward_traversal_tag
84 > {
85   int f;
86  public:
87   void increment() {
88     //printf("incrementing f=%03x..",f);
89     f += oei_edge_delta[f>>ESHIFT];
90     //printf("%03x\n",f);
91   }
92   bool equal(OutEdgeIterator const& other) const { return f == other.f; }
93   int const& dereference() const { return f; }
94   OutEdgeIterator() { }
95   OutEdgeIterator(int _f) : f(_f) { }
96   OutEdgeIterator(int v, int e) : f(e<<ESHIFT | v) {
97     //printf("constructed v=%02x e=%x f=%03x\n",v,e,f);
98   }
100   static int voe_min(int _v) { return (_v & YMASK) ? 2 : 3; }
101   static int voe_max(int _v) { return (_v & YMASK)==(Y-1) ? V6 : 4; }
102   static int voe_degree(int _v) { return RIM_VERTEX_P(_v) ? 4 : V6; }
103 };
105 typedef counting_iterator<int> VertexIterator;
107 namespace boost {
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 { };
118   struct graph_traits<Graph> {
119     // Concept Graph:
120     typedef int vertex_descriptor; /* vertex number, -1 => none */
121     typedef int edge_descriptor; /* see above */
122     typedef undirected_tag directed_category;
123     typedef disallow_parallel_edge_tag edge_parallel_category;
124     typedef layout_graph_traversal_category traversal_category;
126     // Concept IncidenceGraph:
127     typedef OutEdgeIterator out_edge_iterator;
128     typedef unsigned degree_size_type;
130     // Concept VertexListGraph:
131     typedef VertexIterator vertex_iterator;
132     typedef unsigned vertices_size_type;
133   };
135   // Concept Graph:
136   inline int null_vertex() { return -1; }
138   // Concept IncidenceGraph:
139   inline int source(int f, const Graph&) { return f&VMASK; }
140   inline int target(int f, const Graph&) {
141     int v2= EDGE_END2(f&VMASK, f>>ESHIFT);
142     //printf("traversed %03x..%02x\n",f,v2);
143     return v2;
144   }
145   inline std::pair<OutEdgeIterator,OutEdgeIterator>
146   out_edges(int v, const Graph&) {
147     return std::make_pair(OutEdgeIterator(v, OutEdgeIterator::voe_min(v)),
148                           OutEdgeIterator(v, OutEdgeIterator::voe_max(v)));
149   }
150   inline unsigned out_degree(int v, const Graph&) {
151     return OutEdgeIterator::voe_degree(v);
152   }
154   // Concept VertexListGraph:
155   inline std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
156     return std::make_pair(VertexIterator(0), VertexIterator(N));
157   }
158   inline unsigned num_vertices(const Graph&) { return N; }
159 };
161 static void single_source_shortest_paths(int v1,
162                                          const double edge_weights[/*f*/],
163                                          double vertex_distances[/*v*/]) {
164   Graph g;
166   dijkstra_shortest_paths(g, v1,
167      weight_map(edge_weights).
168      vertex_index_map(identity_property_map()).
169      distance_map(vertex_distances));
170 }
172 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
173   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
174    * along edges, and actual distance d_ij = |vi-vj|.
175    *
176    * We will also use the `vertex areas': for each vertex vi the
177    * vertex area a_vi is the mean area of the incident triangles.
178    * This is computed elsewhere.
179    *
180    * Energy contribution is proportional to
181    *
182    *               -4          2
183    *    a  a   .  d   . [ (s/d)  - 1 ]
184    *     vi vj
185    *
186    * (In practice we compute d^2+epsilon and use it for the
187    *  divisions, to avoid division by zero.)
188    */
189   static const double d2_epsilon= 1e-6;
191   double edge_weights[N*V6], vertex_distances[N], total_cost=0;
192   int v1,v2,e,f;
194   FOR_VERTEX(v1)
195     FOR_VEDGE_X(v1,e,v2,
196                 f= v1 | e << ESHIFT,
197                 edge_weights[f]= NAN)
198       edge_weights[f]= hypotD(v[v1], v[v2]);
200   FOR_VERTEX(v1) {
201     double a1= vertex_areas[v1];
202     single_source_shortest_paths(v1, edge_weights, vertex_distances);
203     FOR_VERTEX(v2) {
204       if (v1 == v2) continue;
205       double a2= vertex_areas[v2];
206       double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
207       double s= vertex_distances[v2];
208       double s2= s*s + d2_epsilon;
209       double sd2= s2 / d2;
210       double cost_contrib= a1*a2 * (sd2 - 1) / (d2*d2);
211       if (cost_contrib < -1e-4) {
212         printf("layout %03x..%03x (a=%g,%g) s=%g s2=%g d2=%g sd2=%g"
213                " cost+=%g\n", v1,v2, a1,a2, s,s2,d2,sd2, cost_contrib);
214         abort();
215       }
216       total_cost += cost_contrib;
217     }
218   }
219   return total_cost;
220 }