do not actually compute vertex areas

[moebius2.git] / energy.c

diff --git a/energy.c b/energy.c
index 2bbf0492a5a41312d6819907995857bce9d1e180..e9b93fd2577c1daca16b850baebab3d99e410eb6 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -23,14 +23,15 @@ double compute_energy(const struct Vertices *vs) {
   compute_vertex_areas(vs->a);
   energy= 0;
 
-  printing= printing_check(pr_cost);
+  printing= printing_check(pr_cost,0);
 
   if (printing) printf("cost > energy |");
 
-  COST(1e2, edgewise_vertex_displacement_cost(vs->a));
-  COST(1e2, graph_layout_cost(vs->a));
+  COST(3e2, edgewise_vertex_displacement_cost(vs->a));
   COST(1e3, edge_length_variation_cost(vs->a));
-//  COST(1e6, noncircular_rim_cost(vs->a));
+  COST(0.2e3, rim_proximity_cost(vs->a));
+//  COST(1e2, graph_layout_cost(vs->a));
+  COST(1e8, noncircular_rim_cost(vs->a));
 
   if (printing) printf("| total %# e |", energy);
 
@@ -57,7 +58,7 @@ double compute_energy(const struct Vertices *vs) {
 
 static void addcost(double *energy, double tweight, double tcost, int pr) {
   double tenergy= tweight * tcost;
-  if (pr) printf(" %# e > %# e |", tcost, tenergy);
+  if (pr) printf(" %# e > %# e* |", tcost, tenergy);
   *energy += tenergy;
 }
 
@@ -71,7 +72,8 @@ void compute_edge_lengths(const Vertices vertices) {
 }
 
 void compute_vertex_areas(const Vertices vertices) {
-  int v0,v1,v2, e1,e2, k;
+  int v0,v1,v2, e1,e2;
+//  int k;
 
   FOR_VERTEX(v0) {
     double total= 0.0, edges_total=0;
@@ -84,13 +86,13 @@ void compute_vertex_areas(const Vertices vertices) {
 
       edges_total += edge_lengths[v0][e1];
 
-      double e1v[D3], e2v[D3], av[D3];
-      K {
-       e1v[k]= vertices[v1][k] - vertices[v0][k];
-       e2v[k]= vertices[v2][k] - vertices[v0][k];
-      }
-      xprod(av, e1v, e2v);
-      total += magnD(av);
+//      double e1v[D3], e2v[D3], av[D3];
+//      K {
+//     e1v[k]= vertices[v1][k] - vertices[v0][k];
+//     e2v[k]= vertices[v2][k] - vertices[v0][k];
+//      }
+//      xprod(av, e1v, e2v);
+//      total += magnD(av);
       
       count++;
     }
@@ -181,22 +183,48 @@ double edge_length_variation_cost(const Vertices vertices) {
   return cost;
 }    
 
+/*---------- rim proximity cost ----------*/
+
+static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) {
+  /* By symmetry, nearest point on circle is the one with
+   * the same angle subtended at the z axis. */
+  oncircle[0]= p[0];
+  oncircle[1]= p[1];
+  oncircle[2]= 0;
+  double mult= 1.0/ magnD(oncircle);
+  oncircle[0] *= mult;
+  oncircle[1] *= mult;
+}
+
+double rim_proximity_cost(const Vertices vertices) {
+  double oncircle[3], cost=0;
+  int v;
+
+  FOR_VERTEX(v) {
+    int y= v >> YSHIFT;
+    int nominal_edge_distance= y <= Y/2 ? y : Y-1-y;
+    if (nominal_edge_distance==0) continue;
+
+    find_nearest_oncircle(oncircle, vertices[v]);
+
+    cost +=
+      vertex_mean_edge_lengths[v] *
+      (nominal_edge_distance*nominal_edge_distance) /
+      (hypotD2(vertices[v], oncircle) + 1e-6);
+  }
+  return cost;
+}
+
 /*---------- noncircular rim cost ----------*/
 
 double noncircular_rim_cost(const Vertices vertices) {
   int vy,vx,v;
   double cost= 0.0;
+  double oncircle[3];
 
   FOR_RIM_VERTEX(vy,vx,v) {
-    double oncircle[3];
-    /* By symmetry, nearest point on circle is the one with
-     * the same angle subtended at the z axis. */
-    oncircle[0]= vertices[v][0];
-    oncircle[1]= vertices[v][1];
-    oncircle[2]= 0;
-    double mult= 1.0/ magnD(oncircle);
-    oncircle[0] *= mult;
-    oncircle[1] *= mult;
+    find_nearest_oncircle(oncircle, vertices[v]);
+    
     double d2= hypotD2(vertices[v], oncircle);
     cost += d2*d2;
   }