From 18da2eccaac6795d7a241959a9f2ebaf35f597e9 Mon Sep 17 00:00:00 2001 From: Ian Jackson Date: Sun, 19 Oct 2008 22:37:09 +0100 Subject: [PATCH] wip energy functions --- energy.c | 219 ++++++++++++++++++++++++++++++++++++++++------------- minimise.h | 2 + 2 files changed, 167 insertions(+), 54 deletions(-) diff --git a/energy.c b/energy.c index 19f32e2..72c6fd3 100644 --- a/energy.c +++ b/energy.c @@ -45,30 +45,42 @@ static const CostContribution costs[]= { #endif #if XBITS==4 -#define STOP_EPSILON 1.2e-4 +#define STOP_EPSILON 1.2e-3 COST( 3e3, vertex_displacement_cost) + COST( 3e3, vertex_edgewise_displ_cost) COST( 0.2e3, rim_proximity_cost) + COST( 1e12, noncircular_rim_cost) + COST( 10e1, nonequilateral_triangles_cost) +// COST( 1e1, small_triangles_cost) // COST( 1e6, edge_angle_cost) #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7) -// COST( 1e1, small_triangles_cost) - COST( 1e12, noncircular_rim_cost) #endif #if XBITS==5 #define STOP_EPSILON 1.2e-4 - COST( 3e7, line_bending_cost) - COST( 10e2, prop_edge_length_variation_cost) - COST( 9.0e3, rim_proximity_cost) // 5e1 is too much - // 2.5e1 is too little - // 0.2e1 grows compared to previous ? - // 0.6e0 shrinks compared to previous ? + COST( 3e3, vertex_displacement_cost) + COST( 3e3, vertex_edgewise_displ_cost) + COST( 2e-1, rim_proximity_cost) + COST( 1e12, noncircular_rim_cost) + COST( 10e1, nonequilateral_triangles_cost) +// COST( 1e1, small_triangles_cost) +// COST( 1e6, edge_angle_cost) + #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7) +#endif -// COST( 1e12, edge_angle_cost) - #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.3) - COST( 1e18, noncircular_rim_cost) +#if XBITS==6 +#define STOP_EPSILON 1.2e-4 + COST( 3e3, vertex_displacement_cost) + COST( 3e3, vertex_edgewise_displ_cost) + COST( 0.02e0, rim_proximity_cost) + COST( 1e12, noncircular_rim_cost) + COST( 10e1, nonequilateral_triangles_cost) +// COST( 1e1, small_triangles_cost) +// COST( 1e6, edge_angle_cost) + #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7) #endif -#if XBITS>=6 /* nonsense follows but never mind */ +#if XBITS>=7 /* nonsense follows but never mind */ #define STOP_EPSILON 1e-6 COST( 3e5, line_bending_cost) COST( 10e2, edge_length_variation_cost) @@ -170,7 +182,7 @@ double compute_energy(const struct Vertices *vs) { static void addcost(double *energy, double tweight, double tcost, int pr) { double tenergy= tweight * tcost; - if (pr) printf(" %# e x %g > %# e* |", tcost, tweight, tenergy); + if (pr) printf(" %# e > %# e* |", tcost, tenergy); *energy += tenergy; } @@ -213,46 +225,117 @@ void compute_vertex_areas(const Vertices vertices, int section) { } } -/*---------- combined vertex displacement ----------*/ +/*---------- displacement of vertices across a midpoint ----------*/ /* - * At vertex Q considering edge direction e to R - * and corresponding opposite edge to P. + * Subroutine used where we have + * + * R - - - - - - - M . - - - - R' + * ` . + * ` . + * S * - * Let R' = Q + PQ + * and wish to say that the vector RM should be similar to MS + * or to put it another way S = M + RM + * + * NB this is not symmetric wrt R and S since it divides by + * |SM| but not |RM| so you probably want to call it twice. + * + * Details: + * + * Let R' = M + SM * D = R' - R - * delta - * [ -1 ] - * cost = [ lambda . ( D . PQ/|PQ| ) + | D x PQ/|PQ| | ] - * Q,e [ ------------------------------------------- ] - * [ |PQ| ] + * + * Then the (1/delta)th power of the cost is + * proportional to |D|, and + * inversely proportional to |SM| + * except that |D| is measured in a wierd way which counts + * distance in the same direction as SM 1/lambda times as much + * ie the equipotential surfaces are ellipsoids around R', + * lengthened by lambda in the direction of RM. + * + * So + * delta + * [ -1 ] + * cost = [ lambda . ( D . SM/|SM| ) + | D x SM/|SM| | ] + * R,S,M [ ------------------------------------------- ] + * [ |SM| ] + * + */ + +static double vertex_one_displ_cost(const double r[D3], const double s[D3], + const double midpoint[D3], + double delta, double inv_lambda) { + const double smlen2_epsilon= 1e-12; + double sm[D3], d[D3], ddot, dcross[D3]; + int k; + + K sm[k]= -s[k] + midpoint[k]; + K d[k]= midpoint[k] + sm[k] - r[k]; + ddot= dotprod(d,sm); + xprod(dcross, d,sm); + double smlen2= magnD2(sm); + double cost_basis= inv_lambda * ddot + magnD(dcross); + double cost= pow(cost_basis / (smlen2 + smlen2_epsilon), delta); + + return cost; +} + +/*---------- displacement of vertices opposite at a vertex ----------*/ + + /* + * At vertex Q considering edge direction e to R + * and corresponding opposite edge to S. + * + * This is vertex displacement as above with M=Q */ double vertex_displacement_cost(const Vertices vertices, int section) { const double inv_lambda= 1.0/1; //2; const double delta= 4; - const double pqlen2_epsilon= 1e-12; - int pi,e,qi,ri, k; - double pq[D3], d[D3], ddot, dcross[D3]; + int si,e,qi,ri; double total_cost= 0; FOR_EDGE(qi,e,ri, OUTER) { - pi= EDGE_END2(qi,(e+3)%V6); if (pi<0) continue; + si= EDGE_END2(qi,(e+3)%V6); if (si<0) continue; - K pq[k]= -vertices[pi][k] + vertices[qi][k]; - K d[k]= vertices[qi][k] + pq[k] - vertices[ri][k]; - ddot= dotprod(d,pq); - xprod(dcross, d,pq); - double pqlen2= magnD2(pq); - double cost_basis= inv_lambda * ddot + magnD(dcross); - double cost= pow(cost_basis / (pqlen2 + pqlen2_epsilon), delta); + total_cost += vertex_one_displ_cost(vertices[ri], vertices[si], vertices[qi], + delta, inv_lambda); + } + return total_cost; +} - total_cost += cost; +/*---------- displacement of vertices opposite at an edge ----------*/ + + /* + * At edge PQ considering vertices R and S (see diagram + * below for overly sharp edge cost). + * + * Let M = midpoint of PQ + */ + +double vertex_edgewise_displ_cost(const Vertices vertices, int section) { + const double inv_lambda= 1.0/1; //2; + const double delta= 4; + + int pi,e,qi,ri,si, k; + double m[D3]; + double total_cost= 0; + + FOR_EDGE(pi,e,qi, OUTER) { + si= EDGE_END2(pi,(e+V6-1)%V6); if (si<0) continue; + ri= EDGE_END2(pi,(e +1)%V6); if (ri<0) continue; + + K m[k]= 0.5 * (vertices[pi][k] + vertices[qi][k]); + + total_cost += vertex_one_displ_cost(vertices[ri], vertices[si], m, + delta, inv_lambda); } return total_cost; } + /*---------- at-vertex edge angles ----------*/ /* @@ -488,26 +571,9 @@ double edge_angle_cost(const Vertices vertices, int section) { /*---------- small triangles cost ----------*/ /* + * Consider a triangle PQS * - * Q `-_ - * / | `-_ - * / | `-. - * / | S - * / | _,-' - * / | _,-' - * / , P ' - * / ,-' - * /,-' - * /' - * R - * - * Let delta = angle between two triangles' normals - * - * Giving energy contribution: - * - * 2 - * E = F . delta - * vd, edge PQ vd + * Cost is 1/( area^2 ) */ double small_triangles_cost(const Vertices vertices, int section) { @@ -538,3 +604,48 @@ double small_triangles_cost(const Vertices vertices, int section) { return total_cost; } + +/*---------- nonequilateral triangles cost ----------*/ + + /* + * Consider a triangle PQR + * + * let edge lengths a=|PQ| b=|QR| c=|RP| + * + * predicted edge length p = 1/3 * (a+b+c) + * + * compute cost for each x in {a,b,c} + * + * + * cost = (x-p)^2 / p^2 + * PQR,x + */ + +double nonequilateral_triangles_cost(const Vertices vertices, int section) { + double pr[D3], abc[3]; + int pi,e0,e1,qi,ri, k,i; + double our_epsilon=1e-6; + double total_cost= 0; + + FOR_EDGE(pi,e0,qi, OUTER) { + e1= (e0+V6-1)%V6; + ri= EDGE_END2(pi,e1); if (ri<0) continue; + + K pr[k]= -vertices[pi][k] + vertices[ri][k]; + + abc[0]= edge_lengths[pi][e0]; /* PQ */ + abc[1]= edge_lengths[qi][e1]; /* QR */ + abc[2]= magnD(pr); + + double p= (1/3.0) * (abc[0]+abc[1]+abc[2]); + double p_inv2= 1/(p*p + our_epsilon); + + for (i=0; i<3; i++) { + double diff= (abc[i] - p); + double cost= diff*diff * p_inv2; + total_cost += cost; + } + } + + return total_cost; +} diff --git a/minimise.h b/minimise.h index f75bc57..f4eddfb 100644 --- a/minimise.h +++ b/minimise.h @@ -21,6 +21,7 @@ extern double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6]; extern const double edge_angle_cost_circcircrat; double vertex_displacement_cost(const Vertices vertices, int section); +double vertex_edgewise_displ_cost(const Vertices vertices, int section); double line_bending_cost(const Vertices vertices, int section); double noncircular_rim_cost(const Vertices vertices, int section); double edge_length_variation_cost(const Vertices vertices, int section); @@ -28,6 +29,7 @@ double prop_edge_length_variation_cost(const Vertices vertices, int section); double rim_proximity_cost(const Vertices vertices, int section); double edge_angle_cost(const Vertices vertices, int section); double small_triangles_cost(const Vertices vertices, int section); +double nonequilateral_triangles_cost(const Vertices vertices, int section); extern const char *input_file, *best_file; extern char *best_file_tmp; -- 2.30.2