chiark / gitweb /
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  */
46 #define L 10
47 #define W 10
48 #define N (L*W)
50 static double energy_function(const double vertices[N]) {
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    */
112 static gsl_multimin_fminimizer *minimiser;
116 {
117   gsl_multimin_fminimizer_alloc(