chiark / gitweb /
done graph layout cost function core, now reorganise
[moebius2.git] / bgl.cpp
1 extern "C" {
2 #include "mgraph.h"
3 }
4
5 /*
6  * edge descriptor f =   00 | e | y     | x
7  *                            3  YBITS   XBITS
8  *
9  * e is 0..5.  The edge is edge e out of vertex (x,y).
10  *
11  * BGL expects an undirected graph's edges to have two descriptors
12  * each, one in each direction.
13  */
14
15 /*
16  * We use BGL's implementation of Dijkstra's single source shortest
17  * paths.  We really want all pairs shortest paths, so Johnson All
18  * Pairs Shortest Paths would seem sensible.  But actually Johnson's
19  * algorithm is just a wrapper around Dijkstra's; the extra
20  * functionality is just to deal with -ve edge weights, which we don't
21  * have.  So we can use Dijkstra directly and save some cpu (and some
22  * code: we don't have to supply all of the machinery needed for
23  * Johnson's invocation of Bellman-Ford).  The overall time cost is
24  * O(VE log V); I think the space used is O(E).
25  */
26
27 #define VMASK (YMASK|XMASK)
28 #define ESHIFT (YBITS|XBITS)
29
30 namespace boost {
31   // We make Layout a model of various BGL Graph concepts.
32   // This mainly means that graph_traits<Layout> has lots of stuff.
33
34   // First, some definitions used later:
35   
36   struct layout_graph_traversal_category : 
37     public virtual incidence_graph_tag,
38     public virtual vertex_list_graph_tag,
39     public virtual edge_list_graph_tag { };
40   
41   struct OutEdgeIncrable {
42     int f;
43     OutEdgeIncrable& operator++() { f += 1<<ESHIFT; return self; }
44     OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
45   };
46  
47   struct graph_traits<Layout> {
48
49     // Concept Graph:
50   
51     typedef int vertex_descriptor; /* vertex number, -1 => none */
52     typedef int edge_descriptor; /* see above */
53     typedef undirected_tag directed_category;
54     typedef disallow_parallel_ege edge_parallel_category;
55     typedef layout_graph_traversal_category traversal_category;
56     inline int null_vertex() { return -1; }
57
58     // Concept IncidenceGraph:
59     
60     typedef counting_iterator<OutEdgeIncrable,
61       forward_iterator_tag> out_edge_iterator;
62     typedef int degree_size_type;
63     
64     inline int source(int f, const Layout&) { return f&VMASK; }
65     inline int target(int f, const Layout&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
66     inline std::pair<out_edge_iterator,out_edge_iterator>
67     out_edges(int v, const Layout&) {
68       return std::make_pair(out_edge_iterator(OutEdgeIncrable(v, VE_MIN(v))),
69                             out_edge_iterator(OutEdgeIncrable(v, VE_MAX(v))));
70     }
71     out_degree(int v, const Layout&) { return VE_MAX(v) - VE_MIN(v); }
72
73     // Concept VertexListGraph:
74     typedef counting_iterator<int> vertex_iterator;
75     typedef unsigned vertices_size_type;
76     inline std::pair<vertex_iterator,vertex_iterator>
77     vertices(const Layout&) {
78       return std::make_pair(vertex_iterator(0), vertex_iterator(N));
79     }
80     inline unsigned num_vertices(const Layout&) { return N; }
81
82 }
83
84 static void single_source_shortest_paths(int v1,
85                                          const double edge_weights[/*f*/],
86                                          double vertex_distances[/*v*/]) {
87   boost::dijkstra_shortest_paths
88     (g, v1,
89      weight_map(edge_weights).
90      vertex_index_map(identity_property_map()).
91      distance_map(vertex_distances));
92 }
93     
94 double graph_layout_cost(const Layout *g, const double vertex_areas[N]) {
95   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
96    * along edges, and actual distance d_ij = |vi-vj|.
97    *
98    * We will also use the `vertex areas': for each vertex vi the
99    * vertex area a_vi is the mean area of the incident triangles.
100    * This is computed elsewhere.
101    *
102    * Energy contribution is proportional to
103    *
104    *               -4          2
105    *    a  a   .  d   . [ (s/d)  - 1 ]
106    *     vi vj
107    *
108    * (In practice we compute d^2+epsilon and use it for the
109    *  divisions, to avoid division by zero.)
110    */
111   double edge_weights[N*V6], vertex_distances[N], total_cost;
112   int v1, e, f;
113
114   FOR_VEDGE_X(v1,e,
115               f= v1 | e << ESHIFT,
116               edge_weights[f]= NaN)
117     edge_weights[f]= hypotD(g.v[v1], g.v[v2]);
118
119   FOR_VERTEX(v1) {
120     double a1= vertex_areas[v1];
121     single_source_shortest_paths(v1, edge_weights, vertex_distances);
122     FOR_VERTEX(v2) {
123       double a2= vertex_areas[v2];
124       double d2= hypotD2plus(g->v[v1],g->v[v2], d2_epsilon);
125       double sd= vertex_distances[v2] / d2;
126       double sd2= sd*sd;
127       total_cost= fma_fast(a1*a2, (sd2 - 1.0)/(d2*d2), total_cost);
128     }
129   }
130   return total_cost;
131 }