chiark / gitweb /
wip energy functions
authorIan Jackson <ian@turbine>
Sun, 19 Oct 2008 21:37:09 +0000 (22:37 +0100)
committerIan Jackson <ian@turbine>
Sun, 19 Oct 2008 21:37:09 +0000 (22:37 +0100)
energy.c
minimise.h

index 19f32e2..72c6fd3 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -45,30 +45,42 @@ static const CostContribution costs[]= {
 #endif
 
 #if XBITS==4
-#define STOP_EPSILON 1.2e-4
+#define STOP_EPSILON 1.2e-3
     COST(  3e3,   vertex_displacement_cost)
+    COST(  3e3,   vertex_edgewise_displ_cost)
     COST( 0.2e3,  rim_proximity_cost)
+    COST(  1e12,   noncircular_rim_cost)
+    COST(  10e1,   nonequilateral_triangles_cost)
+//    COST(  1e1,   small_triangles_cost)
 //    COST(  1e6,   edge_angle_cost)
                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
-//    COST(  1e1,   small_triangles_cost)
-    COST(  1e12,   noncircular_rim_cost)
 #endif
 
 #if XBITS==5
 #define STOP_EPSILON 1.2e-4
-    COST(  3e7,   line_bending_cost)
-    COST( 10e2,   prop_edge_length_variation_cost)
-    COST( 9.0e3,  rim_proximity_cost) // 5e1 is too much
-                                     // 2.5e1 is too little
-    // 0.2e1 grows compared to previous ?
-    // 0.6e0 shrinks compared to previous ?
+    COST(  3e3,   vertex_displacement_cost)
+    COST(  3e3,   vertex_edgewise_displ_cost)
+    COST(  2e-1,  rim_proximity_cost)
+    COST(  1e12,   noncircular_rim_cost)
+    COST(  10e1,   nonequilateral_triangles_cost)
+//    COST(  1e1,   small_triangles_cost)
+//    COST(  1e6,   edge_angle_cost)
+                  #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
+#endif
 
-//    COST(  1e12,   edge_angle_cost)
-                  #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.3)
-    COST(  1e18,   noncircular_rim_cost)
+#if XBITS==6
+#define STOP_EPSILON 1.2e-4
+    COST(  3e3,   vertex_displacement_cost)
+    COST(  3e3,   vertex_edgewise_displ_cost)
+    COST( 0.02e0,  rim_proximity_cost)
+    COST(  1e12,   noncircular_rim_cost)
+    COST(  10e1,   nonequilateral_triangles_cost)
+//    COST(  1e1,   small_triangles_cost)
+//    COST(  1e6,   edge_angle_cost)
+                  #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
 #endif
 
-#if XBITS>=6 /* nonsense follows but never mind */
+#if XBITS>=7 /* nonsense follows but never mind */
 #define STOP_EPSILON 1e-6
     COST(  3e5,   line_bending_cost)
     COST( 10e2,   edge_length_variation_cost)
@@ -170,7 +182,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 x %g > %# e* |", tcost, tweight, tenergy);
+  if (pr) printf(" %# e > %# e* |", tcost, tenergy);
   *energy += tenergy;
 }
 
@@ -213,46 +225,117 @@ void compute_vertex_areas(const Vertices vertices, int section) {
   }
 }
 
-/*---------- combined vertex displacement ----------*/
+/*---------- displacement of vertices across a midpoint ----------*/
 
   /*
-   *   At vertex Q considering edge direction e to R
-   *   and corresponding opposite edge to P.
+   * Subroutine used where we have
+   *
+   *        R - - - - - - - M . -  -  -  -  R'
+   *                            ` .
+   *                                ` .
+   *                                    S
    *
-   *   Let  R' = Q + PQ
+   * and wish to say that the vector RM should be similar to MS
+   * or to put it another way S = M + RM
+   *
+   * NB this is not symmetric wrt R and S since it divides by
+   * |SM| but not |RM| so you probably want to call it twice.
+   *
+   * Details:
+   *
+   *   Let  R' = M + SM
    *        D  = R' - R
-   *                                                             delta
-   *              [       -1                                    ]
-   *   cost    =  [ lambda  . ( D . PQ/|PQ| ) + | D x PQ/|PQ| | ]
-   *       Q,e    [ ------------------------------------------- ]
-   *              [                      |PQ|                   ]
+   *
+   * Then the (1/delta)th power of the cost is
+   *      proportional to |D|, and
+   *      inversely proportional to |SM|
+   * except that |D| is measured in a wierd way which counts
+   * distance in the same direction as SM 1/lambda times as much
+   * ie the equipotential surfaces are ellipsoids around R',
+   * lengthened by lambda in the direction of RM.
+   *
+   * So
+   *                                                               delta
+   *                [       -1                                    ]
+   *   cost      =  [ lambda  . ( D . SM/|SM| ) + | D x SM/|SM| | ]
+   *       R,S,M    [ ------------------------------------------- ]
+   *                [                      |SM|                   ]
+   *
+   */
+
+static double vertex_one_displ_cost(const double r[D3], const double s[D3],
+                                   const double midpoint[D3],
+                                   double delta, double inv_lambda) {
+  const double smlen2_epsilon= 1e-12;
+  double sm[D3], d[D3], ddot, dcross[D3];
+  int k;
+
+  K sm[k]= -s[k] + midpoint[k];
+  K d[k]= midpoint[k] + sm[k] - r[k];
+  ddot= dotprod(d,sm);
+  xprod(dcross, d,sm);
+  double smlen2= magnD2(sm);
+  double cost_basis= inv_lambda * ddot + magnD(dcross);
+  double cost= pow(cost_basis / (smlen2 + smlen2_epsilon), delta);
+
+  return cost;
+}
+
+/*---------- displacement of vertices opposite at a vertex ----------*/
+
+  /*
+   *   At vertex Q considering edge direction e to R
+   *   and corresponding opposite edge to S.
+   *
+   *   This is vertex displacement as above with M=Q
    */
 
 double vertex_displacement_cost(const Vertices vertices, int section) {
   const double inv_lambda= 1.0/1; //2;
   const double delta= 4;
-  const double pqlen2_epsilon= 1e-12;
 
-  int pi,e,qi,ri, k;
-  double pq[D3], d[D3], ddot, dcross[D3];
+  int si,e,qi,ri;
   double total_cost= 0;
 
   FOR_EDGE(qi,e,ri, OUTER) {
-    pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
+    si= EDGE_END2(qi,(e+3)%V6); if (si<0) continue;
 
-    K pq[k]= -vertices[pi][k] + vertices[qi][k];
-    K d[k]= vertices[qi][k] + pq[k] - vertices[ri][k];
-    ddot= dotprod(d,pq);
-    xprod(dcross, d,pq);
-    double pqlen2= magnD2(pq);
-    double cost_basis= inv_lambda * ddot + magnD(dcross);
-    double cost= pow(cost_basis / (pqlen2 + pqlen2_epsilon), delta);
+    total_cost += vertex_one_displ_cost(vertices[ri], vertices[si], vertices[qi],
+                                       delta, inv_lambda);
+  }
+  return total_cost;
+}
 
-    total_cost += cost;
+/*---------- displacement of vertices opposite at an edge ----------*/
+
+  /*
+   *   At edge PQ considering vertices R and S (see diagram
+   *   below for overly sharp edge cost).
+   *
+   *   Let  M  = midpoint of PQ
+   */
+
+double vertex_edgewise_displ_cost(const Vertices vertices, int section) {
+  const double inv_lambda= 1.0/1; //2;
+  const double delta= 4;
+
+  int pi,e,qi,ri,si, k;
+  double m[D3];
+  double total_cost= 0;
+
+  FOR_EDGE(pi,e,qi, OUTER) {
+    si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
+    ri= EDGE_END2(pi,(e   +1)%V6);  if (ri<0) continue;
+
+    K m[k]= 0.5 * (vertices[pi][k] + vertices[qi][k]);
+    
+    total_cost += vertex_one_displ_cost(vertices[ri], vertices[si], m,
+                                       delta, inv_lambda);
   }
   return total_cost;
 }
 
+
 /*---------- at-vertex edge angles ----------*/
 
   /*
@@ -488,26 +571,9 @@ double edge_angle_cost(const Vertices vertices, int section) {
 /*---------- small triangles cost ----------*/
 
   /*
+   * Consider a triangle PQS
    *
-   *                       Q `-_
-   *              / |    `-_
-   *             /  |       `-.
-   *            /   |           S
-   *           /    |      _,-'
-   *          /     |  _,-'
-   *         /    , P '
-   *        /  ,-'
-   *       /,-'
-   *      /'
-   *            R
-   *
-   *  Let delta =  angle between two triangles' normals
-   *
-   *  Giving energy contribution:
-   *
-   *                                   2
-   *    E             =  F   .  delta
-   *     vd, edge PQ      vd
+   * Cost is 1/( area^2 )
    */
 
 double small_triangles_cost(const Vertices vertices, int section) {
@@ -538,3 +604,48 @@ double small_triangles_cost(const Vertices vertices, int section) {
 
   return total_cost;
 }
+
+/*---------- nonequilateral triangles cost ----------*/
+
+  /*
+   * Consider a triangle PQR
+   *
+   * let edge lengths a=|PQ| b=|QR| c=|RP|
+   *
+   * predicted edge length p = 1/3 * (a+b+c)
+   *
+   * compute cost for each x in {a,b,c}
+   *
+   *
+   *      cost      = (x-p)^2 / p^2
+   *          PQR,x
+   */
+
+double nonequilateral_triangles_cost(const Vertices vertices, int section) {
+  double pr[D3], abc[3];
+  int pi,e0,e1,qi,ri, k,i;
+  double our_epsilon=1e-6;
+  double total_cost= 0;
+
+  FOR_EDGE(pi,e0,qi, OUTER) {
+    e1= (e0+V6-1)%V6;
+    ri= EDGE_END2(pi,e1);  if (ri<0) continue;
+
+    K pr[k]= -vertices[pi][k] + vertices[ri][k];
+
+    abc[0]= edge_lengths[pi][e0]; /* PQ */
+    abc[1]= edge_lengths[qi][e1]; /* QR */
+    abc[2]= magnD(pr);
+
+    double p= (1/3.0) * (abc[0]+abc[1]+abc[2]);
+    double p_inv2= 1/(p*p + our_epsilon);
+
+    for (i=0; i<3; i++) {
+      double diff= (abc[i] - p);
+      double cost= diff*diff * p_inv2;
+      total_cost += cost;
+    }
+  }
+
+  return total_cost;
+}
index f75bc57..f4eddfb 100644 (file)
@@ -21,6 +21,7 @@ extern double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
 extern const double edge_angle_cost_circcircrat;
 
 double vertex_displacement_cost(const Vertices vertices, int section);
+double vertex_edgewise_displ_cost(const Vertices vertices, int section);
 double line_bending_cost(const Vertices vertices, int section);
 double noncircular_rim_cost(const Vertices vertices, int section);
 double edge_length_variation_cost(const Vertices vertices, int section);
@@ -28,6 +29,7 @@ double prop_edge_length_variation_cost(const Vertices vertices, int section);
 double rim_proximity_cost(const Vertices vertices, int section);
 double edge_angle_cost(const Vertices vertices, int section);
 double small_triangles_cost(const Vertices vertices, int section);
+double nonequilateral_triangles_cost(const Vertices vertices, int section);
 
 extern const char *input_file, *best_file;
 extern char *best_file_tmp;