energy.c

   1 /*
   2  * We try to find an optimal triangle grid
   3  */
   4
   5 #include "common.h"
   6 #include "minimise.h"
   7 #include "mgraph.h"
   8 #include "parallel.h"
   9
  10 double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6];
  11
  12 static double best_energy= DBL_MAX;
  13
  14 static void addcost(double *energy, double tweight, double tcost, int pr);
  15
  16 /*---------- main energy computation, weights, etc. ----------*/
  17
  18 typedef double CostComputation(const Vertices vertices, int section);
  19 typedef void PreComputation(const Vertices vertices, int section);
  20
  21 typedef struct {
  22   double weight;
  23   CostComputation *fn;
  24 } CostContribution;
  25
  26 #define NPRECOMPS ((sizeof(precomps)/sizeof(precomps[0])))
  27 #define NCOSTS ((sizeof(costs)/sizeof(costs[0])))
  28 #define COST(weight, compute) { (weight),(compute) },
  29
  30 static PreComputation *const precomps[]= {
  31   compute_edge_lengths,
  32   compute_vertex_areas
  33 };
  34
  35 static const CostContribution costs[]= {
  36
  37 #if XBITS==3
  38 #define STOP_EPSILON 1e-6
  39     COST(  3e3,   vertex_displacement_cost)
  40     COST( 0.4e3,  rim_proximity_cost)
  41     COST(  1e7,   edge_angle_cost)
  42                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.2/1.7)
  43     COST(  1e2,   small_triangles_cost)
  44     COST(  1e12,   noncircular_rim_cost)
  45 #endif
  46
  47 #if XBITS==4
  48 #define STOP_EPSILON 1.2e-3
  49     COST(  3e3,   vertex_displacement_cost)
  50     COST(  3e3,   vertex_edgewise_displ_cost)
  51     COST( 0.2e3,  rim_proximity_cost)
  52     COST(  1e12,   noncircular_rim_cost)
  53     COST(  10e1,   nonequilateral_triangles_cost)
  54 //    COST(  1e1,   small_triangles_cost)
  55 //    COST(  1e6,   edge_angle_cost)
  56                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
  57 #endif
  58
  59 #if XBITS==5
  60 #define STOP_EPSILON 1.2e-4
  61     COST(  3e3,   vertex_displacement_cost)
  62     COST(  3e3,   vertex_edgewise_displ_cost)
  63     COST(  2e-1,  rim_proximity_cost)
  64     COST(  1e12,   noncircular_rim_cost)
  65     COST(  10e1,   nonequilateral_triangles_cost)
  66 //    COST(  1e1,   small_triangles_cost)
  67 //    COST(  1e6,   edge_angle_cost)
  68                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
  69 #endif
  70
  71 #if XBITS==6
  72 #define STOP_EPSILON 1.2e-4
  73     COST(  3e3,   vertex_displacement_cost)
  74     COST(  3e3,   vertex_edgewise_displ_cost)
  75     COST( 0.02e0,  rim_proximity_cost)
  76     COST(  1e12,   noncircular_rim_cost)
  77     COST(  10e1,   nonequilateral_triangles_cost)
  78 //    COST(  1e1,   small_triangles_cost)
  79 //    COST(  1e6,   edge_angle_cost)
  80                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
  81 #endif
  82
  83 #if XBITS>=7 /* nonsense follows but never mind */
  84 #define STOP_EPSILON 1e-6
  85     COST(  3e5,   line_bending_cost)
  86     COST( 10e2,   edge_length_variation_cost)
  87     COST( 9.0e1,  rim_proximity_cost) // 5e1 is too much
  88                                       // 2.5e1 is too little
  89     // 0.2e1 grows compared to previous ?
  90     // 0.6e0 shrinks compared to previous ?
  91
  92     COST(  1e12,   edge_angle_cost)
  93                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.3)
  94     COST(  1e18,   noncircular_rim_cost)
  95 #endif
  96
  97 };
  98
  99 const double edge_angle_cost_circcircrat= EDGE_ANGLE_COST_CIRCCIRCRAT;
 100
 101 void energy_init(void) {
 102   stop_epsilon= STOP_EPSILON;
 103 }
 104
 105 /*---------- energy computation machinery ----------*/
 106
 107 void compute_energy_separately(const struct Vertices *vs,
 108                          int section, void *energies_v, void *totals_v) {
 109   double *energies= energies_v;
 110   int ci;
 111
 112   for (ci=0; ci<NPRECOMPS; ci++) {
 113     precomps[ci](vs->a, section);
 114     inparallel_barrier();
 115   }
 116   for (ci=0; ci<NCOSTS; ci++)
 117     energies[ci]= costs[ci].fn(vs->a, section);
 118 }
 119
 120 void compute_energy_combine(const struct Vertices *vertices,
 121                          int section, void *energies_v, void *totals_v) {
 122   int ci;
 123   double *energies= energies_v;
 124   double *totals= totals_v;
 125
 126   for (ci=0; ci<NCOSTS; ci++)
 127     totals[ci] += energies[ci];
 128 }
 129
 130 double compute_energy(const struct Vertices *vs) {
 131   static int bests_unprinted;
 132
 133   double totals[NCOSTS], energy;
 134   int ci, printing;
 135
 136   printing= printing_check(pr_cost,0);
 137
 138   if (printing) printf("%15lld c>e |", evaluations);
 139
 140   for (ci=0; ci<NCOSTS; ci++)
 141     totals[ci]= 0;
 142
 143   inparallel(vs,
 144              compute_energy_separately,
 145              compute_energy_combine,
 146              sizeof(totals) /* really, size of energies */,
 147              totals);
 148
 149   energy= 0;
 150   for (ci=0; ci<NCOSTS; ci++)
 151     addcost(&energy, costs[ci].weight, totals[ci], printing);
 152
 153   if (printing) printf("| total %# e |", energy);
 154
 155   if (energy < best_energy) {
 156     FILE *best_f;
 157     int r;
 158
 159     if (printing) {
 160       printf(" BEST");
 161       if (bests_unprinted) printf(" [%4d]",bests_unprinted);
 162       bests_unprinted= 0;
 163     } else {
 164       bests_unprinted++;
 165     }
 166
 167     best_f= fopen(best_file_tmp,"wb");  if (!best_f) diee("fopen new out");
 168     r= fwrite(vs->a,sizeof(vs->a),1,best_f); if (r!=1) diee("fwrite");
 169     if (fclose(best_f)) diee("fclose new best");
 170     if (rename(best_file_tmp,best_file)) diee("rename install new best");
 171
 172     best_energy= energy;
 173   }
 174   if (printing) {
 175     putchar('\n');
 176     flushoutput();
 177   }
 178
 179   evaluations++;
 180   return energy;
 181 }
 182
 183 static void addcost(double *energy, double tweight, double tcost, int pr) {
 184   double tenergy= tweight * tcost;
 185   if (pr) printf(" %# e > %# e* |", tcost, tenergy);
 186   *energy += tenergy;
 187 }
 188
 189 /*---------- Precomputations ----------*/
 190
 191 void compute_edge_lengths(const Vertices vertices, int section) {
 192   int v1,e,v2;
 193
 194   FOR_EDGE(v1,e,v2, OUTER)
 195     edge_lengths[v1][e]= hypotD(vertices[v1],vertices[v2]);
 196 }
 197
 198 void compute_vertex_areas(const Vertices vertices, int section) {
 199   int v0,v1,v2, e1,e2;
 200 //  int k;
 201
 202   FOR_VERTEX(v0, OUTER) {
 203     double total= 0.0, edges_total=0;
 204     int count= 0;
 205
 206     FOR_VEDGE(v0,e1,v1) {
 207       e2= (e1+1) % V6;
 208       v2= EDGE_END2(v0,e2);
 209       if (v2<0) continue;
 210
 211       edges_total += edge_lengths[v0][e1];
 212
 213 //      double e1v[D3], e2v[D3], av[D3];
 214 //      K {
 215 //      e1v[k]= vertices[v1][k] - vertices[v0][k];
 216 //      e2v[k]= vertices[v2][k] - vertices[v0][k];
 217 //      }
 218 //      xprod(av, e1v, e2v);
 219 //      total += magnD(av);
 220
 221       count++;
 222     }
 223     vertex_areas[v0]= total / count;
 224     vertex_mean_edge_lengths[v0]= edges_total / count;
 225   }
 226 }
 227
 228 /*---------- displacement of vertices across a midpoint ----------*/
 229
 230   /*
 231    * Subroutine used where we have
 232    *
 233    *        R - - - - - - - M . -  -  -  -  R'
 234    *                            ` .
 235    *                                ` .
 236    *                                    S
 237    *
 238    * and wish to say that the vector RM should be similar to MS
 239    * or to put it another way S = M + RM
 240    *
 241    * NB this is not symmetric wrt R and S since it divides by
 242    * |SM| but not |RM| so you probably want to call it twice.
 243    *
 244    * Details:
 245    *
 246    *   Let  R' = M + SM
 247    *        D  = R' - R
 248    *
 249    * Then the (1/delta)th power of the cost is
 250    *      proportional to |D|, and
 251    *      inversely proportional to |SM|
 252    * except that |D| is measured in a wierd way which counts
 253    * distance in the same direction as SM 1/lambda times as much
 254    * ie the equipotential surfaces are ellipsoids around R',
 255    * lengthened by lambda in the direction of RM.
 256    *
 257    * So
 258    *                                                               delta
 259    *                [       -1                                    ]
 260    *   cost      =  [ lambda  . ( D . SM/|SM| ) + | D x SM/|SM| | ]
 261    *       R,S,M    [ ------------------------------------------- ]
 262    *                [                      |SM|                   ]
 263    *
 264    */
 265
 266 static double vertex_one_displ_cost(const double r[D3], const double s[D3],
 267                                     const double midpoint[D3],
 268                                     double delta, double inv_lambda) {
 269   const double smlen2_epsilon= 1e-12;
 270   double sm[D3], d[D3], ddot, dcross[D3];
 271   int k;
 272
 273   K sm[k]= -s[k] + midpoint[k];
 274   K d[k]= midpoint[k] + sm[k] - r[k];
 275   ddot= dotprod(d,sm);
 276   xprod(dcross, d,sm);
 277   double smlen2= magnD2(sm);
 278   double cost_basis= inv_lambda * ddot + magnD(dcross);
 279   double cost= pow(cost_basis / (smlen2 + smlen2_epsilon), delta);
 280
 281   return cost;
 282 }
 283
 284 /*---------- displacement of vertices opposite at a vertex ----------*/
 285
 286   /*
 287    *   At vertex Q considering edge direction e to R
 288    *   and corresponding opposite edge to S.
 289    *
 290    *   This is vertex displacement as above with M=Q
 291    */
 292
 293 double vertex_displacement_cost(const Vertices vertices, int section) {
 294   const double inv_lambda= 1.0/1; //2;
 295   const double delta= 4;
 296
 297   int si,e,qi,ri;
 298   double total_cost= 0;
 299
 300   FOR_EDGE(qi,e,ri, OUTER) {
 301     si= EDGE_END2(qi,(e+3)%V6); if (si<0) continue;
 302
 303     total_cost += vertex_one_displ_cost(vertices[ri], vertices[si], vertices[qi],
 304                                         delta, inv_lambda);
 305   }
 306   return total_cost;
 307 }
 308
 309 /*---------- displacement of vertices opposite at an edge ----------*/
 310
 311   /*
 312    *   At edge PQ considering vertices R and S (see diagram
 313    *   below for overly sharp edge cost).
 314    *
 315    *   Let  M  = midpoint of PQ
 316    */
 317
 318 double vertex_edgewise_displ_cost(const Vertices vertices, int section) {
 319   const double inv_lambda= 1.0/1; //2;
 320   const double delta= 4;
 321
 322   int pi,e,qi,ri,si, k;
 323   double m[D3];
 324   double total_cost= 0;
 325
 326   FOR_EDGE(pi,e,qi, OUTER) {
 327     si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
 328     ri= EDGE_END2(pi,(e   +1)%V6);  if (ri<0) continue;
 329
 330     K m[k]= 0.5 * (vertices[pi][k] + vertices[qi][k]);
 331
 332     total_cost += vertex_one_displ_cost(vertices[ri], vertices[si], m,
 333                                         delta, inv_lambda);
 334   }
 335   return total_cost;
 336 }
 337
 338
 339 /*---------- at-vertex edge angles ----------*/
 340
 341   /*
 342    * Definition:
 343    *
 344    *    At each vertex Q, in each direction e:
 345    *
 346    *                                  e
 347    *                           Q ----->----- R
 348    *                      _,-'\__/
 349    *                  _,-'       delta
 350    *               P '
 351    *
 352    *                      r
 353    *       cost    = delta          (we use r=3)
 354    *           Q,e
 355    *
 356    *
 357    * Calculation:
 358    *
 359    *      Let vector A = PQ
 360    *                 B = QR
 361    *
 362    *                   -1   A . B
 363    *      delta =  tan     -------
 364    *                      | A x B |
 365    *
 366    *      which is always in the range 0..pi because the denominator
 367    *      is nonnegative.  We add epsilon to |AxB| to avoid division
 368    *      by zero.
 369    *
 370    *                     r
 371    *      cost    = delta
 372    *          Q,e
 373    */
 374
 375 double line_bending_cost(const Vertices vertices, int section) {
 376   static const double axb_epsilon= 1e-6;
 377   static const double exponent_r= 4;
 378
 379   int pi,e,qi,ri, k;
 380   double  a[D3], b[D3], axb[D3];
 381   double total_cost= 0;
 382
 383   FOR_EDGE(qi,e,ri, OUTER) {
 384     pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
 385
 386 //if (!(qi&XMASK)) fprintf(stderr,"%02x-%02x-%02x (%d)\n",pi,qi,ri,e);
 387
 388     K a[k]= -vertices[pi][k] + vertices[qi][k];
 389     K b[k]= -vertices[qi][k] + vertices[ri][k];
 390
 391     xprod(axb,a,b);
 392
 393     double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b));
 394     double cost= pow(delta,exponent_r);
 395
 396     total_cost += cost;
 397   }
 398   return total_cost;
 399 }
 400
 401 /*---------- edge length variation ----------*/
 402
 403   /*
 404    * Definition:
 405    *
 406    *    See the diagram above.
 407    *                                r
 408    *       cost    = ( |PQ| - |QR| )
 409    *           Q,e
 410    */
 411
 412 double edge_length_variation_cost(const Vertices vertices, int section) {
 413   double diff, cost= 0, exponent_r= 2;
 414   int q, e,r, eback;
 415
 416   FOR_EDGE(q,e,r, OUTER) {
 417     eback= edge_reverse(q,e);
 418     diff= edge_lengths[q][e] - edge_lengths[q][eback];
 419     cost += pow(diff,exponent_r);
 420   }
 421   return cost;
 422 }
 423
 424 /*---------- proportional edge length variation ----------*/
 425
 426   /*
 427    * Definition:
 428    *
 429    *    See the diagram above.
 430    *                                r
 431    *       cost    = ( |PQ| - |QR| )
 432    *           Q,e
 433    */
 434
 435 double prop_edge_length_variation_cost(const Vertices vertices, int section) {
 436   const double num_epsilon= 1e-6;
 437
 438   double cost= 0, exponent_r= 2;
 439   int q, e,r, eback;
 440
 441   FOR_EDGE(q,e,r, OUTER) {
 442     eback= edge_reverse(q,e);
 443     double le= edge_lengths[q][e];
 444     double leback= edge_lengths[q][eback];
 445     double diff= le - leback;
 446     double num= MIN(le, leback);
 447     cost += pow(diff / (num + num_epsilon), exponent_r);
 448   }
 449   return cost;
 450 }
 451
 452 /*---------- rim proximity cost ----------*/
 453
 454 static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) {
 455   /* By symmetry, nearest point on circle is the one with
 456    * the same angle subtended at the z axis. */
 457   oncircle[0]= p[0];
 458   oncircle[1]= p[1];
 459   oncircle[2]= 0;
 460   double mult= 1.0/ magnD(oncircle);
 461   oncircle[0] *= mult;
 462   oncircle[1] *= mult;
 463 }
 464
 465 double rim_proximity_cost(const Vertices vertices, int section) {
 466   double oncircle[3], cost=0;
 467   int v;
 468
 469   FOR_VERTEX(v, OUTER) {
 470     int y= v >> YSHIFT;
 471     int nominal_edge_distance= y <= Y/2 ? y : Y-1-y;
 472     if (nominal_edge_distance==0) continue;
 473
 474     find_nearest_oncircle(oncircle, vertices[v]);
 475
 476     cost +=
 477       vertex_mean_edge_lengths[v] *
 478       (nominal_edge_distance*nominal_edge_distance) /
 479       (hypotD2(vertices[v], oncircle) + 1e-6);
 480   }
 481   return cost;
 482 }
 483
 484 /*---------- noncircular rim cost ----------*/
 485
 486 double noncircular_rim_cost(const Vertices vertices, int section) {
 487   int vy,vx,v;
 488   double cost= 0.0;
 489   double oncircle[3];
 490
 491   FOR_RIM_VERTEX(vy,vx,v, OUTER) {
 492     find_nearest_oncircle(oncircle, vertices[v]);
 493
 494     double d2= hypotD2(vertices[v], oncircle);
 495     cost += d2*d2;
 496   }
 497   return cost;
 498 }
 499
 500 /*---------- overly sharp edge cost ----------*/
 501
 502   /*
 503    *
 504    *                Q `-_
 505    *              / |    `-_                           P'Q' ------ S'
 506    *             /  |       `-.                      _,' `.       .
 507    *            /   |           S                 _,'     :      .
 508    *           /    |      _,-'                _,'        :r    .r
 509    *          /     |  _,-'                R' '           `.   .
 510    *         /    , P '                        ` .   r     :  .
 511    *        /  ,-'                                  `  .   :
 512    *       /,-'                                          ` C'
 513    *      /'
 514    *     R
 515    *
 516    *
 517    *
 518    *  Let delta =  angle between two triangles' normals
 519    *
 520    *  Giving energy contribution:
 521    *
 522    *                                   2
 523    *    E             =  F   .  delta
 524    *     vd, edge PQ      vd
 525    */
 526
 527 double edge_angle_cost(const Vertices vertices, int section) {
 528   double pq1[D3], rp[D3], ps[D3], rp_2d[D3], ps_2d[D3], rs_2d[D3];
 529   double a,b,c,s,r;
 530   const double minradius_base= 0.2;
 531
 532   int pi,e,qi,ri,si, k;
 533 //  double our_epsilon=1e-6;
 534   double total_cost= 0;
 535
 536   FOR_EDGE(pi,e,qi, OUTER) {
 537 //    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
 538
 539     si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
 540     ri= EDGE_END2(pi,(e   +1)%V6);  if (ri<0) continue;
 541
 542     K {
 543       pq1[k]= -vertices[pi][k] + vertices[qi][k];
 544       rp[k]=  -vertices[ri][k] + vertices[pi][k];
 545       ps[k]=  -vertices[pi][k] + vertices[si][k];
 546     }
 547
 548     normalise(pq1,1,1e-6);
 549     xprod(rp_2d, rp,pq1); /* projects RP into plane normal to PQ */
 550     xprod(ps_2d, ps,pq1); /* likewise PS */
 551     K rs_2d[k]= rp_2d[k] + ps_2d[k];
 552     /* radius of circumcircle of R'P'S' from Wikipedia
 553      * `Circumscribed circle' */
 554     a= magnD(rp_2d);
 555     b= magnD(ps_2d);
 556     c= magnD(rs_2d);
 557     s= 0.5*(a+b+c);
 558     r= a*b*c / sqrt((a+b+c)*(a-b+c)*(b-c+a)*(c-a+b) + 1e-6);
 559
 560     double minradius= minradius_base + edge_angle_cost_circcircrat*(a+b);
 561     double deficit= minradius - r;
 562     if (deficit < 0) continue;
 563     double cost= deficit*deficit;
 564
 565     total_cost += cost;
 566   }
 567
 568   return total_cost;
 569 }
 570
 571 /*---------- small triangles cost ----------*/
 572
 573   /*
 574    * Consider a triangle PQS
 575    *
 576    * Cost is 1/( area^2 )
 577    */
 578
 579 double small_triangles_cost(const Vertices vertices, int section) {
 580   double pq[D3], ps[D3];
 581   double x[D3];
 582   int pi,e,qi,si, k;
 583 //  double our_epsilon=1e-6;
 584   double total_cost= 0;
 585
 586   FOR_EDGE(pi,e,qi, OUTER) {
 587 //    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
 588
 589     si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
 590
 591     K {
 592       pq[k]= vertices[qi][k] - vertices[pi][k];
 593       ps[k]= vertices[si][k] - vertices[pi][k];
 594     }
 595     xprod(x, pq,ps);
 596
 597     double cost= 1/(magnD2(x) + 0.01);
 598
 599 //double cost= pow(magnD(spqxpqr), 3);
 600 //assert(dot>=-1 && dot <=1);
 601 //double cost= 1-dot;
 602     total_cost += cost;
 603   }
 604
 605   return total_cost;
 606 }
 607
 608 /*---------- nonequilateral triangles cost ----------*/
 609
 610   /*
 611    * Consider a triangle PQR
 612    *
 613    * let edge lengths a=|PQ| b=|QR| c=|RP|
 614    *
 615    * predicted edge length p = 1/3 * (a+b+c)
 616    *
 617    * compute cost for each x in {a,b,c}
 618    *
 619    *
 620    *      cost      = (x-p)^2 / p^2
 621    *          PQR,x
 622    */
 623
 624 double nonequilateral_triangles_cost(const Vertices vertices, int section) {
 625   double pr[D3], abc[3];
 626   int pi,e0,e1,qi,ri, k,i;
 627   double our_epsilon=1e-6;
 628   double total_cost= 0;
 629
 630   FOR_EDGE(pi,e0,qi, OUTER) {
 631     e1= (e0+V6-1)%V6;
 632     ri= EDGE_END2(pi,e1);  if (ri<0) continue;
 633
 634     K pr[k]= -vertices[pi][k] + vertices[ri][k];
 635
 636     abc[0]= edge_lengths[pi][e0]; /* PQ */
 637     abc[1]= edge_lengths[qi][e1]; /* QR */
 638     abc[2]= magnD(pr);
 639
 640     double p= (1/3.0) * (abc[0]+abc[1]+abc[2]);
 641     double p_inv2= 1/(p*p + our_epsilon);
 642
 643     for (i=0; i<3; i++) {
 644       double diff= (abc[i] - p);
 645       double cost= diff*diff * p_inv2;
 646       total_cost += cost;
 647     }
 648   }
 649
 650   return total_cost;
 651 }