quite nice

diff --git a/energy.c b/energy.c
index 49d4e434e9243b037905b98afb5a2ae975a2fc7f..4723e0ce408a5510d29644432c8cb096f897123a 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -27,9 +27,10 @@ double compute_energy(const struct Vertices *vs) {
 
   if (printing) printf("cost > energy |");
 
 
   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(1e3, edge_length_variation_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);
   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;
 
 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;
 }
 
   *energy += tenergy;
 }
 


@@ -181,22 +182,48 @@ double edge_length_variation_cost(const Vertices vertices) {
   return cost;
 }    
 
   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;
 /*---------- 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) {
 
   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;
   }
     double d2= hypotD2(vertices[v], oncircle);
     cost += d2*d2;
   }

diff --git a/minimise.h b/minimise.h
index 2043d16249eef6ce6458a0de18030ca20437a52e..f1da76f1103458e54f32af6107511f42d947956c 100644 (file)
--- a/minimise.h
+++ b/minimise.h
@@ -19,6 +19,7 @@ extern double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
 double noncircular_rim_cost(const Vertices vertices);
 double edgewise_vertex_displacement_cost(const Vertices vertices);
 double edge_length_variation_cost(const Vertices vertices);
 double noncircular_rim_cost(const Vertices vertices);
 double edgewise_vertex_displacement_cost(const Vertices vertices);
 double edge_length_variation_cost(const Vertices vertices);

+double rim_proximity_cost(const Vertices vertices);

 
 extern const char *input_file, *output_file;
 extern char *output_file_tmp;
 
 extern const char *input_file, *output_file;
 extern char *output_file_tmp;