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];
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;
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);
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 {
}
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;
+}