chiark / gitweb /
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 =   00 | e | y     | x
26  *                            3  YBITS   XBITS
27  *
28  * e is 0..5.  The edge is edge e out of vertex (x,y).
29  *
30  * BGL expects an undirected graph's edges to have two descriptors
31  * each, one in each direction (otherwise e would be just 0..2).
32  */
34 /*
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).
44  */
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 /*
54  * We use the following alternative numbering for iterating over edges:
55  *
56  *      ix:         \0   /1
57  *                   \  /
58  *                ___ 0   __
59  *                2    1   3
60  *                   /  \
61  *                 4/   5\
62  *
63  *
64  * This numbering permits the order-4 nodes at the strip's edge
65  * to have a contiguous edge iterator values 2..5 or 0..3.
66  */
67 static const int oee_edgemap[V6]=
68   { 2<<ESHIFT, 1<<ESHIFT, 3<<ESHIFT, 2<<ESHIFT, 4<<ESHIFT, 5<<ESHIFT };
70 class OutEdgeIterator :
71   public iterator_facade<
72     OutEdgeIterator,
73     int const,
74     forward_traversal_tag
75 > {
76   int v, ix, f;
78   void setf() { f= v | oee_edgemap[ix]; }
79  public:
80   void increment() { ix++; setf(); }
81   bool equal(OutEdgeIterator const& other) const { return ix == other.ix; }
82   int const& dereference() const { return f; }
83   OutEdgeIterator() { }
84   OutEdgeIterator(int _v, int _ix) : v(_v), ix(_ix) { setf(); }
86   static int voe_min(int _v) { return _v & YMASK ? 0 : 2; }
87   static int voe_max(int _v) { return ~_v & YMASK ? V6 : 4; }
88 };
90 typedef counting_iterator<int> VertexIterator;
92 namespace boost {
93   // We make Graph a model of various BGL Graph concepts.
94   // This mainly means that graph_traits<Graph> has lots of stuff.
96   // First, some definitions used later:
98   struct layout_graph_traversal_category :
99     public virtual incidence_graph_tag,
100     public virtual vertex_list_graph_tag,
101     public virtual edge_list_graph_tag { };
103   struct graph_traits<Graph> {
104     // Concept Graph:
105     typedef int vertex_descriptor; /* vertex number, -1 => none */
106     typedef int edge_descriptor; /* see above */
107     typedef undirected_tag directed_category;
108     typedef disallow_parallel_edge_tag edge_parallel_category;
109     typedef layout_graph_traversal_category traversal_category;
111     // Concept IncidenceGraph:
112     typedef OutEdgeIterator out_edge_iterator;
113     typedef unsigned degree_size_type;
115     // Concept VertexListGraph:
116     typedef VertexIterator vertex_iterator;
117     typedef unsigned vertices_size_type;
118   };
120   // Concept Graph:
121   inline int null_vertex() { return -1; }
123   // Concept IncidenceGraph:
124   inline int source(int f, const Graph&) { return f&VMASK; }
125   inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
126   inline std::pair<OutEdgeIterator,OutEdgeIterator>
127   out_edges(int v, const Graph&) {
128     return std::make_pair(OutEdgeIterator(v, OutEdgeIterator::voe_min(v)),
129                           OutEdgeIterator(v, OutEdgeIterator::voe_max(v)));
130   }
131   inline unsigned out_degree(int v, const Graph&) {
132     return OutEdgeIterator::voe_max(v) - OutEdgeIterator::voe_min(v);
133   }
135   // Concept VertexListGraph:
136   inline std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
137     return std::make_pair(VertexIterator(0), VertexIterator(N));
138   }
139   inline unsigned num_vertices(const Graph&) { return N; }
140 };
142 static void single_source_shortest_paths(int v1,
143                                          const double edge_weights[/*f*/],
144                                          double vertex_distances[/*v*/]) {
145   Graph g;
147   dijkstra_shortest_paths(g, v1,
148      weight_map(edge_weights).
149      vertex_index_map(identity_property_map()).
150      distance_map(vertex_distances));
151 }
153 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
154   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
155    * along edges, and actual distance d_ij = |vi-vj|.
156    *
157    * We will also use the `vertex areas': for each vertex vi the
158    * vertex area a_vi is the mean area of the incident triangles.
159    * This is computed elsewhere.
160    *
161    * Energy contribution is proportional to
162    *
163    *               -4          2
164    *    a  a   .  d   . [ (s/d)  - 1 ]
165    *     vi vj
166    *
167    * (In practice we compute d^2+epsilon and use it for the
168    *  divisions, to avoid division by zero.)
169    */
170   static const double d2_epsilon= 1e-6;
172   double edge_weights[N*V6], vertex_distances[N], total_cost=0;
173   int v1,v2,e,f;
175   FOR_VERTEX(v1)
176     FOR_VEDGE_X(v1,e,v2,
177                 f= v1 | e << ESHIFT,
178                 edge_weights[f]= NAN)
179       edge_weights[f]= hypotD(v[v1], v[v2]);
181   FOR_VERTEX(v1) {
182     double a1= vertex_areas[v1];
183     single_source_shortest_paths(v1, edge_weights, vertex_distances);
184     FOR_VERTEX(v2) {
185       double a2= vertex_areas[v2];
186       double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
187       double sd= vertex_distances[v2] / d2;
188       double sd2= sd*sd;
189       total_cost += a1*a2 * (sd2 - 1) / (d2*d2);
190     }
191   }
192   return total_cost;
193 }