done circumcircle but it doesn't work so well

[moebius2.git] / energy.c

diff --git a/energy.c b/energy.c
index cca777f8cff22ded6570402e06cfbb3fe2fe5b59..3459a726351724226f4fe7994c3fdc75d5b4d909 100644 (file)
--- a/energy.c
+++ b/energy.c
@@ -20,7 +20,7 @@ void energy_init(void) {
 
 double compute_energy(const struct Vertices *vs) {
   static int bests_unprinted;
-  
+
   double energy;
   int printing;
 
@@ -263,22 +263,24 @@ double noncircular_rim_cost(const Vertices vertices) {
   return cost;
 }
 
-/*---------- triangle bad normals cost ----------*/
+/*---------- overly sharp edge cost ----------*/
 
   /*
    *
    *                       Q `-_
-   *              / |    `-_
-   *             /  |       `-.
-   *            /   |           S
-   *           /    |      _,-'
-   *          /     |  _,-'
-   *         /    , P '
-   *        /  ,-'
-   *       /,-'
+   *              / |    `-_                          P'Q' ------ S'
+   *             /  |       `-.                             _,' `.       .
+   *            /   |           S                _,'     :      .
+   *           /    |      _,-'                _,'        :r    .r
+   *          /     |  _,-'               R' '           `.   .
+   *         /    , P '                               ` .   r     :  .
+   *        /  ,-'                                 `  .   : 
+   *       /,-'                                                 ` C'
    *      /'
    *            R
    *
+   *
+   *
    *  Let delta =  angle between two triangles' normals
    *
    *  Giving energy contribution:
@@ -289,12 +291,14 @@ double noncircular_rim_cost(const Vertices vertices) {
    */
 
 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];
+  double pq1[D3], rp[D3], ps[D3], rp_2d[D3], ps_2d[D3], rs_2d[D3];
+  double a,b,c,s,r;
+  const double minradius= 0.5;
+
   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;
 
@@ -302,30 +306,31 @@ double edge_angle_cost(const Vertices vertices) {
     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];
+      pq1[k]= -vertices[pi][k] + vertices[qi][k];
+      rp[k]=  -vertices[ri][k] + vertices[pi][k];
+      ps[k]=  -vertices[pi][k] + vertices[si][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);
+    normalise(pq1,1,1e-6);
+    xprod(rp_2d, rp,pq1); /* projects RP into plane normal to PQ */
+    xprod(ps_2d, ps,pq1); /* likewise PS */
+    K rs_2d[k]= rp_2d[k] + ps_2d[k];
+    /* radius of circumcircle of R'P'S' from Wikipedia
+     * `Circumscribed circle' */
+    a= magnD(rp_2d);
+    b= magnD(ps_2d);
+    c= magnD(rs_2d);
+    s= 0.5*(a+b+c);
+    r= a*b*c / sqrt((a+b+c)*(a-b+c)*(b-c+a)*(c-a+b) + 1e-6);
+//fprintf(stderr,"triangle a=%g b=%g c=%g -> s=%g r=%g\n",a,b,c,s,r);
+
+    double deficit= minradius - r;
+    if (deficit < 0) continue;
+    double cost= deficit*deficit;
 
-    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;
 }
 
@@ -360,7 +365,7 @@ double small_triangles_cost(const Vertices vertices) {
   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;
 
@@ -379,6 +384,6 @@ double small_triangles_cost(const Vertices vertices) {
 //double cost= 1-dot;
     total_cost += cost;
   }
-  
+
   return total_cost;
 }