*
* Vertices in strip are numbered as follows:
*
- * ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 ___ 3 ___ 4 __
- * W-1 W-1 0 0 0 0 0
+ * ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 ___ 3 ___ 4 __
+ * Y-1 Y-1 0 0 0 0 0
* / \ / \ / \ / \ / \ / \ / \
* / \ / \ / \ / \ / \ / \ / \
- * L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 ___ 3 ___ 4
- * W-2 W-2 W-2 1 1 1 1 1
+ * X-3 ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 ___ 3 ___ 4
+ * Y-2 Y-2 Y-2 1 1 1 1 1
* \ / \ / \ / \ / \ / \ / \ /
* \ / \ / \ / \ / \ / \ / \ /
- * ___ L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 ___ 3 __
- * W-3 W-3 W-3 2 2 2 2
- * / \ / \ / \ / \ / \ / \ / \
- * / \ / \ / \ / \ / \ / \ / \
- * . . . . . . . . . . . . . . . .
+ * ___ X-3 ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 ___ 3 __
+ * Y-3 Y-3 Y-3 2 2 2 2
+ *
+ * . . . . . . . . . . . . . . .
*
- * _ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 ___ 3
- * 1 1 1 1 W-2 W-2 W-2 W-2
- * / \ / \ / \ / \ / \ / \ / \ /
- * / \ / \ / \ / \ / \ / \ / \ /
- * L-5 ___ L-4 ___ L-3 ___ L-2 ___ L-1 ___ 0 ___ 1 ___ 2 __
- * 0 0 0 0 0 W-1 W-1 W-1
+ * X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 ___ 3
+ * 1 1 1 1 Y-2 Y-2 Y-2 Y-2
+ * \ / \ / \ / \ / \ / \ / \ /
+ * \ / \ / \ / \ / \ / \ / \ /
+ * ___ X-4 ___ X-3 ___ X-2 ___ X-1 ___ 0 ___ 1 ___ 2 __
+ * 0 0 0 0 Y-1 Y-1 Y-1
*
* Node x,y for
- * 0 <= x < L x = distance along L = length
- * 0 <= y < W y = distance across W = width
+ * 0 <= x < X x = distance along
+ * 0 <= y < Y y = distance across
*
* Vertices are in reading order from diagram above ie x varies fastest.
*
- * Where necessary, edges are counted anticlockwise from to-the-right:
+ * Y must be even. The actual location opposite (0,0) is (X-(Y-1)/2,0),
+ * and likewise opposite (0,Y-1) is ((Y-1)/2,0).
+ *
+ * To label edges, we counte anticlockwise[*] from to-the-right:
*
* \2 /1
* \ /
* / \
* 4/ 5\
*
+ * [*] That is, in the direction represented as anticlockwise for
+ * the vertices (0,*)..(4,*) in the diagram above; and of course
+ * that is clockwise for the vertices (X-5,*)..(X-1,*). The
+ * numbering has an actual discontinuity between (X-1,*) and (0,*).
+ *
* When we iterate over edges, we iterate first over vertices and then
* over edges 0 to 2, disregarding edges 3 to 5.
*/
-#define L 10
-#define W 10
-#define N (L*W)
+#define XBITS 4
+#define X (1<<XBITS)
+#define YBITS 4
+#define Y (1<<YBITS)
+
+/* vertex number: 0000 | y | x
+ * YBITS XBITS
+ */
+
+#define N (X*Y)
+#define XMASK (X-1)
+#define YSHIFT XBITS
+#define Y1 (1 << YSHIFT)
+#define YMASK (Y-1 << YSHIFT)
-static double energy_function(const double vertices[N][3]) {
/* Energy has the following parts right now:
- * - worst curvature radius
- * - total angle subtended
*/
/*
+ * Edgewise expected vertex displacement
*
- * Q
- * / | `-.
- * / | `-.
- * / P ------ S
- * /.-'
- * /'
- * R
- */
- /*
- * Curvature radius:
*
- * First flatten into plane _|_ PQ:
+ * Q `-_
+ * / | `-_
+ * R' - _ _ _/_ | `-.
+ * . / M - - - - - S
+ * . / | _,-'
+ * . / | _,-'
+ * . / , P '
+ * . / ,-'
+ * . /,-'
+ * . /'
+ * R
*
- * R' = (R-P) x (Q-P)/|Q-P|
- * S' = (S-P) x (Q-P)/|Q-P|
- * Q' = (Q-P) x (Q-P)/|Q-P| = 0 = P'
*
*
- * R'.
- * \.
- * \.
- * 0
- * |
- * |
- * |
- * S'
+ * Find R', the `expected' location of R, by
+ * reflecting S in M (the midpoint of QP).
*
- * Now radius formula (wikipedia `Circle'), P1=R', P2=0, P3=S':
+ * Let 2d = |RR'|
+ * b = |PQ|
+ * l = |RS|
*
- * |R'| * |S'| * | R' - S' |
- * r = -------------------------
- * 2 * | R' x S' |
+ * Giving energy contribution:
*
- * Energy
- * 2
- * E = C min r
- * r r j in edges j
+ * 2
+ * b d
+ * E = F . ----
+ * vd, edge PQ vd l
*
- */
- /*
- * Angle subtended:
+ * By symmetry, this calculation gives the same answer
+ * with R and S exchanged. Looking at the projection in
+ * the RMS plane:
*
- * -1 | R' . S' |
- * gamma = cos -----------
- * |R'| * |S'|
*
- * l = | P - Q |
+ * S'
+ * ,'
+ * ,'
+ * R' ,' 2d" = |SS'| = |RR'| = 2d
+ * `-._ ,'
+ * `-._ ,' By congruent triangles,
+ * ` M with M' = midpoint of RS,
+ * ,' `-._ |MM'| = |RR'|/2 = d
+ * ,' `-._
+ * ,' ` S So use
+ * ,' M' _ , - ' d = |MM'|
+ * ,' _ , - '
+ * R - '
*
- * Energy
+ * We choose this value for l (rather than |RM|+|MS|, say, or |RM|)
+ * because we want this symmetry and because we're happy to punish
+ * bending more than uneveness in the metric.
*
- * E = C sum gamma * l
- * gamma gamma j in edges j j
+ * In practice to avoid division by zero we'll add epsilon to l and
+ * the huge energy ought then to be sufficient for the model to
+ * avoid being close to R=S.
+ */
+
+#define V6 6
+
+static int edge_end2(unsigned v1, int e) {
+ /* The topology is equivalent to that of a square lattice with only
+ * half of the diagonals. Ie, the result of shearing the triangular
+ * lattice to make the lines of constant x vertical. This gives
+ * these six directions:
+ *
+ * 2 1
+ * | /
+ * |/
+ * 3--*--0
+ * /|
+ * / |
+ * 4 5
*
+ * This also handily makes vertical the numbering discontinuity,
+ * where the join happens.
*/
+ static const unsigned dx[V6]= { +1, +1, 0, -1, -1, 0 },
+ dy[V6]= { 0, +Y1, +Y1, 0, -Y1, -Y1 };
+ unsigned x, y;
+
+ y= (v1 & YMASK) + dy[e];
+ if (y & ~YMASK) return -1;
+
+ x= (v1 & XMASK) + dx[e];
+ if (x & ~XMASK) {
+ y= (Y-1)*Y1 - y;
+ x &= XMASK;;
+ }
+ return x | y;
+}
+
+#define D3 3
+
+#define FOR_VERTEX(v) \
+ for ((v)=0; (v)<N; (v)++)
+
+#define FOR_VPEDGE(v,e) \
+ for ((e)=0; (e)<V6; (e)++)
+
+#define EDGE_END2 edge_end2
+
+#define FOR_VEDGE(v1,e,v2) \
+ FOR_VPEDGE((v1),(e))
+ if (((v2)= EDGE_END2((v1),(e))) < 0) ; else
+
+#define FOR_EDGE(v1,e,v2) \
+ FOR_VERTEX((v1)) \
+ FOR_VEDGE((v1),(e),(v2))
+
+#define FOR_COORD(k) \
+ for ((k)=0; (k)<D3; (k)++)
+
+#define K FOR_COORD(k)
+
+static double hypotD(const double p[], const double q[]) {
+ int k;
+ double pq[D3];
+ gsl_vector v;
+
+ K pq[_k]= p[k] - q[k];
+ v.size= D3;
+ v.stride= 1;
+ v.vector.data= pq;
+ /* owner and block ought not to be used */
+
+ return gsl_blas_snrm2(&v);
+}
+
+#ifdef FP_FAST_FMA
+# define ffsqa(factor,term) fma((factor),(factor),(term))
+#else
+# define ffsqu(factor,term) ((factor)*(factor)+(term))
+#endif
+
+static double energy_function(const double vertices[N][D3]) {
+ int pi,e,qi,ri,si, k;
+ double m[D3], mprime[D3], b, d2, l, sigma_bd2_l;
+
+ FOR_EDGE(pi,e,qi) {
+ ri= EDGE_END2(pi,(e+1)%V6); if (r<0) continue;
+ si= EDGE_END2(pi,(e+5)%V6); if (s<0) continue;
+ assert(ri == EDGE_END2(qi,(e+2)%V6));
+ assert(si == EDGE_END2(qi,(e+4)%V6));
+
+ K m[k]= (vertices[pi][k] + vertices[qi][k]) * 0.5;
+ K mprime[k]= (vertices[ri][k] + vertices[si][k]) * 0.5;
+ b= hypotD(vertices[pi], vertices[qi]);
+ d2= 0; K d2= ffsqa(m[k] - mprime[k], d2);
+ l= hypotD(vertices[ri][k] - vertices[si][k]);
+
+ sigma_bd2_l += b * d2 / l;
+ }
+
+
static gsl_multimin_fminimizer *minimiser;
{
gsl_multimin_fminimizer_alloc(
+
~ian / moebius2.git / blobdiff