anneal.c

   1 /*
   2  * We try to find an optimal triangle grid
   3  *
   4  * Vertices in strip are numbered as follows:
   5  *
   6  *     ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   ___ 4   __
   7  *         W-1     W-1      0       0       0       0       0
   8  *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
   9  *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
  10  *     L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___  3  ___ 4
  11  *     W-2     W-2     W-2      1       1       1       1       1
  12  *       \    /  \    /  \    /  \    /  \    /  \    /  \    /
  13  *        \  /    \  /    \  /    \  /    \  /    \  /    \  /
  14  *     ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3   __
  15  *         W-3     W-3     W-3      2       2       2       2
  16  *        /  \    /  \    /  \    /  \    /  \    /  \    /  \
  17  *       /    \  /    \  /    \  /    \  /    \  /    \  /    \
  18  *    .   .   .   .   .   .   .   .   .   .   .   .   .   .   .   .
  19  *
  20  *      _ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   ___ 3
  21  *         1       1       1       1      W-2     W-2     W-2     W-2
  22  *       /  \    /  \    /  \    /  \    /  \    /  \    /  \    /
  23  *      /    \  /    \  /    \  /    \  /    \  /    \  /    \  /
  24  *    L-5 ___ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0   ___ 1   ___ 2   __
  25  *     0       0       0       0       0      W-1     W-1     W-1
  26  *
  27  * Node x,y for
  28  *   0 <= x < L     x = distance along    L = length
  29  *   0 <= y < W     y = distance across   W = width
  30  *
  31  * Vertices are in reading order from diagram above ie x varies fastest.
  32  *
  33  * Where necessary, edges are counted anticlockwise from to-the-right:
  34  *
  35  *                          \2   /1
  36  *                           \  /
  37  *                        ___ 0   __
  38  *                        3    1   0
  39  *                           /  \
  40  *                         4/   5\
  41  *
  42  * When we iterate over edges, we iterate first over vertices and then
  43  * over edges 0 to 2, disregarding edges 3 to 5.
  44  */
  45
  46 #define L 10
  47 #define W 10
  48 #define N (L*W)
  49
  50 static double energy_function(const double vertices[N][3]) {
  51   /* Energy has the following parts right now:
  52    *  - worst curvature radius
  53    *  - total angle subtended
  54    */
  55   /*
  56    *
  57    *                Q
  58    *              / | `-.
  59    *             /  |    `-.
  60    *            /   P ------ S
  61    *           /.-'
  62    *          /'
  63    *         R
  64    */
  65   /*
  66    * Curvature radius:
  67    *
  68    *  First flatten into plane _|_ PQ:
  69    *
  70    *   R' = (R-P) x (Q-P)/|Q-P|
  71    *   S' = (S-P) x (Q-P)/|Q-P|
  72    *   Q' = (Q-P) x (Q-P)/|Q-P| = 0 = P'
  73    *
  74    *
  75    *          R'.
  76    *             \.
  77    *               \.
  78    *                 0
  79    *                 |
  80    *                 |
  81    *                 |
  82    *                 S'
  83    *
  84    *  Now radius formula (wikipedia `Circle'), P1=R', P2=0, P3=S':
  85    *
  86    *          |R'| * |S'| * | R' - S' |
  87    *     r =  -------------------------
  88    *               2 * | R' x S' |
  89    *
  90    *  Energy
  91    *                            2
  92    *     E  = C  min           r
  93    *      r    r   j in edges   j
  94    *
  95    */
  96   /*
  97    * Angle subtended:
  98    *
  99    *                -1   | R' . S' |
 100    *    gamma =  cos     -----------
 101    *                     |R'| * |S'|
 102    *
 103    *    l = | P - Q |
 104    *
 105    *  Energy
 106    *
 107    *    E      = C      sum           gamma  * l
 108    *     gamma    gamma   j in edges       j    j
 109    *
 110    */
 111
 112 static gsl_multimin_fminimizer *minimiser;
 113
 114
 115
 116 {
 117   gsl_multimin_fminimizer_alloc(