chiark / gitweb /
compiles and most output facets are right; needs debugging
[moebius2.git] / graph.c
1 /*
2  * graph layout energy
3  */
4
5 #include "mgraph.h"
6 #include "minimise.h"
7
8 static int sqdistances[N][N];
9
10 static double alpha, beta, beta_prime;
11
12 static void breadth_first_search(int start, int sqdistances_r[N]) {
13   int d[N], buffer[N], *buf_pop=buffer, *buf_push=buffer;
14   int v,e, current, future, dfuture;
15
16   buf_push= buf_pop= buffer;
17   FOR_VERTEX(v) d[v]= -1;
18
19   d[start]= 0;
20   *buf_push++= start;
21
22   while (buf_pop < buf_push) {
23     current= *buf_pop++;
24     dfuture= d[current] + 1;
25     FOR_VEDGE(current,e,future) {
26       if (d[future] >= 0) continue; /* already found this one */
27       d[future]= dfuture;
28       *buf_push++= future;
29     }
30   }
31   assert(buf_pop==buf_push);
32   assert(buf_push <= buffer+sizeof(buffer)/sizeof(buffer[0]));
33
34   FOR_VERTEX(v) {
35     assert(d[v] >= 0);
36     sqdistances_r[v]= d[v] * d[v];
37   }
38
39 //  FOR_VERTEX(v) {
40 //    
41 }
42
43 void graph_layout_prepare() {
44   int v1;
45   
46   FOR_VERTEX(v1)
47     breadth_first_search(v1, sqdistances[v1]);
48
49   alpha= 2;
50   beta= log(10)/log(alpha);
51   beta_prime= (1-beta)/2;
52   printf("alpha=%g beta=%g beta'=%g\n", alpha,beta,beta_prime);
53 }
54
55
56 double graph_layout_cost(const Vertices v) {
57   /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
58    * along edges, and actual distance d_ij = |vi-vj|.
59    *
60    * We will also use the `vertex areas': for each vertex vi the
61    * vertex area a_vi is the mean area of the incident triangles.
62    * This is computed elsewhere.
63    *
64    * Energy contribution is proportional to
65    *
66    *               -4          2
67    *    a  a   .  d   . [ (s/d)  - 1 ]
68    *     vi vj
69    *
70    * (In practice we compute d^2+epsilon and use it for the
71    *  divisions, to avoid division by zero.)
72    */
73   //static const double d2_epsilon= 1e-6;
74
75   double total_cost=0;
76   int v1,v2,e, nedges=0;
77   double totaledgelength=0, meanedgelength, meanedgelength2;
78
79   FOR_EDGE(v1,e,v2) {
80     totaledgelength += hypotD(v[v1], v[v2]);
81     nedges++;
82   }
83
84   meanedgelength= totaledgelength / nedges;
85   meanedgelength2= meanedgelength * meanedgelength;
86 //  printf("mean=%g mean^2=%g\n", meanedgelength, meanedgelength2);
87     
88   FOR_VERTEX(v1) {
89     FOR_VERTEX(v2) {
90       if (v1 == v2) continue;
91
92       double d2= hypotD2(v[v1],v[v2]);
93       
94       int dist2= sqdistances[v1][v2];
95       assert(dist2>0);
96
97       double s2= dist2 * meanedgelength2;
98
99       /* energy = (d/s)^(1-beta)     where beta is log\_{alpha}(10)
100        * energy = ((d/s)^2) ^ (1-beta)/2
101        *          let beta' = (1-beta)/2
102        */
103
104       double cost= (vertex_mean_edge_lengths[v1]) * pow(d2/s2, beta_prime);
105       
106       //printf("layout %03x..%03x dist^2=%d s^2=%g d^2=%g "
107                //" cost+=%g\n", v1,v2, dist2,
108                //            s2,d2, cost);
109       total_cost += cost;
110     }
111   }
112   return total_cost/meanedgelength;
113 }