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-4
  49     COST(  3e3,   vertex_displacement_cost)
  50     COST( 0.2e3,  rim_proximity_cost)
  51 //    COST(  1e6,   edge_angle_cost)
  52                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7)
  53 //    COST(  1e1,   small_triangles_cost)
  54     COST(  1e12,   noncircular_rim_cost)
  55 #endif
  56
  57 #if XBITS==5
  58 #define STOP_EPSILON 1.2e-4
  59     COST(  3e7,   line_bending_cost)
  60     COST( 10e2,   prop_edge_length_variation_cost)
  61     COST( 9.0e3,  rim_proximity_cost) // 5e1 is too much
  62                                       // 2.5e1 is too little
  63     // 0.2e1 grows compared to previous ?
  64     // 0.6e0 shrinks compared to previous ?
  65
  66 //    COST(  1e12,   edge_angle_cost)
  67                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.3)
  68     COST(  1e18,   noncircular_rim_cost)
  69 #endif
  70
  71 #if XBITS>=6 /* nonsense follows but never mind */
  72 #define STOP_EPSILON 1e-6
  73     COST(  3e5,   line_bending_cost)
  74     COST( 10e2,   edge_length_variation_cost)
  75     COST( 9.0e1,  rim_proximity_cost) // 5e1 is too much
  76                                       // 2.5e1 is too little
  77     // 0.2e1 grows compared to previous ?
  78     // 0.6e0 shrinks compared to previous ?
  79
  80     COST(  1e12,   edge_angle_cost)
  81                   #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.3)
  82     COST(  1e18,   noncircular_rim_cost)
  83 #endif
  84
  85 };
  86
  87 const double edge_angle_cost_circcircrat= EDGE_ANGLE_COST_CIRCCIRCRAT;
  88
  89 void energy_init(void) {
  90   stop_epsilon= STOP_EPSILON;
  91 }
  92
  93 /*---------- energy computation machinery ----------*/
  94
  95 void compute_energy_separately(const struct Vertices *vs,
  96                          int section, void *energies_v, void *totals_v) {
  97   double *energies= energies_v;
  98   int ci;
  99
 100   for (ci=0; ci<NPRECOMPS; ci++) {
 101     precomps[ci](vs->a, section);
 102     inparallel_barrier();
 103   }
 104   for (ci=0; ci<NCOSTS; ci++)
 105     energies[ci]= costs[ci].fn(vs->a, section);
 106 }
 107
 108 void compute_energy_combine(const struct Vertices *vertices,
 109                          int section, void *energies_v, void *totals_v) {
 110   int ci;
 111   double *energies= energies_v;
 112   double *totals= totals_v;
 113
 114   for (ci=0; ci<NCOSTS; ci++)
 115     totals[ci] += energies[ci];
 116 }
 117
 118 double compute_energy(const struct Vertices *vs) {
 119   static int bests_unprinted;
 120
 121   double totals[NCOSTS], energy;
 122   int ci, printing;
 123
 124   printing= printing_check(pr_cost,0);
 125
 126   if (printing) printf("%15lld c>e |", evaluations);
 127
 128   for (ci=0; ci<NCOSTS; ci++)
 129     totals[ci]= 0;
 130
 131   inparallel(vs,
 132              compute_energy_separately,
 133              compute_energy_combine,
 134              sizeof(totals) /* really, size of energies */,
 135              totals);
 136
 137   energy= 0;
 138   for (ci=0; ci<NCOSTS; ci++)
 139     addcost(&energy, costs[ci].weight, totals[ci], printing);
 140
 141   if (printing) printf("| total %# e |", energy);
 142
 143   if (energy < best_energy) {
 144     FILE *best_f;
 145     int r;
 146
 147     if (printing) {
 148       printf(" BEST");
 149       if (bests_unprinted) printf(" [%4d]",bests_unprinted);
 150       bests_unprinted= 0;
 151     } else {
 152       bests_unprinted++;
 153     }
 154
 155     best_f= fopen(best_file_tmp,"wb");  if (!best_f) diee("fopen new out");
 156     r= fwrite(vs->a,sizeof(vs->a),1,best_f); if (r!=1) diee("fwrite");
 157     if (fclose(best_f)) diee("fclose new best");
 158     if (rename(best_file_tmp,best_file)) diee("rename install new best");
 159
 160     best_energy= energy;
 161   }
 162   if (printing) {
 163     putchar('\n');
 164     flushoutput();
 165   }
 166
 167   evaluations++;
 168   return energy;
 169 }
 170
 171 static void addcost(double *energy, double tweight, double tcost, int pr) {
 172   double tenergy= tweight * tcost;
 173   if (pr) printf(" %# e x %g > %# e* |", tcost, tweight, tenergy);
 174   *energy += tenergy;
 175 }
 176
 177 /*---------- Precomputations ----------*/
 178
 179 void compute_edge_lengths(const Vertices vertices, int section) {
 180   int v1,e,v2;
 181
 182   FOR_EDGE(v1,e,v2, OUTER)
 183     edge_lengths[v1][e]= hypotD(vertices[v1],vertices[v2]);
 184 }
 185
 186 void compute_vertex_areas(const Vertices vertices, int section) {
 187   int v0,v1,v2, e1,e2;
 188 //  int k;
 189
 190   FOR_VERTEX(v0, OUTER) {
 191     double total= 0.0, edges_total=0;
 192     int count= 0;
 193
 194     FOR_VEDGE(v0,e1,v1) {
 195       e2= (e1+1) % V6;
 196       v2= EDGE_END2(v0,e2);
 197       if (v2<0) continue;
 198
 199       edges_total += edge_lengths[v0][e1];
 200
 201 //      double e1v[D3], e2v[D3], av[D3];
 202 //      K {
 203 //      e1v[k]= vertices[v1][k] - vertices[v0][k];
 204 //      e2v[k]= vertices[v2][k] - vertices[v0][k];
 205 //      }
 206 //      xprod(av, e1v, e2v);
 207 //      total += magnD(av);
 208
 209       count++;
 210     }
 211     vertex_areas[v0]= total / count;
 212     vertex_mean_edge_lengths[v0]= edges_total / count;
 213   }
 214 }
 215
 216 /*---------- combined vertex displacement ----------*/
 217
 218   /*
 219    *   At vertex Q considering edge direction e to R
 220    *   and corresponding opposite edge to P.
 221    *
 222    *   Let  R' = Q + PQ
 223    *        D  = R' - R
 224    *                                                             delta
 225    *              [       -1                                    ]
 226    *   cost    =  [ lambda  . ( D . PQ/|PQ| ) + | D x PQ/|PQ| | ]
 227    *       Q,e    [ ------------------------------------------- ]
 228    *              [                      |PQ|                   ]
 229    */
 230
 231 double vertex_displacement_cost(const Vertices vertices, int section) {
 232   const double inv_lambda= 1.0/1; //2;
 233   const double delta= 4;
 234   const double pqlen2_epsilon= 1e-12;
 235
 236   int pi,e,qi,ri, k;
 237   double pq[D3], d[D3], ddot, dcross[D3];
 238   double total_cost= 0;
 239
 240   FOR_EDGE(qi,e,ri, OUTER) {
 241     pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
 242
 243     K pq[k]= -vertices[pi][k] + vertices[qi][k];
 244     K d[k]= vertices[qi][k] + pq[k] - vertices[ri][k];
 245     ddot= dotprod(d,pq);
 246     xprod(dcross, d,pq);
 247     double pqlen2= magnD2(pq);
 248     double cost_basis= inv_lambda * ddot + magnD(dcross);
 249     double cost= pow(cost_basis / (pqlen2 + pqlen2_epsilon), delta);
 250
 251     total_cost += cost;
 252   }
 253   return total_cost;
 254 }
 255
 256 /*---------- at-vertex edge angles ----------*/
 257
 258   /*
 259    * Definition:
 260    *
 261    *    At each vertex Q, in each direction e:
 262    *
 263    *                                  e
 264    *                           Q ----->----- R
 265    *                      _,-'\__/
 266    *                  _,-'       delta
 267    *               P '
 268    *
 269    *                      r
 270    *       cost    = delta          (we use r=3)
 271    *           Q,e
 272    *
 273    *
 274    * Calculation:
 275    *
 276    *      Let vector A = PQ
 277    *                 B = QR
 278    *
 279    *                   -1   A . B
 280    *      delta =  tan     -------
 281    *                      | A x B |
 282    *
 283    *      which is always in the range 0..pi because the denominator
 284    *      is nonnegative.  We add epsilon to |AxB| to avoid division
 285    *      by zero.
 286    *
 287    *                     r
 288    *      cost    = delta
 289    *          Q,e
 290    */
 291
 292 double line_bending_cost(const Vertices vertices, int section) {
 293   static const double axb_epsilon= 1e-6;
 294   static const double exponent_r= 4;
 295
 296   int pi,e,qi,ri, k;
 297   double  a[D3], b[D3], axb[D3];
 298   double total_cost= 0;
 299
 300   FOR_EDGE(qi,e,ri, OUTER) {
 301     pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
 302
 303 //if (!(qi&XMASK)) fprintf(stderr,"%02x-%02x-%02x (%d)\n",pi,qi,ri,e);
 304
 305     K a[k]= -vertices[pi][k] + vertices[qi][k];
 306     K b[k]= -vertices[qi][k] + vertices[ri][k];
 307
 308     xprod(axb,a,b);
 309
 310     double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b));
 311     double cost= pow(delta,exponent_r);
 312
 313     total_cost += cost;
 314   }
 315   return total_cost;
 316 }
 317
 318 /*---------- edge length variation ----------*/
 319
 320   /*
 321    * Definition:
 322    *
 323    *    See the diagram above.
 324    *                                r
 325    *       cost    = ( |PQ| - |QR| )
 326    *           Q,e
 327    */
 328
 329 double edge_length_variation_cost(const Vertices vertices, int section) {
 330   double diff, cost= 0, exponent_r= 2;
 331   int q, e,r, eback;
 332
 333   FOR_EDGE(q,e,r, OUTER) {
 334     eback= edge_reverse(q,e);
 335     diff= edge_lengths[q][e] - edge_lengths[q][eback];
 336     cost += pow(diff,exponent_r);
 337   }
 338   return cost;
 339 }
 340
 341 /*---------- proportional edge length variation ----------*/
 342
 343   /*
 344    * Definition:
 345    *
 346    *    See the diagram above.
 347    *                                r
 348    *       cost    = ( |PQ| - |QR| )
 349    *           Q,e
 350    */
 351
 352 double prop_edge_length_variation_cost(const Vertices vertices, int section) {
 353   const double num_epsilon= 1e-6;
 354
 355   double cost= 0, exponent_r= 2;
 356   int q, e,r, eback;
 357
 358   FOR_EDGE(q,e,r, OUTER) {
 359     eback= edge_reverse(q,e);
 360     double le= edge_lengths[q][e];
 361     double leback= edge_lengths[q][eback];
 362     double diff= le - leback;
 363     double num= MIN(le, leback);
 364     cost += pow(diff / (num + num_epsilon), exponent_r);
 365   }
 366   return cost;
 367 }
 368
 369 /*---------- rim proximity cost ----------*/
 370
 371 static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) {
 372   /* By symmetry, nearest point on circle is the one with
 373    * the same angle subtended at the z axis. */
 374   oncircle[0]= p[0];
 375   oncircle[1]= p[1];
 376   oncircle[2]= 0;
 377   double mult= 1.0/ magnD(oncircle);
 378   oncircle[0] *= mult;
 379   oncircle[1] *= mult;
 380 }
 381
 382 double rim_proximity_cost(const Vertices vertices, int section) {
 383   double oncircle[3], cost=0;
 384   int v;
 385
 386   FOR_VERTEX(v, OUTER) {
 387     int y= v >> YSHIFT;
 388     int nominal_edge_distance= y <= Y/2 ? y : Y-1-y;
 389     if (nominal_edge_distance==0) continue;
 390
 391     find_nearest_oncircle(oncircle, vertices[v]);
 392
 393     cost +=
 394       vertex_mean_edge_lengths[v] *
 395       (nominal_edge_distance*nominal_edge_distance) /
 396       (hypotD2(vertices[v], oncircle) + 1e-6);
 397   }
 398   return cost;
 399 }
 400
 401 /*---------- noncircular rim cost ----------*/
 402
 403 double noncircular_rim_cost(const Vertices vertices, int section) {
 404   int vy,vx,v;
 405   double cost= 0.0;
 406   double oncircle[3];
 407
 408   FOR_RIM_VERTEX(vy,vx,v, OUTER) {
 409     find_nearest_oncircle(oncircle, vertices[v]);
 410
 411     double d2= hypotD2(vertices[v], oncircle);
 412     cost += d2*d2;
 413   }
 414   return cost;
 415 }
 416
 417 /*---------- overly sharp edge cost ----------*/
 418
 419   /*
 420    *
 421    *                Q `-_
 422    *              / |    `-_                           P'Q' ------ S'
 423    *             /  |       `-.                      _,' `.       .
 424    *            /   |           S                 _,'     :      .
 425    *           /    |      _,-'                _,'        :r    .r
 426    *          /     |  _,-'                R' '           `.   .
 427    *         /    , P '                        ` .   r     :  .
 428    *        /  ,-'                                  `  .   :
 429    *       /,-'                                          ` C'
 430    *      /'
 431    *     R
 432    *
 433    *
 434    *
 435    *  Let delta =  angle between two triangles' normals
 436    *
 437    *  Giving energy contribution:
 438    *
 439    *                                   2
 440    *    E             =  F   .  delta
 441    *     vd, edge PQ      vd
 442    */
 443
 444 double edge_angle_cost(const Vertices vertices, int section) {
 445   double pq1[D3], rp[D3], ps[D3], rp_2d[D3], ps_2d[D3], rs_2d[D3];
 446   double a,b,c,s,r;
 447   const double minradius_base= 0.2;
 448
 449   int pi,e,qi,ri,si, k;
 450 //  double our_epsilon=1e-6;
 451   double total_cost= 0;
 452
 453   FOR_EDGE(pi,e,qi, OUTER) {
 454 //    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
 455
 456     si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
 457     ri= EDGE_END2(pi,(e   +1)%V6);  if (ri<0) continue;
 458
 459     K {
 460       pq1[k]= -vertices[pi][k] + vertices[qi][k];
 461       rp[k]=  -vertices[ri][k] + vertices[pi][k];
 462       ps[k]=  -vertices[pi][k] + vertices[si][k];
 463     }
 464
 465     normalise(pq1,1,1e-6);
 466     xprod(rp_2d, rp,pq1); /* projects RP into plane normal to PQ */
 467     xprod(ps_2d, ps,pq1); /* likewise PS */
 468     K rs_2d[k]= rp_2d[k] + ps_2d[k];
 469     /* radius of circumcircle of R'P'S' from Wikipedia
 470      * `Circumscribed circle' */
 471     a= magnD(rp_2d);
 472     b= magnD(ps_2d);
 473     c= magnD(rs_2d);
 474     s= 0.5*(a+b+c);
 475     r= a*b*c / sqrt((a+b+c)*(a-b+c)*(b-c+a)*(c-a+b) + 1e-6);
 476
 477     double minradius= minradius_base + edge_angle_cost_circcircrat*(a+b);
 478     double deficit= minradius - r;
 479     if (deficit < 0) continue;
 480     double cost= deficit*deficit;
 481
 482     total_cost += cost;
 483   }
 484
 485   return total_cost;
 486 }
 487
 488 /*---------- small triangles cost ----------*/
 489
 490   /*
 491    *
 492    *                Q `-_
 493    *              / |    `-_
 494    *             /  |       `-.
 495    *            /   |           S
 496    *           /    |      _,-'
 497    *          /     |  _,-'
 498    *         /    , P '
 499    *        /  ,-'
 500    *       /,-'
 501    *      /'
 502    *     R
 503    *
 504    *  Let delta =  angle between two triangles' normals
 505    *
 506    *  Giving energy contribution:
 507    *
 508    *                                   2
 509    *    E             =  F   .  delta
 510    *     vd, edge PQ      vd
 511    */
 512
 513 double small_triangles_cost(const Vertices vertices, int section) {
 514   double pq[D3], ps[D3];
 515   double x[D3];
 516   int pi,e,qi,si, k;
 517 //  double our_epsilon=1e-6;
 518   double total_cost= 0;
 519
 520   FOR_EDGE(pi,e,qi, OUTER) {
 521 //    if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue;
 522
 523     si= EDGE_END2(pi,(e+V6-1)%V6);  if (si<0) continue;
 524
 525     K {
 526       pq[k]= vertices[qi][k] - vertices[pi][k];
 527       ps[k]= vertices[si][k] - vertices[pi][k];
 528     }
 529     xprod(x, pq,ps);
 530
 531     double cost= 1/(magnD2(x) + 0.01);
 532
 533 //double cost= pow(magnD(spqxpqr), 3);
 534 //assert(dot>=-1 && dot <=1);
 535 //double cost= 1-dot;
 536     total_cost += cost;
 537   }
 538
 539   return total_cost;
 540 }