chiark / gitweb /
allow pausing updates
[moebius2.git] / 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 }