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 (otherwise e would be just 0..2).
  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 using namespace boost;
  52
  53 /*
  54  * We iterate over edges in the following order:
  55  *
  56  *          \#0  /1#
  57  *           \  /
  58  *        ___ 0   __
  59  *       #2    1   #3
  60  *           /  \
  61  *        #4/  #5\     and finally #6 is V6
  62  *
  63  *
  64  * This ordering permits the order-4 nodes at the strip's edge
  65  * to have a contiguous edge iterator values.  The iterator
  66  * starts at #0 which is edge 2 (see mgraph.h), or #2 (edge 3).
  67  */
  68 static const int oei_edge_delta[V6]=
  69   /*     0       1       2       3       4       5      initial e
  70    *    #3      #1      #0      #2      #4      #5      initial ix
  71    *    #4      #2      #1      #3      #5      #6      next ix
  72    *     4       3       1       0       5      V6      next e
  73    */ {
  74      4<<ESHIFT, 2<<ESHIFT, -1<<ESHIFT,
  75                               -3<<ESHIFT, 1<<ESHIFT, (V6-5)<<ESHIFT
  76 };
  77
  78 class OutEdgeIterator :
  79   public iterator_facade<
  80     OutEdgeIterator,
  81     int const,
  82     forward_traversal_tag
  83 > {
  84   int f;
  85  public:
  86   void increment() {
  87     //printf("incrementing f=%03x..",f);
  88     f += oei_edge_delta[f>>ESHIFT];
  89     //printf("%03x\n",f);
  90   }
  91   bool equal(OutEdgeIterator const& other) const { return f == other.f; }
  92   int const& dereference() const { return f; }
  93   OutEdgeIterator() { }
  94   OutEdgeIterator(int _f) : f(_f) { }
  95   OutEdgeIterator(int v, int e) : f(e<<ESHIFT | v) {
  96     //printf("constructed v=%x e=%x f=%03x\n",v,e,f);
  97   }
  98
  99   static int voe_min(int _v) { return (_v & YMASK) ? 2 : 3; }
 100   static int voe_max(int _v) { return (~_v & YMASK) ? V6 : 4; }
 101   static int voe_degree(int _v) { return (_v & YMASK | ~_v & YMASK) ? 4 : V6; }
 102 };
 103
 104 typedef counting_iterator<int> VertexIterator;
 105
 106 namespace boost {
 107   // We make Graph a model of various BGL Graph concepts.
 108   // This mainly means that graph_traits<Graph> has lots of stuff.
 109
 110   // First, some definitions used later:
 111
 112   struct layout_graph_traversal_category :
 113     public virtual incidence_graph_tag,
 114     public virtual vertex_list_graph_tag,
 115     public virtual edge_list_graph_tag { };
 116
 117   struct graph_traits<Graph> {
 118     // Concept Graph:
 119     typedef int vertex_descriptor; /* vertex number, -1 => none */
 120     typedef int edge_descriptor; /* see above */
 121     typedef undirected_tag directed_category;
 122     typedef disallow_parallel_edge_tag edge_parallel_category;
 123     typedef layout_graph_traversal_category traversal_category;
 124
 125     // Concept IncidenceGraph:
 126     typedef OutEdgeIterator out_edge_iterator;
 127     typedef unsigned degree_size_type;
 128
 129     // Concept VertexListGraph:
 130     typedef VertexIterator vertex_iterator;
 131     typedef unsigned vertices_size_type;
 132   };
 133
 134   // Concept Graph:
 135   inline int null_vertex() { return -1; }
 136
 137   // Concept IncidenceGraph:
 138   inline int source(int f, const Graph&) { return f&VMASK; }
 139   inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
 140   inline std::pair<OutEdgeIterator,OutEdgeIterator>
 141   out_edges(int v, const Graph&) {
 142     return std::make_pair(OutEdgeIterator(v, OutEdgeIterator::voe_min(v)),
 143                           OutEdgeIterator(v, OutEdgeIterator::voe_max(v)));
 144   }
 145   inline unsigned out_degree(int v, const Graph&) {
 146     return OutEdgeIterator::voe_degree(v);
 147   }
 148
 149   // Concept VertexListGraph:
 150   inline std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
 151     return std::make_pair(VertexIterator(0), VertexIterator(N));
 152   }
 153   inline unsigned num_vertices(const Graph&) { return N; }
 154 };
 155
 156 static void single_source_shortest_paths(int v1,
 157                                          const double edge_weights[/*f*/],
 158                                          double vertex_distances[/*v*/]) {
 159   Graph g;
 160
 161   dijkstra_shortest_paths(g, v1,
 162      weight_map(edge_weights).
 163      vertex_index_map(identity_property_map()).
 164      distance_map(vertex_distances));
 165 }
 166
 167 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
 168   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
 169    * along edges, and actual distance d_ij = |vi-vj|.
 170    *
 171    * We will also use the `vertex areas': for each vertex vi the
 172    * vertex area a_vi is the mean area of the incident triangles.
 173    * This is computed elsewhere.
 174    *
 175    * Energy contribution is proportional to
 176    *
 177    *               -4          2
 178    *    a  a   .  d   . [ (s/d)  - 1 ]
 179    *     vi vj
 180    *
 181    * (In practice we compute d^2+epsilon and use it for the
 182    *  divisions, to avoid division by zero.)
 183    */
 184   static const double d2_epsilon= 1e-6;
 185
 186   double edge_weights[N*V6], vertex_distances[N], total_cost=0;
 187   int v1,v2,e,f;
 188
 189   FOR_VERTEX(v1)
 190     FOR_VEDGE_X(v1,e,v2,
 191                 f= v1 | e << ESHIFT,
 192                 edge_weights[f]= NAN)
 193       edge_weights[f]= hypotD(v[v1], v[v2]);
 194
 195   FOR_VERTEX(v1) {
 196     double a1= vertex_areas[v1];
 197     single_source_shortest_paths(v1, edge_weights, vertex_distances);
 198     FOR_VERTEX(v2) {
 199       double a2= vertex_areas[v2];
 200       double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
 201       double sd= vertex_distances[v2] / d2;
 202       double sd2= sd*sd;
 203       total_cost += a1*a2 * (sd2 - 1) / (d2*d2);
 204     }
 205   }
 206   return total_cost;
 207 }