wip new energy functions

[moebius2.git] / energy.c

diff --git a/energy.c b/energy.c
index 80ac7e885f989e627b2c01edf3f0ba02a6bf5c73..cca777f8cff22ded6570402e06cfbb3fe2fe5b59 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -35,16 +35,19 @@ double compute_energy(const struct Vertices *vs) {
   if (XBITS==3) {
     COST(  3e2,   line_bending_cost(vs->a));
     COST(  1e3,   edge_length_variation_cost(vs->a));
-    COST( 0.2e3,  rim_proximity_cost(vs->a));
-    COST(  1e8,   noncircular_rim_cost(vs->a));
+    COST( 0.4e3,  rim_proximity_cost(vs->a));
+    COST(  1e6,   edge_angle_cost(vs->a));
+//    COST(  1e1,   small_triangles_cost(vs->a));
+    COST(  1e12,   noncircular_rim_cost(vs->a));
     stop_epsilon= 1e-6;
   } else if (XBITS==4) {
     COST(  3e2,   line_bending_cost(vs->a));
     COST(  3e3,   edge_length_variation_cost(vs->a));
-    COST( 3.8e1,  rim_proximity_cost(vs->a)); // 5e1 is too much
+    COST( 9.0e1,  rim_proximity_cost(vs->a)); // 5e1 is too much
                                                 // 2.5e1 is too little
     // 0.2e1 grows compared to previous ?
     // 0.6e0 shrinks compared to previous ?
+    COST(  1e12,   edge_angle_cost(vs->a));
     COST(  1e12,   noncircular_rim_cost(vs->a));
     stop_epsilon= 1e-5;
   } else {
@@ -259,3 +262,123 @@ double noncircular_rim_cost(const Vertices vertices) {
   }
   return cost;
 }
+
+/*---------- triangle bad normals cost ----------*/
+
+  /*
+   *
+   *                       Q `-_
+   *              / |    `-_
+   *             /  |       `-.
+   *            /   |           S
+   *           /    |      _,-'
+   *          /     |  _,-'
+   *         /    , P '
+   *        /  ,-'
+   *       /,-'
+   *      /'
+   *            R
+   *
+   *  Let delta =  angle between two triangles' normals
+   *
+   *  Giving energy contribution:
+   *
+   *                                   2
+   *    E             =  F   .  delta
+   *     vd, edge PQ      vd
+   */
+
+double edge_angle_cost(const Vertices vertices) {
+  double sq[D3], pq[D3], qr[D3], sqp[D3], pqr[D3], rs[D3];
+  double x[D3], y[D3];
+  int pi,e,qi,ri,si, k;
+//  double our_epsilon=1e-6;
+  double total_cost= 0;
+  
+  FOR_EDGE(pi,e,qi) {
+//    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
+
+    si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
+    ri= EDGE_END2(pi,(e   +1)%V6);  if (ri<0) continue;
+
+    K {
+      sq[k]= vertices[si][k] - vertices[qi][k];
+      pq[k]= vertices[pi][k] - vertices[qi][k];
+      qr[k]= vertices[qi][k] - vertices[ri][k];
+    }
+    xprod(sqp, pq,sq);
+    xprod(pqr, pq,qr);
+
+    double dot= dotprod(sqp,pqr);
+    xprod(x,sqp,pqr);
+
+    K rs[k]= -vertices[ri][k] + vertices[si][k];
+    xprod(y, pq,rs);
+
+    double delta=
+      pow(atan2(magnD(x), dot), 2.0) * pow(magnD2(pq), 2.0) /
+      (pow(magnD(y), 0.3) + 1e-6);
+    double cost= pow(delta, 2.0);
+
+//double cost= pow(magnD(spqxpqr), 3);
+//assert(dot>=-1 && dot <=1);
+//double cost= 1-dot;
+    total_cost += cost;
+  }
+  
+  return total_cost;
+}
+
+/*---------- small triangles cost ----------*/
+
+  /*
+   *
+   *                       Q `-_
+   *              / |    `-_
+   *             /  |       `-.
+   *            /   |           S
+   *           /    |      _,-'
+   *          /     |  _,-'
+   *         /    , P '
+   *        /  ,-'
+   *       /,-'
+   *      /'
+   *            R
+   *
+   *  Let delta =  angle between two triangles' normals
+   *
+   *  Giving energy contribution:
+   *
+   *                                   2
+   *    E             =  F   .  delta
+   *     vd, edge PQ      vd
+   */
+
+double small_triangles_cost(const Vertices vertices) {
+  double pq[D3], ps[D3];
+  double x[D3];
+  int pi,e,qi,si, k;
+//  double our_epsilon=1e-6;
+  double total_cost= 0;
+  
+  FOR_EDGE(pi,e,qi) {
+//    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
+
+    si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
+
+    K {
+      pq[k]= vertices[qi][k] - vertices[pi][k];
+      ps[k]= vertices[si][k] - vertices[pi][k];
+    }
+    xprod(x, pq,ps);
+
+    double cost= 1/(magnD2(x) + 0.01);
+
+//double cost= pow(magnD(spqxpqr), 3);
+//assert(dot>=-1 && dot <=1);
+//double cost= 1-dot;
+    total_cost += cost;
+  }
+  
+  return total_cost;
+}