From: Ian Jackson <ian@davenant.greenend.org.uk>
Date: Thu, 21 Feb 2008 19:41:27 +0000 (+0000)
Subject: done circumcircle but it doesn't work so well
X-Git-Url: https://www.chiark.greenend.org.uk/ucgi/~ian/git?a=commitdiff_plain;h=2f8c87aa65ec594687eac88d4311475e289d8e51;p=moebius2.git

done circumcircle but it doesn't work so well
---

diff --git a/common.c b/common.c
index 3270181..cc72ed3 100644
--- a/common.c
+++ b/common.c
@@ -51,12 +51,16 @@ 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) {
+void normalise(double v[D3], double one, double absepsilon) {
   int k;
+  double multby= one/(magnD(v) + absepsilon);
+  K v[k] *= multby;
+}
+
+void xprod_norm(double r[D3], const double a[D3], const double b[D3],
+		double one, double absepsilon) {
   xprod(r,a,b);
-  double multby= one/(magnD(r) + absepsilon);
-  K r[k] *= multby;
+  normalise(r, absepsilon, one);
 }
 
 double dotprod(const double a[D3], const double b[D3]) {
diff --git a/common.h b/common.h
index a11f566..58bc1cc 100644
--- a/common.h
+++ b/common.h
@@ -40,7 +40,8 @@ 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);
+		double one, double absepsilon);
+void normalise(double v[D3], double one, double absepsilon);
 
 void flushoutput(void);
 void diee(const char *what);
diff --git a/energy.c b/energy.c
index cca777f..3459a72 100644
--- 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;
 }