-static double edgewise_vertex_displacement_cost(const Vertices vertices) {
- static const l3_epsison= 1e-6;
-
- int pi,e,qi,ri,si, k;
- double m[D3], mprime[D3], b, d2, l, sigma_bd2_l3;
-
- 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= hypotD2(m, mprime);
- l= hypotD(vertices[ri][k] - vertices[si][k]);
- l3 = l*l*l + l3_epsilon;
-
- sigma_bd2_l3 += b * d2 / l3;
+double line_bending_adjcost(const Vertices vertices) {
+ static const double axb_epsilon= 1e-6;
+ static const double exponent_r= 3;
+
+ int pi,e,qi,ri, k;
+ double a[D3], b[D3], axb[D3];
+ double total_cost= 0;
+
+ FOR_EDGE(qi,e,ri) {
+ pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue;
+
+ K a[k]= -vertices[pi][k] + vertices[qi][k];
+ K b[k]= -vertices[qi][k] + vertices[ri][k];
+
+ xprod(axb,a,b);
+
+ double delta= atan2(magnD(axb) + axb_epsilon, dotprod(a,b));
+ double cost= pow(delta,exponent_r);
+
+ if (!e && !(qi & YMASK))
+ cost *= 10;
+
+ total_cost += cost;
+ }
+ return total_cost / (N / density);
+}
+
+/*---------- edge length variation ----------*/
+
+double edge_length_variation_cost(const Vertices vertices) {
+ double diff, cost= 0;
+ int v0, efwd,vfwd, eback;
+
+ FOR_EDGE(v0,efwd,vfwd) {
+ eback= edge_reverse(v0,efwd);
+ diff= edge_lengths[v0][efwd] - edge_lengths[v0][eback];
+ cost += diff*diff;
+ }
+ 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;
+
+ find_nearest_oncircle(oncircle, vertices[v]);
+
+ cost +=
+ vertex_mean_edge_lengths[v] *
+ (nominal_edge_distance*nominal_edge_distance) /
+ (hypotD2(vertices[v], oncircle) + 1e-6);
+ }
+ return cost;
+}
+
+/*---------- noncircular rim cost ----------*/
+
+double noncircular_rim_cost(const Vertices vertices) {
+ int vy,vx,v;
+ double cost= 0.0;
+ double oncircle[3];
+
+ FOR_RIM_VERTEX(vy,vx,v) {
+ find_nearest_oncircle(oncircle, vertices[v]);
+
+ double d2= hypotD2(vertices[v], oncircle);
+ cost += d2*d2;