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 =   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  */
  34
  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  */
  46
  47 #define VMASK (YMASK|XMASK)
  48 #define ESHIFT (YBITS+XBITS)
  49
  50 class Graph { }; // this is a dummy as our graph has no actual representation
  51
  52 using namespace boost;
  53
  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 };
  78
  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   }
  99
 100   static int voe_min(int _v) { return (_v & YMASK) ? 2 : 3; }
 101   static int voe_max(int _v) { return (~_v & YMASK) ? V6 : 4; }
 102   static int voe_degree(int _v) { return (_v & YMASK | ~_v & YMASK) ? 4 : V6; }
 103 };
 104
 105 typedef counting_iterator<int> VertexIterator;
 106
 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.
 110
 111   // First, some definitions used later:
 112
 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 { };
 117
 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;
 125
 126     // Concept IncidenceGraph:
 127     typedef OutEdgeIterator out_edge_iterator;
 128     typedef unsigned degree_size_type;
 129
 130     // Concept VertexListGraph:
 131     typedef VertexIterator vertex_iterator;
 132     typedef unsigned vertices_size_type;
 133   };
 134
 135   // Concept Graph:
 136   inline int null_vertex() { return -1; }
 137
 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   }
 153
 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 };
 160
 161 static void single_source_shortest_paths(int v1,
 162                                          const double edge_weights[/*f*/],
 163                                          double vertex_distances[/*v*/]) {
 164   Graph g;
 165
 166   dijkstra_shortest_paths(g, v1,
 167      weight_map(edge_weights).
 168      vertex_index_map(identity_property_map()).
 169      distance_map(vertex_distances));
 170 }
 171
 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;
 190
 191   double edge_weights[N*V6], vertex_distances[N], total_cost=0;
 192   int v1,v2,e,f;
 193
 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]);
 199
 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       double a2= vertex_areas[v2];
 205       double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
 206       double s= vertex_distances[v2];
 207       double sd= s / d2;
 208       double sd2= sd*sd;
 209       double cost_contrib= a1*a2 * (sd2 - 1) / (d2*d2);
 210       //printf("layout %03x..%03x (a=%g,%g) s=%g d2=%g cost+=%g\n",
 211       //             v1,v2, a1,a2, s,d2, cost_contrib);
 212       total_cost += cost_contrib;
 213     }
 214   }
 215   return total_cost;
 216 }