From: Ian Jackson Date: Thu, 21 Feb 2008 19:12:17 +0000 (+0000) Subject: wip new energy functions X-Git-Url: http://www.chiark.greenend.org.uk/ucgi/~ian/git?p=moebius2.git;a=commitdiff_plain;h=ca6581b23c843e3117d3450a4ecf13bd6b71b8a5 wip new energy functions --- diff --git a/common.c b/common.c index 62b981e..3270181 100644 --- a/common.c +++ b/common.c @@ -15,6 +15,14 @@ double magnD(const double pq[D3]) { return gsl_blas_dnrm2(&v); } +double magnD2(const double pq[D3]) { + int k; + double d2= 0; + + K d2= ffsqa(pq[k], d2); + return d2; +} + double hypotD(const double p[D3], const double q[D3]) { int k; double pq[D3]; @@ -43,6 +51,14 @@ void xprod(double r[D3], const double a[D3], const double b[D3]) { r[2]= a[0]*b[1] - a[1]*b[0]; } +void xprod_norm(double r[D3], const double a[D3], const double b[D3], + double absepsilon, double one) { + int k; + xprod(r,a,b); + double multby= one/(magnD(r) + absepsilon); + K r[k] *= multby; +} + double dotprod(const double a[D3], const double b[D3]) { int k; double result= 0; diff --git a/common.h b/common.h index 7e70a75..a11f566 100644 --- a/common.h +++ b/common.h @@ -36,8 +36,11 @@ double hypotD2(const double p[D3], const double q[D3]); double hypotD2plus(const double p[D3], const double q[D3], double add); double magnD(const double pq[D3]); +double magnD2(const double pq[D3]); void xprod(double r[D3], const double a[D3], const double b[D3]); double dotprod(const double a[D3], const double b[D3]); +void xprod_norm(double r[D3], const double a[D3], const double b[D3], + double absepsilon, double one); void flushoutput(void); void diee(const char *what); diff --git a/energy.c b/energy.c index 80ac7e8..cca777f 100644 --- a/energy.c +++ b/energy.c @@ -35,16 +35,19 @@ double compute_energy(const struct Vertices *vs) { if (XBITS==3) { COST( 3e2, line_bending_cost(vs->a)); COST( 1e3, edge_length_variation_cost(vs->a)); - COST( 0.2e3, rim_proximity_cost(vs->a)); - COST( 1e8, noncircular_rim_cost(vs->a)); + COST( 0.4e3, rim_proximity_cost(vs->a)); + COST( 1e6, edge_angle_cost(vs->a)); +// COST( 1e1, small_triangles_cost(vs->a)); + COST( 1e12, noncircular_rim_cost(vs->a)); stop_epsilon= 1e-6; } else if (XBITS==4) { COST( 3e2, line_bending_cost(vs->a)); COST( 3e3, edge_length_variation_cost(vs->a)); - COST( 3.8e1, rim_proximity_cost(vs->a)); // 5e1 is too much + COST( 9.0e1, rim_proximity_cost(vs->a)); // 5e1 is too much // 2.5e1 is too little // 0.2e1 grows compared to previous ? // 0.6e0 shrinks compared to previous ? + COST( 1e12, edge_angle_cost(vs->a)); COST( 1e12, noncircular_rim_cost(vs->a)); stop_epsilon= 1e-5; } else { @@ -259,3 +262,123 @@ double noncircular_rim_cost(const Vertices vertices) { } return cost; } + +/*---------- triangle bad normals cost ----------*/ + + /* + * + * Q `-_ + * / | `-_ + * / | `-. + * / | S + * / | _,-' + * / | _,-' + * / , P ' + * / ,-' + * /,-' + * /' + * R + * + * Let delta = angle between two triangles' normals + * + * Giving energy contribution: + * + * 2 + * E = F . delta + * vd, edge PQ vd + */ + +double edge_angle_cost(const Vertices vertices) { + double sq[D3], pq[D3], qr[D3], sqp[D3], pqr[D3], rs[D3]; + double x[D3], y[D3]; + int pi,e,qi,ri,si, k; +// double our_epsilon=1e-6; + double total_cost= 0; + + FOR_EDGE(pi,e,qi) { +// if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue; + + si= EDGE_END2(pi,(e+V6-1)%V6); if (si<0) continue; + ri= EDGE_END2(pi,(e +1)%V6); if (ri<0) continue; + + K { + sq[k]= vertices[si][k] - vertices[qi][k]; + pq[k]= vertices[pi][k] - vertices[qi][k]; + qr[k]= vertices[qi][k] - vertices[ri][k]; + } + xprod(sqp, pq,sq); + xprod(pqr, pq,qr); + + double dot= dotprod(sqp,pqr); + xprod(x,sqp,pqr); + + K rs[k]= -vertices[ri][k] + vertices[si][k]; + xprod(y, pq,rs); + + double delta= + pow(atan2(magnD(x), dot), 2.0) * pow(magnD2(pq), 2.0) / + (pow(magnD(y), 0.3) + 1e-6); + double cost= pow(delta, 2.0); + +//double cost= pow(magnD(spqxpqr), 3); +//assert(dot>=-1 && dot <=1); +//double cost= 1-dot; + total_cost += cost; + } + + return total_cost; +} + +/*---------- small triangles cost ----------*/ + + /* + * + * Q `-_ + * / | `-_ + * / | `-. + * / | S + * / | _,-' + * / | _,-' + * / , P ' + * / ,-' + * /,-' + * /' + * R + * + * Let delta = angle between two triangles' normals + * + * Giving energy contribution: + * + * 2 + * E = F . delta + * vd, edge PQ vd + */ + +double small_triangles_cost(const Vertices vertices) { + double pq[D3], ps[D3]; + double x[D3]; + int pi,e,qi,si, k; +// double our_epsilon=1e-6; + double total_cost= 0; + + FOR_EDGE(pi,e,qi) { +// if (!(RIM_VERTEX_P(pi) || RIM_VERTEX_P(qi))) continue; + + si= EDGE_END2(pi,(e+V6-1)%V6); if (si<0) continue; + + K { + pq[k]= vertices[qi][k] - vertices[pi][k]; + ps[k]= vertices[si][k] - vertices[pi][k]; + } + xprod(x, pq,ps); + + double cost= 1/(magnD2(x) + 0.01); + +//double cost= pow(magnD(spqxpqr), 3); +//assert(dot>=-1 && dot <=1); +//double cost= 1-dot; + total_cost += cost; + } + + return total_cost; +} diff --git a/minimise.h b/minimise.h index a076125..60efb3e 100644 --- a/minimise.h +++ b/minimise.h @@ -21,6 +21,8 @@ double line_bending_cost(const Vertices vertices); double noncircular_rim_cost(const Vertices vertices); double edge_length_variation_cost(const Vertices vertices); double rim_proximity_cost(const Vertices vertices); +double edge_angle_cost(const Vertices vertices); +double small_triangles_cost(const Vertices vertices); extern const char *input_file, *best_file; extern char *best_file_tmp;