+/*
+ * 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"
}
/*
- * edge descriptor f = 00 | e | y | x
- * 3 YBITS XBITS
+ * edge descriptor f = 0000 | e | y | x
+ * 3 YBITS XBITS
*
- * e is 0..5. The edge is edge e out of vertex (x,y).
+ * e is 0..6. The edge is edge e out of vertex (x,y), or if
+ * e==6 it's the `at end' value for the out edge iterator.
*
* BGL expects an undirected graph's edges to have two descriptors
- * each, one in each direction.
+ * each, one in each direction (otherwise e would be just 0..2).
*/
/*
*/
#define VMASK (YMASK|XMASK)
-#define ESHIFT (YBITS|XBITS)
+#define ESHIFT (YBITS+XBITS)
+
+using namespace boost;
+
+/*
+ * We iterate over edges in the following order:
+ *
+ * \#0 /1#
+ * \ /
+ * ___ 0 __
+ * #2 1 #3
+ * / \
+ * #4/ #5\ and finally #6 is V6
+ *
+ *
+ * This ordering permits the order-4 nodes at the strip's edge
+ * to have a contiguous edge iterator values. The iterator
+ * starts at #0 which is edge 2 (see mgraph.h), or #2 (edge 3).
+ */
+static const int oei_edge_delta[V6]=
+ /* 0 1 2 3 4 5 initial e
+ * #3 #1 #0 #2 #4 #5 initial ix
+ * #4 #2 #1 #3 #5 #6 next ix
+ * 4 3 1 0 5 V6 next e
+ */ {
+ 4<<ESHIFT, 2<<ESHIFT, -1<<ESHIFT,
+ -3<<ESHIFT, 1<<ESHIFT, (V6-5)<<ESHIFT
+};
+
+class OutEdgeIterator :
+ public iterator_facade<
+ OutEdgeIterator,
+ int const,
+ forward_traversal_tag
+> {
+ int f;
+ public:
+ void increment() {
+ //printf("incrementing f=%03x..",f);
+ f += oei_edge_delta[f>>ESHIFT];
+ //printf("%03x\n",f);
+ }
+ 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) {
+ //printf("constructed v=%02x e=%x f=%03x\n",v,e,f);
+ }
+
+ static int voe_min(int _v) { return (_v & YMASK) ? 2 : 3; }
+ static int voe_max(int _v) { return (_v & YMASK)==(Y-1) ? V6 : 4; }
+ static int voe_degree(int _v) { return RIM_VERTEX_P(_v) ? 4 : V6; }
+};
+
+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.
+ class Graph { }; // this is a dummy as our graph has no actual representation
+
+ // 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:
public virtual incidence_graph_tag,
public virtual vertex_list_graph_tag,
public virtual edge_list_graph_tag { };
-
- struct OutEdgeIncrable {
- int f;
- OutEdgeIncrable& operator++() { f += 1<<ESHIFT; return self; }
- OutEdgeIncrable(int v, int e) : f(v | (e << ESHIFT)) { }
- };
-
- struct graph_traits<Layout> {
+ template<>
+ 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&) { return f&VMASK; }
- inline int target(int f, const Layout&) { return EDGE_END2(f&VMASK, f>>ESHIFT); }
- inline std::pair<out_edge_iterator,out_edge_iterator>
- out_edges(int v, const Layout&) {
- 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&) { 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;
- inline std::pair<vertex_iterator,vertex_iterator>
- vertices(const Layout&) {
- return std::make_pair(vertex_iterator(0), vertex_iterator(N));
- }
- inline unsigned num_vertices(const Layout&) { return N; }
+ };
+
+ // Concept Graph:
+ inline int null_vertex() { return -1; }
-}
+ // Concept IncidenceGraph:
+ inline int source(int f, const Graph&) { return f&VMASK; }
+ inline int target(int f, const Graph&) {
+ int v2= EDGE_END2(f&VMASK, f>>ESHIFT);
+ //printf("traversed %03x..%02x\n",f,v2);
+ return v2;
+ }
+ inline std::pair<OutEdgeIterator,OutEdgeIterator>
+ out_edges(int v, const Graph&) {
+ return std::make_pair(OutEdgeIterator(v, OutEdgeIterator::voe_min(v)),
+ OutEdgeIterator(v, OutEdgeIterator::voe_max(v)));
+ }
+ inline unsigned out_degree(int v, const Graph&) {
+ return OutEdgeIterator::voe_degree(v);
+ }
-struct VertexIndexMap;
+ // 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; }
+};
-namespace boost {
- struct property_traits<VertexIndexMap> {
- // Concept Readable Property Map:
- typedef int value_type, reference, key_type;
- category
-class Boost
- }};
-
-void single_source_shortest_paths(int v1,
- const double edge_weights[/*f*/],
- ) {
- boost::dijkstra_shortest_paths
- (g, v1,
+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()).
-
-
-void all_pairs_shortest_paths(const Layout *g) {
+ distance_map(vertex_distances));
+}
+
+static int distances[N][N];
+
+void graph_layout_prepare() {
+ Graph g;
+ int v1, v2;
FOR_VERTEX(v1) {
-
- 0);
+ int *d= distances[v1];
+ FOR_VERTEX(v2) d[v2]= -1;
+ d[v1]= 0;
+ breadth_first_search
+ (g, v1,
+ vertex_index_map(identity_property_map()).
+ visitor(make_bfs_visitor(record_distances(d,on_tree_edge()))));
+ printf("%02x:",v1);
+ FOR_VERTEX(v2) printf(" %02x:%d",v2,d[v2]);
+ putchar('\n');
+ }
+ printf("---\n");
+ FOR_VERTEX(v1) {
+ int *d= distances[v1];
+ printf("%02x:",v1);
+ FOR_VERTEX(v2) printf(" %02x:%d",v2,d[v2]);
+ putchar('\n');
+ }
+}
+
+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[V6<<ESHIFT], vertex_distances[N],
+ double total_cost=0;
+ int v1,v2,e, nedges=0;
+ double totaledgelength=0, meanedgelength;
+
+ FOR_EDGE(v1,e,v2) {
+ totaledgelength += hypotD(v[v1], v[v2]);
+ nedges++;
+ }
+
+ meanedgelength= totaledgelength / nedges;
- /* weight_map(). ? */
- /*
+ FOR_VERTEX(v1) {
+ FOR_VERTEX(v2) {
+ if (v1 == v2) continue;
+
+ double d= hypotD(v[v1],v[v2]);
+
+ int dist= distances[v1][v2];
+ assert(dist>=0);
+
+ double s= dist * meanedgelength * 0.03;
-
- predecessor_map().
- distance_map()
+ double enoughdistance= d - s;
+ if (enoughdistance > 1e-6) continue;
+ /* energy = 1/2 stiffness deviation^2
+ * where stiffness = 1/d
+ */
+ double cost= pow(enoughdistance,4);
+
+ //double s2= s*s + d2_epsilon;
+ //double sd2= s2 / d2;
+ //double cost_contrib= a1*a2 * (sd2 - 1) / (d2*d2);
+ //double cost_contrib= sd2;
+ printf("layout %03x..%03x dist=%d mean=%g s=%g d=%g enough=%g"
+ " cost+=%g\n", v1,v2, dist, meanedgelength,
+ s,d, enoughdistance, cost);
+ total_cost += cost;
+ }
+ }
+ return total_cost;
+}
~ian / moebius2.git / blobdiff