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 using namespace boost;
  51
  52 /*
  53  * We iterate over edges in the following order:
  54  *
  55  *          \#0  /1#
  56  *           \  /
  57  *        ___ 0   __
  58  *       #2    1   #3
  59  *           /  \
  60  *        #4/  #5\     and finally #6 is V6
  61  *
  62  *
  63  * This ordering permits the order-4 nodes at the strip's edge
  64  * to have a contiguous edge iterator values.  The iterator
  65  * starts at #0 which is edge 2 (see mgraph.h), or #2 (edge 3).
  66  */
  67 static const int oei_edge_delta[V6]=
  68   /*     0       1       2       3       4       5      initial e
  69    *    #3      #1      #0      #2      #4      #5      initial ix
  70    *    #4      #2      #1      #3      #5      #6      next ix
  71    *     4       3       1       0       5      V6      next e
  72    */ {
  73      4<<ESHIFT, 2<<ESHIFT, -1<<ESHIFT,
  74                               -3<<ESHIFT, 1<<ESHIFT, (V6-5)<<ESHIFT
  75 };
  76
  77 class OutEdgeIterator :
  78   public iterator_facade<
  79     OutEdgeIterator,
  80     int const,
  81     forward_traversal_tag
  82 > {
  83   int f;
  84  public:
  85   void increment() {
  86     //printf("incrementing f=%03x..",f);
  87     f += oei_edge_delta[f>>ESHIFT];
  88     //printf("%03x\n",f);
  89   }
  90   bool equal(OutEdgeIterator const& other) const { return f == other.f; }
  91   int const& dereference() const { return f; }
  92   OutEdgeIterator() { }
  93   OutEdgeIterator(int _f) : f(_f) { }
  94   OutEdgeIterator(int v, int e) : f(e<<ESHIFT | v) {
  95     //printf("constructed v=%02x e=%x f=%03x\n",v,e,f);
  96   }
  97
  98   static int voe_min(int _v) { return (_v & YMASK) ? 2 : 3; }
  99   static int voe_max(int _v) { return (_v & YMASK)==(Y-1) ? V6 : 4; }
 100   static int voe_degree(int _v) { return RIM_VERTEX_P(_v) ? 4 : V6; }
 101 };
 102
 103 typedef counting_iterator<int> VertexIterator;
 104
 105 namespace boost {
 106   class Graph { }; // this is a dummy as our graph has no actual representation
 107
 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   template<>
 119   struct graph_traits<Graph> {
 120     // Concept Graph:
 121     typedef int vertex_descriptor; /* vertex number, -1 => none */
 122     typedef int edge_descriptor; /* see above */
 123     typedef undirected_tag directed_category;
 124     typedef disallow_parallel_edge_tag edge_parallel_category;
 125     typedef layout_graph_traversal_category traversal_category;
 126
 127     // Concept IncidenceGraph:
 128     typedef OutEdgeIterator out_edge_iterator;
 129     typedef unsigned degree_size_type;
 130
 131     // Concept VertexListGraph:
 132     typedef VertexIterator vertex_iterator;
 133     typedef unsigned vertices_size_type;
 134   };
 135
 136   // Concept Graph:
 137   inline int null_vertex() { return -1; }
 138
 139   // Concept IncidenceGraph:
 140   inline int source(int f, const Graph&) { return f&VMASK; }
 141   inline int target(int f, const Graph&) {
 142     int v2= EDGE_END2(f&VMASK, f>>ESHIFT);
 143     //printf("traversed %03x..%02x\n",f,v2);
 144     return v2;
 145   }
 146   inline std::pair<OutEdgeIterator,OutEdgeIterator>
 147   out_edges(int v, const Graph&) {
 148     return std::make_pair(OutEdgeIterator(v, OutEdgeIterator::voe_min(v)),
 149                           OutEdgeIterator(v, OutEdgeIterator::voe_max(v)));
 150   }
 151   inline unsigned out_degree(int v, const Graph&) {
 152     return OutEdgeIterator::voe_degree(v);
 153   }
 154
 155   // Concept VertexListGraph:
 156   inline
 157   std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
 158     return std::make_pair(VertexIterator(0), VertexIterator(N));
 159   }
 160   inline unsigned num_vertices(const Graph&) { return N; }
 161 };
 162
 163 static void single_source_shortest_paths(int v1,
 164                                          const double edge_weights[/*f*/],
 165                                          double vertex_distances[/*v*/]) {
 166   Graph g;
 167
 168   dijkstra_shortest_paths(g, v1,
 169      weight_map(edge_weights).
 170      vertex_index_map(identity_property_map()).
 171      distance_map(vertex_distances));
 172 }
 173
 174 static int distances[N][N];
 175
 176 void graph_layout_prepare() {
 177   Graph g;
 178   int v1, v2;
 179
 180   FOR_VERTEX(v1) {
 181     int *d= distances[v1];
 182     FOR_VERTEX(v2) d[v2]= -1;
 183     d[v1]= 0;
 184     breadth_first_search
 185       (g, v1,
 186        vertex_index_map(identity_property_map()).
 187        visitor(make_bfs_visitor(record_distances(d,on_tree_edge()))));
 188     printf("%02x:",v1);
 189     FOR_VERTEX(v2) printf(" %02x:%d",v2,d[v2]);
 190     putchar('\n');
 191   }
 192   printf("---\n");
 193   FOR_VERTEX(v1) {
 194     int *d= distances[v1];
 195     printf("%02x:",v1);
 196     FOR_VERTEX(v2) printf(" %02x:%d",v2,d[v2]);
 197     putchar('\n');
 198   }
 199 }
 200
 201 double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
 202   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
 203    * along edges, and actual distance d_ij = |vi-vj|.
 204    *
 205    * We will also use the `vertex areas': for each vertex vi the
 206    * vertex area a_vi is the mean area of the incident triangles.
 207    * This is computed elsewhere.
 208    *
 209    * Energy contribution is proportional to
 210    *
 211    *               -4          2
 212    *    a  a   .  d   . [ (s/d)  - 1 ]
 213    *     vi vj
 214    *
 215    * (In practice we compute d^2+epsilon and use it for the
 216    *  divisions, to avoid division by zero.)
 217    */
 218   static const double d2_epsilon= 1e-6;
 219
 220   //  double edge_weights[V6<<ESHIFT], vertex_distances[N],
 221   double total_cost=0;
 222   int v1,v2,e, nedges=0;
 223   double totaledgelength=0, meanedgelength;
 224
 225   FOR_EDGE(v1,e,v2) {
 226     totaledgelength += hypotD(v[v1], v[v2]);
 227     nedges++;
 228   }
 229
 230   meanedgelength= totaledgelength / nedges;
 231
 232   FOR_VERTEX(v1) {
 233     FOR_VERTEX(v2) {
 234       if (v1 == v2) continue;
 235
 236       double d= hypotD(v[v1],v[v2]);
 237
 238       int dist= distances[v1][v2];
 239       assert(dist>=0);
 240
 241       double s= dist * meanedgelength * 0.03;
 242
 243       double enoughdistance= d - s;
 244       if (enoughdistance > 1e-6) continue;
 245
 246       /* energy = 1/2 stiffness deviation^2
 247        * where stiffness = 1/d
 248        */
 249
 250       double cost= pow(enoughdistance,4);
 251
 252       //double s2= s*s + d2_epsilon;
 253       //double sd2= s2 / d2;
 254       //double cost_contrib= a1*a2 * (sd2 - 1) / (d2*d2);
 255       //double cost_contrib= sd2;
 256
 257       printf("layout %03x..%03x dist=%d mean=%g s=%g d=%g enough=%g"
 258                " cost+=%g\n", v1,v2, dist, meanedgelength,
 259              s,d, enoughdistance, cost);
 260       total_cost += cost;
 261     }
 262   }
 263   return total_cost;
 264 }