X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=bgl.cpp;h=c35110b0682be288013d45eca339eb59f1a3f503;hp=774be825a1cc53bc017884e94876ddd49c134984;hb=97ef96dd318b684d829ac950403582788df427d8;hpb=0cd9cfd406d72d3d40941f87937b60f0d10f1d4d

diff --git a/bgl.cpp b/bgl.cpp
index 774be82..c35110b 100644
--- a/bgl.cpp
+++ b/bgl.cpp
@@ -1,15 +1,35 @@
+/*
+ * 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"
 }
 
 /*
- * edge descriptor f =   00 | e | y     | x
- *                            3  YBITS   XBITS
+ * edge descriptor f =   0000 | e | y     | x
+ *                              3  YBITS   XBITS
  *
- * e is 0..5.  The edge is edge e out of vertex (x,y).
+ * e is 0..6.  The edge is edge e out of vertex (x,y), or if
+ * e==6 it's the `at end' value for the out edge iterator.
  *
  * BGL expects an undirected graph's edges to have two descriptors
- * each, one in each direction.
+ * each, one in each direction (otherwise e would be just 0..2).
  */
 
 /*
@@ -25,11 +45,68 @@ extern "C" {
  */
 
 #define VMASK (YMASK|XMASK)
-#define ESHIFT (YBITS|XBITS)
+#define ESHIFT (YBITS+XBITS)
+
+using namespace boost;
+
+/*
+ * We iterate over edges in the following order:
+ *
+ *          \#0  /1#
+ *	     \  /
+ *	  ___ 0   __
+ *	 #2    1   #3
+ *	     /  \
+ *	  #4/  #5\     and finally #6 is V6
+ *
+ *
+ * This ordering permits the order-4 nodes at the strip's edge
+ * to have a contiguous edge iterator values.  The iterator
+ * starts at #0 which is edge 2 (see mgraph.h), or #2 (edge 3).
+ */
+static const int oei_edge_delta[V6]=
+  /*	 0	 1	 2	 3	 4	 5	initial e
+   *	#3	#1	#0	#2	#4	#5	initial ix
+   *	#4	#2	#1	#3	#5	#6	next ix
+   *	 4	 3	 1	 0	 5	V6	next e
+   */ {
+     4<<ESHIFT, 2<<ESHIFT, -1<<ESHIFT,
+                              -3<<ESHIFT, 1<<ESHIFT, (V6-5)<<ESHIFT
+};
+
+class OutEdgeIterator :
+  public iterator_facade<
+    OutEdgeIterator,
+    int const,
+    forward_traversal_tag
+> {
+  int f;
+ public:
+  void increment() {
+    //printf("incrementing f=%03x..",f);
+    f += oei_edge_delta[f>>ESHIFT];
+    //printf("%03x\n",f);
+  }
+  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) {
+    //printf("constructed v=%02x e=%x f=%03x\n",v,e,f);
+  }
+
+  static int voe_min(int _v) { return (_v & YMASK) ? 2 : 3; }
+  static int voe_max(int _v) { return (_v & YMASK)==(Y-1) ? V6 : 4; }
+  static int voe_degree(int _v) { return RIM_VERTEX_P(_v) ? 4 : V6; }
+};
+ 
+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.
+  class Graph { }; // this is a dummy as our graph has no actual representation
+
+  // 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:
   
@@ -37,81 +114,151 @@ namespace boost {
     public virtual incidence_graph_tag,
     public virtual vertex_list_graph_tag,
     public virtual edge_list_graph_tag { };
-  
-  struct OutEdgeIncrable {
-    int f;
-    OutEdgeIncrable& operator++() { f += 1<<ESHIFT; return self; }
-    OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
-  };
- 
-  struct graph_traits<Layout> {
 
+  template<>
+  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&) { return f&VMASK; }
-    inline int target(int f, const Layout&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
-    inline std::pair<out_edge_iterator,out_edge_iterator>
-    out_edges(int v, const Layout&) {
-      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&) { 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;
-    inline std::pair<vertex_iterator,vertex_iterator>
-    vertices(const Layout&) {
-      return std::make_pair(vertex_iterator(0), vertex_iterator(N));
-    }
-    inline unsigned num_vertices(const Layout&) { return N; }
+  };
+    
+  // Concept Graph:
+  inline int null_vertex() { return -1; }
 
-}
+  // Concept IncidenceGraph:
+  inline int source(int f, const Graph&) { return f&VMASK; }
+  inline int target(int f, const Graph&) {
+    int v2= EDGE_END2(f&VMASK, f>>ESHIFT);
+    //printf("traversed %03x..%02x\n",f,v2);
+    return v2;
+  }
+  inline std::pair<OutEdgeIterator,OutEdgeIterator>
+  out_edges(int v, const Graph&) {
+    return std::make_pair(OutEdgeIterator(v, OutEdgeIterator::voe_min(v)),
+			  OutEdgeIterator(v, OutEdgeIterator::voe_max(v)));
+  }
+  inline unsigned out_degree(int v, const Graph&) {
+    return OutEdgeIterator::voe_degree(v);
+  }
 
-struct VertexIndexMap;
+  // 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; }
+};
 
-namespace boost {
-  struct property_traits<VertexIndexMap> {
-    // Concept Readable Property Map:
-    typedef int value_type, reference, key_type;
-    category 
-class Boost
-  }};
-
-void single_source_shortest_paths(int v1,
-				  const double edge_weights[/*f*/],
-				  ) {
-  boost::dijkstra_shortest_paths
-    (g, v1,
+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()).
-     
-    
-void all_pairs_shortest_paths(const Layout *g) {
+     distance_map(vertex_distances));
+}
+
+static int distances[N][N];
+
+void graph_layout_prepare() {
+  Graph g;
+  int v1, v2;
   
   FOR_VERTEX(v1) {
-       
- 0);
+    int *d= distances[v1];
+    FOR_VERTEX(v2) d[v2]= -1;
+    d[v1]= 0;
+    breadth_first_search
+      (g, v1,
+       vertex_index_map(identity_property_map()).
+       visitor(make_bfs_visitor(record_distances(d,on_tree_edge()))));
+    printf("%02x:",v1);
+    FOR_VERTEX(v2) printf(" %02x:%d",v2,d[v2]);
+    putchar('\n');
+  }
+  printf("---\n");
+  FOR_VERTEX(v1) {
+    int *d= distances[v1];
+    printf("%02x:",v1);
+    FOR_VERTEX(v2) printf(" %02x:%d",v2,d[v2]);
+    putchar('\n');
+  }  
+}
+
+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[V6<<ESHIFT], vertex_distances[N],
+  double total_cost=0;
+  int v1,v2,e, nedges=0;
+  double totaledgelength=0, meanedgelength;
+
+  FOR_EDGE(v1,e,v2) {
+    totaledgelength += hypotD(v[v1], v[v2]);
+    nedges++;
+  }
+
+  meanedgelength= totaledgelength / nedges;
     
-			    /* weight_map(). ? */
-			    /* 
+  FOR_VERTEX(v1) {
+    FOR_VERTEX(v2) {
+      if (v1 == v2) continue;
+
+      double d= hypotD(v[v1],v[v2]);
+      
+      int dist= distances[v1][v2];
+      assert(dist>=0);
+
+      double s= dist * meanedgelength * 0.03;
 
-			    
-			    predecessor_map().
-			    distance_map()
+      double enoughdistance= d - s;
+      if (enoughdistance > 1e-6) continue;
 
+      /* energy = 1/2 stiffness deviation^2
+       * where stiffness = 1/d
+       */
 
+      double cost= pow(enoughdistance,4);
+      
+      //double s2= s*s + d2_epsilon;
+      //double sd2= s2 / d2;
+      //double cost_contrib= a1*a2 * (sd2 - 1) / (d2*d2);
+      //double cost_contrib= sd2;
 
+      printf("layout %03x..%03x dist=%d mean=%g s=%g d=%g enough=%g"
+      	       " cost+=%g\n", v1,v2, dist, meanedgelength,
+	     s,d, enoughdistance, cost);
+      total_cost += cost;
+    }
+  }
+  return total_cost;
+}