- 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, OUTER) {
- 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);
+ 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;