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 extern "C" {
   9 #include "bgl.h"
  10 #include "mgraph.h"
  11 #include "common.h"
  12 }
  13
  14 /*
  15  * edge descriptor f =   00 | e | y     | x
  16  *                            3  YBITS   XBITS
  17  *
  18  * e is 0..5.  The edge is edge e out of vertex (x,y).
  19  *
  20  * BGL expects an undirected graph's edges to have two descriptors
  21  * each, one in each direction.
  22  */
  23
  24 /*
  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).
  34  */
  35
  36 #define VMASK (YMASK|XMASK)
  37 #define ESHIFT (YBITS|XBITS)
  38
  39 class Graph { }; // this is a dummy as our graph has no actual representation
  40
  41 namespace boost {
  42   // We make Graph a model of various BGL Graph concepts.
  43   // This mainly means that graph_traits<Graph> has lots of stuff.
  44
  45   // First, some definitions used later:
  46
  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 { };
  51
  52   struct OutEdgeIncrable {
  53     int f;
  54     OutEdgeIncrable& operator++() { f += 1<<ESHIFT; return self; }
  55     OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
  56   };
  57
  58   struct graph_traits<Graph> {
  59
  60     // Concept Graph:
  61
  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; }
  68
  69     // Concept IncidenceGraph:
  70
  71     typedef counting_iterator<OutEdgeIncrable,
  72       forward_iterator_tag> out_edge_iterator;
  73     typedef int degree_size_type;
  74
  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))));
  81     }
  82     inline out_degree(int v, const Graph&) { return VE_MAX(v) - VE_MIN(v); }
  83
  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));
  90     }
  91     inline unsigned num_vertices(const Graph&) { return N; }
  92   };
  93 };
  94
  95 static void single_source_shortest_paths(int v1,
  96                                          const double edge_weights[/*f*/],
  97                                          double vertex_distances[/*v*/]) {
  98   Graph g;
  99
 100   boost::dijkstra_shortest_paths(g, v1,
 101      weight_map(edge_weights).
 102      vertex_index_map(identity_property_map()).
 103      distance_map(vertex_distances));
 104 }
 105
 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|.
 109    *
 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.
 113    *
 114    * Energy contribution is proportional to
 115    *
 116    *               -4          2
 117    *    a  a   .  d   . [ (s/d)  - 1 ]
 118    *     vi vj
 119    *
 120    * (In practice we compute d^2+epsilon and use it for the
 121    *  divisions, to avoid division by zero.)
 122    */
 123   static const d2_epsilon= 1e-6;
 124
 125   double edge_weights[N*V6], vertex_distances[N], total_cost;
 126   int v1,v2,e,f;
 127
 128   FOR_VEDGE_X(v1,e,v2,
 129               f= v1 | e << ESHIFT,
 130               edge_weights[f]= NaN)
 131     edge_weights[f]= hypotD(v[v1], v[v2]);
 132
 133   FOR_VERTEX(v1) {
 134     double a1= vertex_areas[v1];
 135     single_source_shortest_paths(v1, edge_weights, vertex_distances);
 136     FOR_VERTEX(v2) {
 137       double a2= vertex_areas[v2];
 138       double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
 139       double sd= vertex_distances[v2] / d2;
 140       double sd2= sd*sd;
 141       total_cost += a1*a2 * (sd2 - 1) / (d2*d2);
 142     }
 143   }
 144   return total_cost;
 145 }