X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=bgl.cpp;h=ae36ea4abb5d6ecf8a8027a00598aff5048a95c4;hp=3e2912ebb07bdd56c05c3a23f6104c2aa025ed06;hb=bb0e48bffe4373fabc17a954bbae651bf41def2d;hpb=f9088d46488d2241f9a8e4a15e0cf4cecae02785;ds=sidebyside

diff --git a/bgl.cpp b/bgl.cpp
index 3e2912e..ae36ea4 100644
--- a/bgl.cpp
+++ b/bgl.cpp
@@ -1,4 +1,23 @@
+/*
+ * Everything that needs the Boost Graph Library and C++ templates etc.
+ * (and what a crazy set of stuff that all is)
+ */
+
+#include <math.h>
+
+#include <iterator>
+
+#include <boost/config.hpp>
+#include <boost/iterator/iterator_facade.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/graph_concepts.hpp>
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+#include <boost/graph/properties.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+#include <boost/iterator/iterator_categories.hpp>
+
 extern "C" {
+#include "bgl.h"
 #include "mgraph.h"
 }
 
@@ -13,73 +32,140 @@ extern "C" {
  */
 
 /*
- * We use BGL's implementation of Johnson All Pairs Shortest Paths
+ * We use BGL's implementation of Dijkstra's single source shortest
+ * paths.  We really want all pairs shortest paths, so Johnson All
+ * Pairs Shortest Paths would seem sensible.  But actually Johnson's
+ * algorithm is just a wrapper around Dijkstra's; the extra
+ * functionality is just to deal with -ve edge weights, which we don't
+ * have.  So we can use Dijkstra directly and save some cpu (and some
+ * code: we don't have to supply all of the machinery needed for
+ * Johnson's invocation of Bellman-Ford).  The overall time cost is
+ * O(VE log V); I think the space used is O(E).
  */
 
 #define VMASK (YMASK|XMASK)
 #define ESHIFT (YBITS|XBITS)
-#define FMAX ((5 << ESHIFT) | VMASK)
+
+class Graph { }; // this is a dummy as our graph has no actual representation
+
+using namespace boost;
+
+struct OutEdgeIterator :
+  public iterator_facade<
+    OutEdgeIterator,
+    int const,
+    forward_traversal_tag
+> {
+  int f;
+  void increment() { f += 1<<ESHIFT; }
+  bool equal(OutEdgeIterator const& other) const { return f == other.f; }
+  int const& dereference() const { return f; }
+  OutEdgeIterator() { }
+  OutEdgeIterator(int _f) : f(_f) { }
+  OutEdgeIterator(int v, int e) : f(e << ESHIFT | v) { }
+};
+ 
+typedef counting_iterator<int> VertexIterator;
 
 namespace boost {
-  // We make Layout a model of various BGL Graph concepts.
-  // This mainly means that graph_traits<Layout> has lots of stuff.
+  // We make Graph a model of various BGL Graph concepts.
+  // This mainly means that graph_traits<Graph> has lots of stuff.
 
   // First, some definitions used later:
   
   struct layout_graph_traversal_category : 
     public virtual incidence_graph_tag,
-    public virtual edge_list_graph_tag
-    public virtual vertex_list_graph_tag { };
-  
-  struct OutEdgeIncrable {
-    int f;
-    OutEdgeIncrable operator++(OutEdgeIncrable f) { return f + 1<<ESHIFT; }
-    OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
-  }
-    
-  struct graph_traits<Layout> {
+    public virtual vertex_list_graph_tag,
+    public virtual edge_list_graph_tag { };
 
+  struct graph_traits<Graph> {
     // Concept Graph:
-  
     typedef int vertex_descriptor; /* vertex number, -1 => none */
     typedef int edge_descriptor; /* see above */
     typedef undirected_tag directed_category;
-    typedef disallow_parallel_ege edge_parallel_category;
+    typedef disallow_parallel_edge_tag edge_parallel_category;
     typedef layout_graph_traversal_category traversal_category;
-    inline int null_vertex() { return -1; }
 
     // Concept IncidenceGraph:
-    
-    typedef counting_iterator<OutEdgeIncrable,
-      forward_iterator_tag> out_edge_iterator;
-    typedef int degree_size_type;
-    
-    inline int source(int f, const Layout& g) { return f&VMASK; }
-    inline int target(int f, const Layout& g) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
-    inline std::pair<typename graph_traits<Layout>::out_edge_iterator,
-                     typename graph_traits<Layout>::out_edge_iterator>
-    out_edges(int v, const Layout& g) {
-      return std::make_pair(out_edge_iterator(OutEdgeIncrable(v, VE_MIN(v))),
-			    out_edge_iterator(OutEdgeIncrable(v, VE_MAX(v))));
-    }
-    out_degree(int v, const Layout& g) { return VE_MAX(v) - VE_MIN(v); }
+    typedef OutEdgeIterator out_edge_iterator;
+    typedef unsigned degree_size_type;
 
     // Concept VertexListGraph:
-    typedef counting_iterator<int> vertex_iterator;
+    typedef VertexIterator vertex_iterator;
+    typedef unsigned vertices_size_type;
+  };
     
-}
+  // Concept Graph:
+  inline int null_vertex() { return -1; }
 
-void calculate_layout_energy(const Layout*) {
-  
-  FOR_VERTEX(v1) {
-    boost::dijkstra_shortest_paths(g, v1, 0);
-    
-			    /* weight_map(). ? */
-			    /* vertex_index_map(vimap). */
+  // Concept IncidenceGraph:
+  inline int source(int f, const Graph&) { return f&VMASK; }
+  inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
+  inline std::pair<OutEdgeIterator,OutEdgeIterator>
+  out_edges(int v, const Graph&) {
+    return std::make_pair(OutEdgeIterator(v, VE_MIN(v)),
+			  OutEdgeIterator(v, VE_MAX(v)));
+  }
+  inline unsigned out_degree(int v, const Graph&) {
+    return VE_MAX(v) - VE_MIN(v);
+  }
 
-			    
-			    predecessor_map().
-			    distance_map()
+  // Concept VertexListGraph:
+  inline std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
+    return std::make_pair(VertexIterator(0), VertexIterator(N));
+  }
+  inline unsigned num_vertices(const Graph&) { return N; }
+};
 
+static void single_source_shortest_paths(int v1,
+					 const double edge_weights[/*f*/],
+					 double vertex_distances[/*v*/]) {
+  Graph g;
 
+  dijkstra_shortest_paths(g, v1,
+     weight_map(edge_weights).
+     vertex_index_map(identity_property_map()).
+     distance_map(vertex_distances));
+}
+    
+double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
+  /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
+   * along edges, and actual distance d_ij = |vi-vj|.
+   *
+   * We will also use the `vertex areas': for each vertex vi the
+   * vertex area a_vi is the mean area of the incident triangles.
+   * This is computed elsewhere.
+   *
+   * Energy contribution is proportional to
+   *
+   *               -4          2
+   *    a  a   .  d   . [ (s/d)  - 1 ]
+   *     vi vj
+   *
+   * (In practice we compute d^2+epsilon and use it for the
+   *  divisions, to avoid division by zero.)
+   */
+  static const double d2_epsilon= 1e-6;
+  
+  double edge_weights[N*V6], vertex_distances[N], total_cost=0;
+  int v1,v2,e,f;
 
+  FOR_VERTEX(v1)
+    FOR_VEDGE_X(v1,e,v2,
+		f= v1 | e << ESHIFT,
+		edge_weights[f]= NAN)
+      edge_weights[f]= hypotD(v[v1], v[v2]);
+
+  FOR_VERTEX(v1) {
+    double a1= vertex_areas[v1];
+    single_source_shortest_paths(v1, edge_weights, vertex_distances);
+    FOR_VERTEX(v2) {
+      double a2= vertex_areas[v2];
+      double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
+      double sd= vertex_distances[v2] / d2;
+      double sd2= sd*sd;
+      total_cost += a1*a2 * (sd2 - 1) / (d2*d2);
+    }
+  }
+  return total_cost;
+}