- double diff, cost= 0;
- int v0, efwd,vfwd, eback;
+ double diff, cost= 0, exponent_r= 2;
+ int q, e,r, eback;
+
+ FOR_EDGE(q,e,r) {
+ eback= edge_reverse(q,e);
+ diff= edge_lengths[q][e] - edge_lengths[q][eback];
+ cost += pow(diff,exponent_r);
+ }
+ return cost;
+}
+
+/*---------- rim proximity cost ----------*/
+
+static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) {
+ /* By symmetry, nearest point on circle is the one with
+ * the same angle subtended at the z axis. */
+ oncircle[0]= p[0];
+ oncircle[1]= p[1];
+ oncircle[2]= 0;
+ double mult= 1.0/ magnD(oncircle);
+ oncircle[0] *= mult;
+ oncircle[1] *= mult;
+}
+
+double rim_proximity_cost(const Vertices vertices) {
+ double oncircle[3], cost=0;
+ int v;
+
+ FOR_VERTEX(v) {
+ int y= v >> YSHIFT;
+ int nominal_edge_distance= y <= Y/2 ? y : Y-1-y;
+ if (nominal_edge_distance==0) continue;