wip new energy functions

diff --git a/common.c b/common.c
index 62b981efa930fd19bbb3eab3bb94dd0ba0d14d69..32701815aa753bcd1e43a4e71ee5d1cd1c43214b 100644 (file)
--- a/common.c
+++ b/common.c
@@ -15,6 +15,14 @@ double magnD(const double pq[D3]) {
   return gsl_blas_dnrm2(&v);
 }
 
+double magnD2(const double pq[D3]) {
+  int k;
+  double d2= 0;
+  
+  K d2= ffsqa(pq[k], d2);
+  return d2;
+}
+
 double hypotD(const double p[D3], const double q[D3]) {
   int k;
   double pq[D3];
@@ -43,6 +51,14 @@ void xprod(double r[D3], const double a[D3], const double b[D3]) {
   r[2]= a[0]*b[1] - a[1]*b[0];
 }
 
+void xprod_norm(double r[D3], const double a[D3], const double b[D3],
+               double absepsilon, double one) {
+  int k;
+  xprod(r,a,b);
+  double multby= one/(magnD(r) + absepsilon);
+  K r[k] *= multby;
+}
+
 double dotprod(const double a[D3], const double b[D3]) {
   int k;
   double result= 0;

diff --git a/common.h b/common.h
index 7e70a752666e8ad45ee485451bbb13b72d2cd6cb..a11f5660a6b7d4e8d88ef630c67b9521398a27e3 100644 (file)
--- a/common.h
+++ b/common.h
@@ -36,8 +36,11 @@ double hypotD2(const double p[D3], const double q[D3]);
 double hypotD2plus(const double p[D3], const double q[D3], double add);
 
 double magnD(const double pq[D3]);
+double magnD2(const double pq[D3]);
 void xprod(double r[D3], const double a[D3], const double b[D3]);
 double dotprod(const double a[D3], const double b[D3]);
+void xprod_norm(double r[D3], const double a[D3], const double b[D3],
+               double absepsilon, double one);
 
 void flushoutput(void);
 void diee(const char *what);

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;
+}

diff --git a/minimise.h b/minimise.h
index a0761250e28ae5131988c0c10f7bed321bb16a4f..60efb3e53a07d83985cb277e927d400897d4c988 100644 (file)
--- a/minimise.h
+++ b/minimise.h
@@ -21,6 +21,8 @@ double line_bending_cost(const Vertices vertices);
 double noncircular_rim_cost(const Vertices vertices);
 double edge_length_variation_cost(const Vertices vertices);
 double rim_proximity_cost(const Vertices vertices);
+double edge_angle_cost(const Vertices vertices);
+double small_triangles_cost(const Vertices vertices);
 
 extern const char *input_file, *best_file;
 extern char *best_file_tmp;