#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;
static PreComputation *const precomps[]= {
compute_edge_lengths,
- compute_vertex_areas
+ compute_vertex_areas,
+ compute_rim_twist_angles
};
static const CostContribution costs[]= {
#endif
#if XBITS==4
-#define STOP_EPSILON 1.2e-3
- COST( 3e3, vertex_displacement_cost)
- COST( 3e3, vertex_edgewise_displ_cost)
- COST( 0.2e3, rim_proximity_cost)
+#define STOP_EPSILON 5e-3
+ COST( 3e4, vertex_displacement_cost) // NB this is probably wrong now
+ COST( 3e4, vertex_edgewise_displ_cost) // we have changed the power
+ COST( 2e2, rim_proximity_cost)
+ COST( 1e4, rim_twist_cost)
COST( 1e12, noncircular_rim_cost)
COST( 10e1, nonequilateral_triangles_cost)
// COST( 1e1, small_triangles_cost)
#endif
#if XBITS==5
-#define STOP_EPSILON 1.2e-4
- COST( 3e3, vertex_displacement_cost)
- COST( 3e3, vertex_edgewise_displ_cost)
+#define STOP_EPSILON 7e-4
+ COST( 3e4, vertex_displacement_cost)
+ COST( 3e4, 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)
#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( 3e4, vertex_displacement_cost)
+ COST( 3e4, vertex_edgewise_displ_cost)
+ COST( 2e-1, rim_proximity_cost)
+ COST( 1e3, rim_twist_cost)
COST( 1e12, noncircular_rim_cost)
COST( 10e1, nonequilateral_triangles_cost)
// COST( 1e1, small_triangles_cost)
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;
}
double vertex_displacement_cost(const Vertices vertices, int section) {
const double inv_lambda= 1.0/1; //2;
- const double delta= 4;
+ const double delta= 6;
int si,e,qi,ri;
double total_cost= 0;
double vertex_edgewise_displ_cost(const Vertices vertices, int section) {
const double inv_lambda= 1.0/1; //2;
- const double delta= 4;
+ const double delta= 6;
int pi,e,qi,ri,si, k;
double m[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) {
return cost;
}
+/*---------- rim contact angle rotation ----------*/
+
+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) {
+ 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;
+
+ double cost= pow(delta, 4);
+ total_cost += cost;
+ }
+
+ return total_cost;
+}
+
/*---------- overly sharp edge cost ----------*/
/*