wip before increase power of bendings

diff --git a/energy.c b/energy.c
index b3af3faae50a8c705c0d2a7295b73e03b6648845..8bc125af99f53df80d45f60e507cea337cf7efa7 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -7,7 +7,11 @@
 #include "mgraph.h"
 #include "parallel.h"
 
-double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
+double vertex_mean_edge_lengths[N];
+
+static double vertex_areas[N];
+static double edge_lengths[N][V6];
+static double rim_vertex_angles[N];
 
 static double best_energy= DBL_MAX;
 
@@ -29,7 +33,8 @@ typedef struct {
 
 static PreComputation *const precomps[]= {
   compute_edge_lengths,
-  compute_vertex_areas
+  compute_vertex_areas,
+  compute_rim_twist_angles
 };
 
 static const CostContribution costs[]= {
@@ -45,11 +50,11 @@ static const CostContribution costs[]= {
 #endif
 
 #if XBITS==4
-#define STOP_EPSILON 1.2e-3
+#define STOP_EPSILON 5e-3
     COST(  3e3,   vertex_displacement_cost)
     COST(  3e3,   vertex_edgewise_displ_cost)
     COST( 0.2e3,  rim_proximity_cost)
-    COST( 1e6,  rim_twist_cost)
+    COST( 1e4,  rim_twist_cost)
     COST(  1e12,   noncircular_rim_cost)
     COST(  10e1,   nonequilateral_triangles_cost)
 //    COST(  1e1,   small_triangles_cost)
@@ -62,8 +67,9 @@ static const CostContribution costs[]= {
     COST(  3e3,   vertex_displacement_cost)
     COST(  3e3,   vertex_edgewise_displ_cost)
     COST(  2e-1,  rim_proximity_cost)
+    COST(  3e3,  rim_twist_cost)
     COST(  1e12,   noncircular_rim_cost)
-    COST(  10e1,   nonequilateral_triangles_cost)
+    COST(   3e2,   nonequilateral_triangles_cost)
 //    COST(  1e1,   small_triangles_cost)
 //    COST(  1e6,   edge_angle_cost)
                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
@@ -74,6 +80,7 @@ static const CostContribution costs[]= {
     COST(  3e3,   vertex_displacement_cost)
     COST(  3e3,   vertex_edgewise_displ_cost)
     COST( 0.02e0,  rim_proximity_cost)
+    COST(  1e4,  rim_twist_cost)
     COST(  1e12,   noncircular_rim_cost)
     COST(  10e1,   nonequilateral_triangles_cost)
 //    COST(  1e1,   small_triangles_cost)
@@ -183,7 +190,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;
 }
 
@@ -464,7 +471,7 @@ static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) {
 }
 
 double rim_proximity_cost(const Vertices vertices, int section) {
-  double oncircle[3], cost=0;
+  double oncircle[D3], cost=0;
   int v;
 
   FOR_VERTEX(v, OUTER) {
@@ -500,42 +507,33 @@ double noncircular_rim_cost(const Vertices vertices, int section) {
 
 /*---------- rim contact angle rotation ----------*/
 
-/*
- *
- *               vq ---- vs
- *              /       /
- *           e0/       /e1
- *           vp 0--- vr
- *
- */
+void compute_rim_twist_angles(const Vertices vertices, int section) {
+  double oncircle[D3], distance[D3];
+  int vpy,vpx,v,k;
+    
+  FOR_NEAR_RIM_VERTEX(vpy,vpx,v, 1,OUTER) {
+    find_nearest_oncircle(oncircle, vertices[v]);
+    /* we are interested in the angle subtended at the rim, from the
+     * rim's point of view. */
+    K distance[k]= vertices[v][k] - oncircle[k];
+    double distance_positive_z= distance[3];
+    double distance_radial_outwards= dotprod(distance, oncircle);
+    rim_vertex_angles[v]= atan2(distance_positive_z, distance_radial_outwards);
+  }
+}
 
 double rim_twist_cost(const Vertices vertices, int section) {
-  const double our_epsilon=1e-6;
-  int vpy,vpx,vpi, vqi,vri,vsi, e0,e1, k;
-  double total_cost= 0.0;
-  double pq[D3], rs[D3], pq_x_rs[D3];
-
-  FOR_RIM_VERTEX(vpy,vpx,vpi, INNER) {
-    FOR_VPEDGE(e0) {
-      if (e0==0 || e0==3) continue;
-      vqi= EDGE_END2(vpi,e0);  if (vqi<0) continue;
-      vri= EDGE_END2(vpi,0);  assert(vri>=0);
-      e1= !vertices_span_join_p(vpi,vri) ? e0 : V6 - e0;
-      vsi= EDGE_END2(vri,e1);  assert(vsi>=0);
-
-      K pq[k]= -vertices[vpi][k] + vertices[vqi][k];
-      K rs[k]= -vertices[vri][k] + vertices[vsi][k];
-
-      xprod(pq_x_rs, pq,rs);
-      double magndiv= edge_lengths[vpi][e0] * edge_lengths[vri][e1];
-
-      double cost= magnD2(pq_x_rs) / (magndiv*magndiv + our_epsilon);
-
-//fprintf(stderr,"rimtwist P=%03x Q=%03x R=%03x S=%03x e0=%d e1=%d %f\n",
-//     vpi,vqi,vri,vsi,e0,e1, cost);
+  double total_cost= 0;
+  int vpy,vpx,v0,v1;
+  
+  FOR_NEAR_RIM_VERTEX(vpy,vpx,v0, 1,OUTER) {
+    v1= EDGE_END2(v0,0);  assert(v1!=0);
+    double delta= rim_vertex_angles[v0] - rim_vertex_angles[v1];
+    if (delta < M_PI) delta += 2*M_PI;
+    if (delta > M_PI) delta -= 2*M_PI;
 
-      total_cost += cost;
-    }
+    double cost= pow(delta, 4);
+    total_cost += cost;
   }
 
   return total_cost;

diff --git a/mgraph.h b/mgraph.h
index 2267fdbda1ffcf59bbf7a9919c651bee46ffa6b0..63fd5dca82ae60e4a78d865bcbb3cdad71a04527 100644 (file)
--- a/mgraph.h
+++ b/mgraph.h
@@ -129,10 +129,13 @@ int edge_reverse(int v1, int e);
   FOR_VERTEX((v1), loop)                       \
     FOR_VEDGE((v1),(e),(v2))
 
-#define FOR_RIM_VERTEX(vy,vx,v, loop)          \
-  for ((vy)=0; (vy)<Y; (vy)+=Y-1)              \
+#define FOR_NEAR_RIM_VERTEX(vy,vx,v, disttorim, loop)          \
+  for ((vy)=(disttorim); (vy)<Y; (vy)+=Y-1-2*(disttorim))      \
     loop ((vx), 0, X, (v)= (vy)<<YSHIFT | (vx))
 
+#define FOR_RIM_VERTEX(vy,vx,v, loop)          \
+  FOR_NEAR_RIM_VERTEX((vy),(vx),(v), 0, loop)
+
 int vertices_span_join_p(int v0, int v1);
 
 typedef double Vertices[N][D3];

diff --git a/minimise.h b/minimise.h
index 4999fd04790bab482f69acba42fa46d74814addc..9ce07c429318949d5f9ca74d42c011ad2937a2a3 100644 (file)
--- a/minimise.h
+++ b/minimise.h
@@ -15,8 +15,10 @@ void graph_layout_prepare();
 
 void compute_vertex_areas(const Vertices vertices, int section);
 void compute_edge_lengths(const Vertices vertices, int section);
+void compute_rim_twist_angles(const Vertices vertices, int section);
 
-extern double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
+//extern double vertex_areas[N], edge_lengths[N][V6];
+extern double vertex_mean_edge_lengths[N];
 
 extern const double edge_angle_cost_circcircrat;