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:
*
- * 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'
+ * Q `-_
+ * / | `-_
+ * R' - _ _ _/_ | `-.
+ * . / M - - - - - S
+ * . / | _,-'
+ * . / | _,-'
+ * . / , P '
+ * . / ,-'
+ * . /,-'
+ * . /'
+ * R
*
*
- * R'.
- * \.
- * \.
- * 0
- * |
- * |
- * |
- * S'
*
- * Now radius formula (wikipedia `Circle'), P1=R', P2=0, P3=S':
+ * Find R', the `expected' location of R, by
+ * reflecting S in M (the midpoint of QP).
*
- * |R'| * |S'| * | R' - S' |
- * r = -------------------------
- * 2 * | R' x S' |
+ * Let d = |RR'|
+ * m = |PQ|
+ * l = |RS|
*
- * Energy
- * 2
- * E = C min r
- * r r j in edges j
+ * Giving energy contribution:
*
- */
- /*
- * Angle subtended:
+ * 2
+ * m d
+ * E = F . ----
+ * vd, edge PQ vd l
*
- * -1 | R' . S' |
- * gamma = cos -----------
- * |R'| * |S'|
+ * By symmetry, this calculation gives the same answer
+ * with R and S exchanged. Looking at the projection in
+ * the RMS plane:
*
- * l = | P - Q |
*
- * Energy
+ * S'
+ * ,'
+ * ,'
+ * R' ,' d" = |SS'| = |RR'| = d
+ * `-._ ,'
+ * `-._ ,'
+ * ` M
+ * ,' `-._
+ * ,' `-._
+ * ,' ` S
+ * ,'
+ * ,'
+ * R
*
- * E = C sum gamma * l
- * gamma gamma j in edges j j
+ * 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.
*
+ * 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.
*/
static gsl_multimin_fminimizer *minimiser;
{
gsl_multimin_fminimizer_alloc(
+