bgl.cpp

   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  */
   5
   6 #include <math.h>
   7
   8 #include <iterator>
   9
  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>
  18
  19 extern "C" {
  20 #include "bgl.h"
  21 #include "mgraph.h"
  22 }
  23
  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.
  32  */
  33
  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  */
  45
  46 #define VMASK (YMASK|XMASK)
  47 #define ESHIFT (YBITS|XBITS)
  48
  49 class Graph { }; // this is a dummy as our graph has no actual representation
  50
  51 struct OutEdgeIncrable {
  52   int f;
  53   OutEdgeIncrable& operator++() { f += 1<<ESHIFT; return *this; }
  54   OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
  55 };
  56
  57 namespace boost {
  58   // We make Graph a model of various BGL Graph concepts.
  59   // This mainly means that graph_traits<Graph> has lots of stuff.
  60
  61   // First, some definitions used later:
  62
  63   struct layout_graph_traversal_category :
  64     public virtual incidence_graph_tag,
  65     public virtual vertex_list_graph_tag,
  66     public virtual edge_list_graph_tag { };
  67
  68   struct graph_traits<Graph> {
  69
  70     // Concept Graph:
  71
  72     typedef int vertex_descriptor; /* vertex number, -1 => none */
  73     typedef int edge_descriptor; /* see above */
  74     typedef undirected_tag directed_category;
  75     typedef disallow_parallel_edge_tag edge_parallel_category;
  76     typedef layout_graph_traversal_category traversal_category;
  77     inline int null_vertex() { return -1; }
  78
  79     // Concept IncidenceGraph:
  80
  81     typedef counting_iterator<OutEdgeIncrable,
  82       std::forward_iterator_tag> out_edge_iterator;
  83     typedef unsigned degree_size_type;
  84
  85     inline int source(int f, const Graph&) { return f&VMASK; }
  86     inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
  87     inline std::pair<out_edge_iterator,out_edge_iterator>
  88     out_edges(int v, const Graph&) {
  89       return std::make_pair(out_edge_iterator(OutEdgeIncrable(v, VE_MIN(v))),
  90                             out_edge_iterator(OutEdgeIncrable(v, VE_MAX(v))));
  91     }
  92     inline unsigned out_degree(int v, const Graph&) {
  93       return VE_MAX(v) - VE_MIN(v);
  94     }
  95
  96     // Concept VertexListGraph:
  97     typedef counting_iterator<int> vertex_iterator;
  98     typedef unsigned vertices_size_type;
  99     inline std::pair<vertex_iterator,vertex_iterator>
 100     vertices(const Graph&) {
 101       return std::make_pair(vertex_iterator(0), vertex_iterator(N));
 102     }
 103     inline unsigned num_vertices(const Graph&) { return N; }
 104   };
 105 };
 106
 107 static void single_source_shortest_paths(int v1,
 108                                          const double edge_weights[/*f*/],
 109                                          double vertex_distances[/*v*/]) {
 110   Graph g;
 111
 112   boost::dijkstra_shortest_paths(g, v1,
 113      weight_map(edge_weights).
 114      vertex_index_map(identity_property_map()).
 115      distance_map(vertex_distances));
 116 }
 117
 118 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
 119   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
 120    * along edges, and actual distance d_ij = |vi-vj|.
 121    *
 122    * We will also use the `vertex areas': for each vertex vi the
 123    * vertex area a_vi is the mean area of the incident triangles.
 124    * This is computed elsewhere.
 125    *
 126    * Energy contribution is proportional to
 127    *
 128    *               -4          2
 129    *    a  a   .  d   . [ (s/d)  - 1 ]
 130    *     vi vj
 131    *
 132    * (In practice we compute d^2+epsilon and use it for the
 133    *  divisions, to avoid division by zero.)
 134    */
 135   static const d2_epsilon= 1e-6;
 136
 137   double edge_weights[N*V6], vertex_distances[N], total_cost;
 138   int v1,v2,e,f;
 139
 140   FOR_VEDGE_X(v1,e,v2,
 141               f= v1 | e << ESHIFT,
 142               edge_weights[f]= NaN)
 143     edge_weights[f]= hypotD(v[v1], v[v2]);
 144
 145   FOR_VERTEX(v1) {
 146     double a1= vertex_areas[v1];
 147     single_source_shortest_paths(v1, edge_weights, vertex_distances);
 148     FOR_VERTEX(v2) {
 149       double a2= vertex_areas[v2];
 150       double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
 151       double sd= vertex_distances[v2] / d2;
 152       double sd2= sd*sd;
 153       total_cost += a1*a2 * (sd2 - 1) / (d2*d2);
 154     }
 155   }
 156   return total_cost;
 157 }