minor improvements

[moebius2.git] / energy.c

diff --git a/energy.c b/energy.c
index a54cfe2f93bf07ab70d0f9bce02639a868bd9b1b..2bbf0492a5a41312d6819907995857bce9d1e180 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -6,7 +6,8 @@
 #include "minimise.h"
 #include "mgraph.h"
 
-static void compute_vertex_areas(const Vertices vertices, double areas[N]);
+double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
+
 static double best_energy= DBL_MAX;
 
 static void addcost(double *energy, double tweight, double tcost, int pr);
@@ -15,10 +16,11 @@ static void addcost(double *energy, double tweight, double tcost, int pr);
 /*---------- main energy computation and subroutines ----------*/
 
 double compute_energy(const struct Vertices *vs) {
-  double vertex_areas[N], energy;
+  double energy;
   int printing;
 
-  compute_vertex_areas(vs->a, vertex_areas);
+  compute_edge_lengths(vs->a);
+  compute_vertex_areas(vs->a);
   energy= 0;
 
   printing= printing_check(pr_cost);
@@ -26,8 +28,9 @@ double compute_energy(const struct Vertices *vs) {
   if (printing) printf("cost > energy |");
 
   COST(1e2, edgewise_vertex_displacement_cost(vs->a));
-  COST(1e2, graph_layout_cost(vs->a,vertex_areas));
-  COST(1e6, noncircular_rim_cost(vs->a));
+  COST(1e2, graph_layout_cost(vs->a));
+  COST(1e3, edge_length_variation_cost(vs->a));
+//  COST(1e6, noncircular_rim_cost(vs->a));
 
   if (printing) printf("| total %# e |", energy);
 
@@ -58,11 +61,20 @@ static void addcost(double *energy, double tweight, double tcost, int pr) {
   *energy += tenergy;
 }
 
-static void compute_vertex_areas(const Vertices vertices, double areas[N]) {
+/*---------- Precomputations ----------*/
+
+void compute_edge_lengths(const Vertices vertices) {
+  int v1,e,v2;
+  
+  FOR_EDGE(v1,e,v2)
+    edge_lengths[v1][e]= hypotD(vertices[v1],vertices[v2]);
+}
+
+void compute_vertex_areas(const Vertices vertices) {
   int v0,v1,v2, e1,e2, k;
 
   FOR_VERTEX(v0) {
-    double total= 0.0;
+    double total= 0.0, edges_total=0;
     int count= 0;
 
     FOR_VEDGE(v0,e1,v1) {
@@ -70,6 +82,8 @@ static void compute_vertex_areas(const Vertices vertices, double areas[N]) {
       v2= EDGE_END2(v0,e2);
       if (v2<0) continue;
 
+      edges_total += edge_lengths[v0][e1];
+
       double e1v[D3], e2v[D3], av[D3];
       K {
        e1v[k]= vertices[v1][k] - vertices[v0][k];
@@ -77,9 +91,11 @@ static void compute_vertex_areas(const Vertices vertices, double areas[N]) {
       }
       xprod(av, e1v, e2v);
       total += magnD(av);
+      
       count++;
     }
-    areas[v0]= total / count;
+    vertex_areas[v0]= total / count;
+    vertex_mean_edge_lengths[v0]= edges_total / count;
   }
 }
 
@@ -151,6 +167,20 @@ double edgewise_vertex_displacement_cost(const Vertices vertices) {
   return total_cost;
 }
 
+/*---------- edge length variation ----------*/
+
+double edge_length_variation_cost(const Vertices vertices) {
+  double diff, cost= 0;
+  int v0, efwd,vfwd, eback;
+
+  FOR_EDGE(v0,efwd,vfwd) {
+    eback= edge_reverse(v0,efwd);
+    diff= edge_lengths[v0][efwd] - edge_lengths[v0][eback];
+    cost += diff*diff;
+  }
+  return cost;
+}    
+
 /*---------- noncircular rim cost ----------*/
 
 double noncircular_rim_cost(const Vertices vertices) {