From 7bc1fc859f691dbd2194883ed3160411a082bbc8 Mon Sep 17 00:00:00 2001 From: Ian Jackson Date: Mon, 20 Oct 2008 00:38:16 +0100 Subject: [PATCH] wip before increase power of bendings --- energy.c | 78 ++++++++++++++++++++++++++---------------------------- mgraph.h | 7 +++-- minimise.h | 4 ++- 3 files changed, 46 insertions(+), 43 deletions(-) diff --git a/energy.c b/energy.c index b3af3fa..8bc125a 100644 --- a/energy.c +++ b/energy.c @@ -7,7 +7,11 @@ #include "mgraph.h" #include "parallel.h" -double vertex_areas[N], vertex_mean_edge_lengths[N], edge_lengths[N][V6]; +double vertex_mean_edge_lengths[N]; + +static double vertex_areas[N]; +static double edge_lengths[N][V6]; +static double rim_vertex_angles[N]; static double best_energy= DBL_MAX; @@ -29,7 +33,8 @@ typedef struct { static PreComputation *const precomps[]= { compute_edge_lengths, - compute_vertex_areas + compute_vertex_areas, + compute_rim_twist_angles }; static const CostContribution costs[]= { @@ -45,11 +50,11 @@ static const CostContribution costs[]= { #endif #if XBITS==4 -#define STOP_EPSILON 1.2e-3 +#define STOP_EPSILON 5e-3 COST( 3e3, vertex_displacement_cost) COST( 3e3, vertex_edgewise_displ_cost) COST( 0.2e3, rim_proximity_cost) - COST( 1e6, rim_twist_cost) + COST( 1e4, rim_twist_cost) COST( 1e12, noncircular_rim_cost) COST( 10e1, nonequilateral_triangles_cost) // COST( 1e1, small_triangles_cost) @@ -62,8 +67,9 @@ static const CostContribution costs[]= { COST( 3e3, vertex_displacement_cost) COST( 3e3, vertex_edgewise_displ_cost) COST( 2e-1, rim_proximity_cost) + COST( 3e3, rim_twist_cost) COST( 1e12, noncircular_rim_cost) - COST( 10e1, nonequilateral_triangles_cost) + COST( 3e2, nonequilateral_triangles_cost) // COST( 1e1, small_triangles_cost) // COST( 1e6, edge_angle_cost) #define EDGE_ANGLE_COST_CIRCCIRCRAT (0.5/1.7) @@ -74,6 +80,7 @@ static const CostContribution costs[]= { COST( 3e3, vertex_displacement_cost) COST( 3e3, vertex_edgewise_displ_cost) COST( 0.02e0, rim_proximity_cost) + COST( 1e4, rim_twist_cost) COST( 1e12, noncircular_rim_cost) COST( 10e1, nonequilateral_triangles_cost) // COST( 1e1, small_triangles_cost) @@ -183,7 +190,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 > %# e* |", tcost, tenergy); + if (pr) printf(/*" %# e >"*/ " %# e* |", /*tcost,*/ tenergy); *energy += tenergy; } @@ -464,7 +471,7 @@ static void find_nearest_oncircle(double oncircle[D3], const double p[D3]) { } double rim_proximity_cost(const Vertices vertices, int section) { - double oncircle[3], cost=0; + double oncircle[D3], cost=0; int v; FOR_VERTEX(v, OUTER) { @@ -500,42 +507,33 @@ double noncircular_rim_cost(const Vertices vertices, int section) { /*---------- rim contact angle rotation ----------*/ -/* - * - * vq ---- vs - * / / - * e0/ /e1 - * vp 0--- vr - * - */ +void compute_rim_twist_angles(const Vertices vertices, int section) { + double oncircle[D3], distance[D3]; + int vpy,vpx,v,k; + + FOR_NEAR_RIM_VERTEX(vpy,vpx,v, 1,OUTER) { + find_nearest_oncircle(oncircle, vertices[v]); + /* we are interested in the angle subtended at the rim, from the + * rim's point of view. */ + K distance[k]= vertices[v][k] - oncircle[k]; + double distance_positive_z= distance[3]; + double distance_radial_outwards= dotprod(distance, oncircle); + rim_vertex_angles[v]= atan2(distance_positive_z, distance_radial_outwards); + } +} double rim_twist_cost(const Vertices vertices, int section) { - const double our_epsilon=1e-6; - int vpy,vpx,vpi, vqi,vri,vsi, e0,e1, k; - double total_cost= 0.0; - double pq[D3], rs[D3], pq_x_rs[D3]; - - FOR_RIM_VERTEX(vpy,vpx,vpi, INNER) { - FOR_VPEDGE(e0) { - if (e0==0 || e0==3) continue; - vqi= EDGE_END2(vpi,e0); if (vqi<0) continue; - vri= EDGE_END2(vpi,0); assert(vri>=0); - e1= !vertices_span_join_p(vpi,vri) ? e0 : V6 - e0; - vsi= EDGE_END2(vri,e1); assert(vsi>=0); - - K pq[k]= -vertices[vpi][k] + vertices[vqi][k]; - K rs[k]= -vertices[vri][k] + vertices[vsi][k]; - - xprod(pq_x_rs, pq,rs); - double magndiv= edge_lengths[vpi][e0] * edge_lengths[vri][e1]; - - double cost= magnD2(pq_x_rs) / (magndiv*magndiv + our_epsilon); - -//fprintf(stderr,"rimtwist P=%03x Q=%03x R=%03x S=%03x e0=%d e1=%d %f\n", -// vpi,vqi,vri,vsi,e0,e1, cost); + double total_cost= 0; + int vpy,vpx,v0,v1; + + FOR_NEAR_RIM_VERTEX(vpy,vpx,v0, 1,OUTER) { + v1= EDGE_END2(v0,0); assert(v1!=0); + double delta= rim_vertex_angles[v0] - rim_vertex_angles[v1]; + if (delta < M_PI) delta += 2*M_PI; + if (delta > M_PI) delta -= 2*M_PI; - total_cost += cost; - } + double cost= pow(delta, 4); + total_cost += cost; } return total_cost; diff --git a/mgraph.h b/mgraph.h index 2267fdb..63fd5dc 100644 --- a/mgraph.h +++ b/mgraph.h @@ -129,10 +129,13 @@ int edge_reverse(int v1, int e); FOR_VERTEX((v1), loop) \ FOR_VEDGE((v1),(e),(v2)) -#define FOR_RIM_VERTEX(vy,vx,v, loop) \ - for ((vy)=0; (vy)