X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=blobdiff_plain;f=bgl.cpp;h=c35110b0682be288013d45eca339eb59f1a3f503;hp=ae36ea4abb5d6ecf8a8027a00598aff5048a95c4;hb=c3e758476ec10723ccb81a7dd91967d2996665b9;hpb=fe0566cb1f070bcad4d8afa20457d6937137482d

diff --git a/bgl.cpp b/bgl.cpp
index ae36ea4..c35110b 100644
--- a/bgl.cpp
+++ b/bgl.cpp
@@ -22,13 +22,14 @@ extern "C" {
 }
 
 /*
- * 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).
  */
 
 /*
@@ -44,30 +45,66 @@ extern "C" {
  */
 
 #define VMASK (YMASK|XMASK)
-#define ESHIFT (YBITS|XBITS)
-
-class Graph { }; // this is a dummy as our graph has no actual representation
+#define ESHIFT (YBITS+XBITS)
 
 using namespace boost;
 
-struct OutEdgeIterator :
+/*
+ * 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;
-  void increment() { f += 1<<ESHIFT; }
+ 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) { }
+  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 {
+  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.
 
@@ -78,6 +115,7 @@ namespace boost {
     public virtual vertex_list_graph_tag,
     public virtual edge_list_graph_tag { };
 
+  template<>
   struct graph_traits<Graph> {
     // Concept Graph:
     typedef int vertex_descriptor; /* vertex number, -1 => none */
@@ -100,18 +138,23 @@ namespace boost {
 
   // 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 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, VE_MIN(v)),
-			  OutEdgeIterator(v, VE_MAX(v)));
+    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 VE_MAX(v) - VE_MIN(v);
+    return OutEdgeIterator::voe_degree(v);
   }
 
   // Concept VertexListGraph:
-  inline std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
+  inline
+  std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
     return std::make_pair(VertexIterator(0), VertexIterator(N));
   }
   inline unsigned num_vertices(const Graph&) { return N; }
@@ -127,7 +170,34 @@ static void single_source_shortest_paths(int v1,
      vertex_index_map(identity_property_map()).
      distance_map(vertex_distances));
 }
-    
+
+static int distances[N][N];
+
+void graph_layout_prepare() {
+  Graph g;
+  int v1, v2;
+  
+  FOR_VERTEX(v1) {
+    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|.
@@ -146,25 +216,48 @@ double graph_layout_cost(const Vertices v, const double vertex_areas[N]) {
    *  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]);
+  //  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;
+    
   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);
+      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;
+
+      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;