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 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 };
  69
  70 class OutEdgeIterator :
  71   public iterator_facade<
  72     OutEdgeIterator,
  73     int const,
  74     forward_traversal_tag
  75 > {
  76   int v, ix, f;
  77
  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(); }
  85
  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 };
  89
  90 typedef counting_iterator<int> VertexIterator;
  91
  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.
  95
  96   // First, some definitions used later:
  97
  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 { };
 102
 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;
 110
 111     // Concept IncidenceGraph:
 112     typedef OutEdgeIterator out_edge_iterator;
 113     typedef unsigned degree_size_type;
 114
 115     // Concept VertexListGraph:
 116     typedef VertexIterator vertex_iterator;
 117     typedef unsigned vertices_size_type;
 118   };
 119
 120   // Concept Graph:
 121   inline int null_vertex() { return -1; }
 122
 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   }
 134
 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 };
 141
 142 static void single_source_shortest_paths(int v1,
 143                                          const double edge_weights[/*f*/],
 144                                          double vertex_distances[/*v*/]) {
 145   Graph g;
 146
 147   dijkstra_shortest_paths(g, v1,
 148      weight_map(edge_weights).
 149      vertex_index_map(identity_property_map()).
 150      distance_map(vertex_distances));
 151 }
 152
 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;
 171
 172   double edge_weights[N*V6], vertex_distances[N], total_cost=0;
 173   int v1,v2,e,f;
 174
 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]);
 180
 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 }