done graph layout cost function core, now reorganise

[moebius2.git] / bgl.cpp

diff --git a/bgl.cpp b/bgl.cpp
index 775fb12141c5d806bc3d98f36f9db2600eec7564..eb7e096955bf507cf1445446237259b9dd63b423 100644 (file)
--- a/bgl.cpp
+++ b/bgl.cpp
@@ -91,18 +91,24 @@ static void single_source_shortest_paths(int v1,
      distance_map(vertex_distances));
 }
     
-double graph_layout_energy(const Layout *g) {
-  /* For each (vi,vj) computes shortest path pij = vi..vj along edges,
-   * and actual distance dij = vi-vj.
+double graph_layout_cost(const Layout *g, 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
    *
-   *     pij - dij   
-   *     ---------  .  Sigma                                 length
-   *           2        e member (edges(vi) union edges(vj)        e
-   *        dij
+   *               -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.)
    */
-  double edge_weights[N*V6], vertex_distances[N];
+  double edge_weights[N*V6], vertex_distances[N], total_cost;
   int v1, e, f;
 
   FOR_VEDGE_X(v1,e,
@@ -111,23 +117,15 @@ double graph_layout_energy(const Layout *g) {
     edge_weights[f]= hypotD(g.v[v1], g.v[v2]);
 
   FOR_VERTEX(v1) {
+    double a1= vertex_areas[v1];
     single_source_shortest_paths(v1, edge_weights, vertex_distances);
     FOR_VERTEX(v2) {
-      
-            
-      
-    
-    { ) {
-    
-       
- 0);
-    
-                           /* weight_map(). ? */
-                           /* 
-
-                           
-                           predecessor_map().
-                           distance_map()
-
-
-
+      double a2= vertex_areas[v2];
+      double d2= hypotD2plus(g->v[v1],g->v[v2], d2_epsilon);
+      double sd= vertex_distances[v2] / d2;
+      double sd2= sd*sd;
+      total_cost= fma_fast(a1*a2, (sd2 - 1.0)/(d2*d2), total_cost);
+    }
+  }
+  return total_cost;
+}