return d2;
}
+double hypotD2plus(const double p[D3], const double q[D3], double d2) {
+ K d2= ffsqa(p[k] - q[k], d2);
+ return d2;
+}
+
#ifdef FP_FAST_FMA
-# define ffsqa(factor,term) fma((factor),(factor),(term))
+# define fma_fast fma
#else
-# define ffsqa(factor,term) ((factor)*(factor)+(term))
+# define fma_fast(f1,f2,t) ((f1)*(f2)+(t))
#endif
+#define ffsqa(factor,term) fma_fast((factor),(factor),(term))
static const l3_epsison= 1e-6;
distance_map(vertex_distances));
}
-double graph_layout_energy(const Layout *g) {
- /* For each (vi,vj) computes shortest path pij = vi..vj along edges,
- * and actual distance dij = vi-vj.
+double graph_layout_cost(const Layout *g, const double vertex_areas[N]) {
+ /* For each (vi,vj) computes shortest path s_ij = |vi..vj|
+ * along edges, and actual distance d_ij = |vi-vj|.
+ *
+ * We will also use the `vertex areas': for each vertex vi the
+ * vertex area a_vi is the mean area of the incident triangles.
+ * This is computed elsewhere.
*
* Energy contribution is proportional to
*
- * pij - dij
- * --------- . Sigma length
- * 2 e member (edges(vi) union edges(vj) e
- * dij
+ * -4 2
+ * a a . d . [ (s/d) - 1 ]
+ * vi vj
+ *
+ * (In practice we compute d^2+epsilon and use it for the
+ * divisions, to avoid division by zero.)
*/
- double edge_weights[N*V6], vertex_distances[N];
+ double edge_weights[N*V6], vertex_distances[N], total_cost;
int v1, e, f;
FOR_VEDGE_X(v1,e,
edge_weights[f]= hypotD(g.v[v1], g.v[v2]);
FOR_VERTEX(v1) {
+ double a1= vertex_areas[v1];
single_source_shortest_paths(v1, edge_weights, vertex_distances);
FOR_VERTEX(v2) {
-
-
-
-
- { ) {
-
-
- 0);
-
- /* weight_map(). ? */
- /*
-
-
- predecessor_map().
- distance_map()
-
-
-
+ double a2= vertex_areas[v2];
+ double d2= hypotD2plus(g->v[v1],g->v[v2], d2_epsilon);
+ double sd= vertex_distances[v2] / d2;
+ double sd2= sd*sd;
+ total_cost= fma_fast(a1*a2, (sd2 - 1.0)/(d2*d2), total_cost);
+ }
+ }
+ return total_cost;
+}