+/*
+ * Everything that needs the Boost Graph Library and C++ templates etc.
+ * (and what a crazy set of stuff that all is)
+ */
+
+#include <math.h>
+
+#include <iterator>
+
+#include <boost/config.hpp>
+#include <boost/iterator/iterator_facade.hpp>
+#include <boost/graph/graph_traits.hpp>
+#include <boost/graph/graph_concepts.hpp>
+#include <boost/graph/dijkstra_shortest_paths.hpp>
+#include <boost/graph/properties.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+#include <boost/iterator/iterator_categories.hpp>
+
extern "C" {
+#include "bgl.h"
#include "mgraph.h"
}
*/
/*
- * We use BGL's implementation of Johnson All Pairs Shortest Paths
+ * We use BGL's implementation of Dijkstra's single source shortest
+ * paths. We really want all pairs shortest paths, so Johnson All
+ * Pairs Shortest Paths would seem sensible. But actually Johnson's
+ * algorithm is just a wrapper around Dijkstra's; the extra
+ * functionality is just to deal with -ve edge weights, which we don't
+ * have. So we can use Dijkstra directly and save some cpu (and some
+ * code: we don't have to supply all of the machinery needed for
+ * Johnson's invocation of Bellman-Ford). The overall time cost is
+ * O(VE log V); I think the space used is O(E).
*/
#define VMASK (YMASK|XMASK)
#define ESHIFT (YBITS|XBITS)
-#define FMAX ((5 << ESHIFT) | VMASK)
+
+class Graph { }; // this is a dummy as our graph has no actual representation
+
+using namespace boost;
+
+struct OutEdgeIterator :
+ public iterator_facade<
+ OutEdgeIterator,
+ int const,
+ forward_traversal_tag
+> {
+ int f;
+ void increment() { f += 1<<ESHIFT; }
+ bool equal(OutEdgeIterator const& other) const { return f == other.f; }
+ int const& dereference() const { return f; }
+ OutEdgeIterator() { }
+ OutEdgeIterator(int _f) : f(_f) { }
+ OutEdgeIterator(int v, int e) : f(e << ESHIFT | v) { }
+};
+
+typedef counting_iterator<int> VertexIterator;
namespace boost {
- // We make Layout a model of various BGL Graph concepts.
- // This mainly means that graph_traits<Layout> has lots of stuff.
+ // We make Graph a model of various BGL Graph concepts.
+ // This mainly means that graph_traits<Graph> has lots of stuff.
// First, some definitions used later:
struct layout_graph_traversal_category :
public virtual incidence_graph_tag,
- public virtual edge_list_graph_tag
- public virtual vertex_list_graph_tag { };
-
- struct OutEdgeIncrable {
- int f;
- OutEdgeIncrable operator++(OutEdgeIncrable f) { return f + 1<<ESHIFT; }
- OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
- }
-
- struct graph_traits<Layout> {
+ public virtual vertex_list_graph_tag,
+ public virtual edge_list_graph_tag { };
+ struct graph_traits<Graph> {
// Concept Graph:
-
typedef int vertex_descriptor; /* vertex number, -1 => none */
typedef int edge_descriptor; /* see above */
typedef undirected_tag directed_category;
- typedef disallow_parallel_ege edge_parallel_category;
+ typedef disallow_parallel_edge_tag edge_parallel_category;
typedef layout_graph_traversal_category traversal_category;
- inline int null_vertex() { return -1; }
// Concept IncidenceGraph:
-
- typedef counting_iterator<OutEdgeIncrable,
- forward_iterator_tag> out_edge_iterator;
- typedef int degree_size_type;
-
- inline int source(int f, const Layout& g) { return f&VMASK; }
- inline int target(int f, const Layout& g) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
- inline std::pair<typename graph_traits<Layout>::out_edge_iterator,
- typename graph_traits<Layout>::out_edge_iterator>
- out_edges(int v, const Layout& g) {
- return std::make_pair(out_edge_iterator(OutEdgeIncrable(v, VE_MIN(v))),
- out_edge_iterator(OutEdgeIncrable(v, VE_MAX(v))));
- }
- out_degree(int v, const Layout& g) { return VE_MAX(v) - VE_MIN(v); }
+ typedef OutEdgeIterator out_edge_iterator;
+ typedef unsigned degree_size_type;
// Concept VertexListGraph:
- typedef counting_iterator<int> vertex_iterator;
+ typedef VertexIterator vertex_iterator;
+ typedef unsigned vertices_size_type;
+ };
-}
+ // Concept Graph:
+ inline int null_vertex() { return -1; }
-void calculate_layout_energy(const Layout*) {
-
- FOR_VERTEX(v1) {
- boost::dijkstra_shortest_paths(g, v1, 0);
-
- /* weight_map(). ? */
- /* vertex_index_map(vimap). */
+ // Concept IncidenceGraph:
+ inline int source(int f, const Graph&) { return f&VMASK; }
+ inline int target(int f, const Graph&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
+ inline std::pair<OutEdgeIterator,OutEdgeIterator>
+ out_edges(int v, const Graph&) {
+ return std::make_pair(OutEdgeIterator(v, VE_MIN(v)),
+ OutEdgeIterator(v, VE_MAX(v)));
+ }
+ inline unsigned out_degree(int v, const Graph&) {
+ return VE_MAX(v) - VE_MIN(v);
+ }
-
- predecessor_map().
- distance_map()
+ // Concept VertexListGraph:
+ inline std::pair<VertexIterator,VertexIterator> vertices(const Graph&) {
+ return std::make_pair(VertexIterator(0), VertexIterator(N));
+ }
+ inline unsigned num_vertices(const Graph&) { return N; }
+};
+static void single_source_shortest_paths(int v1,
+ const double edge_weights[/*f*/],
+ double vertex_distances[/*v*/]) {
+ Graph g;
+ dijkstra_shortest_paths(g, v1,
+ weight_map(edge_weights).
+ vertex_index_map(identity_property_map()).
+ distance_map(vertex_distances));
+}
+
+double graph_layout_cost(const Vertices v, 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
+ *
+ * -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.)
+ */
+ static const double d2_epsilon= 1e-6;
+
+ double edge_weights[N*V6], vertex_distances[N], total_cost=0;
+ int v1,v2,e,f;
+ FOR_VERTEX(v1)
+ FOR_VEDGE_X(v1,e,v2,
+ f= v1 | e << ESHIFT,
+ edge_weights[f]= NAN)
+ edge_weights[f]= hypotD(v[v1], v[v2]);
+
+ FOR_VERTEX(v1) {
+ double a1= vertex_areas[v1];
+ single_source_shortest_paths(v1, edge_weights, vertex_distances);
+ FOR_VERTEX(v2) {
+ double a2= vertex_areas[v2];
+ double d2= hypotD2plus(v[v1],v[v2], d2_epsilon);
+ double sd= vertex_distances[v2] / d2;
+ double sd2= sd*sd;
+ total_cost += a1*a2 * (sd2 - 1) / (d2*d2);
+ }
+ }
+ return total_cost;
+}