Fix periodic boundary conditions in geometric_graph()

src/graph/generation/graph_geometric.hh

@@ -34,52 +34,50 @@ using namespace std;
 using namespace boost;
 
 
-void get_box(const vector<double>& p, double w, vector<int>& box)
+template <class Point, class Range>
+void get_box(const Point& p, double w, vector<int>& box,
+             const Range& ranges, bool periodic)
 {
     if (box.size() != p.size())
         box.resize(p.size());
     for (size_t i = 0; i < p.size(); ++i)
-        box[i] = int(p[i] / w);
+    {
+        box[i] = int(floor(p[i] / w));
+        if (periodic && p[i] == ranges[i].second)
+            box[i] -= 1;
+    }
 }
 
-template <class Point>
-double get_dist(const Point& p1, const Point& p2)
+template <class Point, class Range>
+double get_dist(const Point& p1, const Point& p2,
+                const Range& ranges, bool periodic)
 {
-    double r = 0;
+    double r = 0, diff, size;
     for (size_t i = 0; i < p1.size(); ++i)
-        r += pow(p1[i] - p2[i], 2);
+    {
+        diff = abs(p1[i] - p2[i]);
+        if (periodic)
+        {
+            size = abs(ranges[i].second - ranges[i].first);
+            diff = min(diff, abs(diff - size));
+        }
+        r += pow(diff, 2);
+    }
     return sqrt(r);
 }
 
-void periodic(int& x, int size)
-{
-    if (x >= int(size))
-        x -= size;
-    if (x < 0)
-        x += size;
-}
-
-
-bool is_boundary(vector<double>& point, double r,
-                 const vector<pair<double, double> > ranges)
+void periodic(vector<int>& box, const vector<pair<int,int> >& ranges)
 {
-    bool boundary = 0;
-    for (size_t j = 0; j < ranges.size(); ++j)
+    for (size_t i = 0; i < box.size(); ++i)
     {
-        if (point[j] - r < ranges[j].first)
-        {
-            point[j] += ranges[j].second - ranges[j].first;
-            boundary = true;
-        }
-        if (point[j] + r >= ranges[j].second)
-        {
-            point[j] -= ranges[j].second - ranges[j].first;
-            boundary = true;
-        }
+        if (box[i] >= ranges[i].second)
+            box[i] = ranges[i].first;
+        if (box[i] < ranges[i].first)
+            box[i] = ranges[i].second - 1;
     }
-    return boundary;
 }
 
+
 template <class Graph>
 bool is_adjacent(typename graph_traits<Graph>::vertex_descriptor v1,
                  typename graph_traits<Graph>::vertex_descriptor v2,
@@ -102,6 +100,18 @@ struct get_geometric
         int N = points.size();
         typename Pos::checked_t pos = upos.get_checked();
         vector<int> box;
+        vector<pair<int, int> > box_ranges;
+
+        double w = 2 * r;
+        if (periodic_boundary)
+        {
+            box_ranges.resize(ranges.size());
+            for (size_t i = 0; i < ranges.size(); ++i)
+            {
+                box_ranges[i].first = floor(ranges[i].first / w);
+                box_ranges[i].second = ceil(ranges[i].second / w);
+            }
+        }
 
         tr1::unordered_multimap<vector<int>,
                                 typename graph_traits<Graph>::vertex_descriptor,
@@ -112,20 +122,8 @@ struct get_geometric
             typename graph_traits<Graph>::vertex_descriptor v = add_vertex(g);
             pos[v].resize(points[i].size());
             copy(points[i].begin(), points[i].end(), pos[v].begin());
-
-            get_box(points[i], r, box);
+            get_box(points[i], w, box, ranges, periodic_boundary);
             boxes.insert(make_pair(box, v));
-
-            if (periodic_boundary)
-            {
-                //insert the same vertex at the opposite box
-                vector<double> point = points[i];
-                if (is_boundary(point, r, ranges))
-                {
-                    get_box(point, r, box);
-                    boxes.insert(make_pair(box, v));
-                }
-            }
         }
 
         #pragma omp parallel for default(shared) private(i, box) \
@@ -133,40 +131,33 @@ struct get_geometric
         for (int i = 0; i < N; ++i)
         {
             typename graph_traits<Graph>::vertex_descriptor v = vertex(i, g);
-            get_box(points[i], r, box);
 
+            get_box(points[i], w, box, ranges, periodic_boundary);
             for (size_t k = 0; k < pow(3, box.size()); ++k)
             {
                 for (size_t j = 0; j < box.size(); ++j)
                     box[j] += (int(k / pow(3, j)) % 3) - 1;
 
+                if (periodic_boundary)
+                    periodic(box, box_ranges);
+
                 typeof(boxes.begin()) iter, end;
                 for (tie(iter, end) = boxes.equal_range(box);
                      iter != end; ++iter)
                 {
                     typename graph_traits<Graph>::vertex_descriptor w =
                         iter->second;
-                    double d = get_dist(pos[v], pos[w]);
-
-                    if (periodic_boundary)
-                    {
-                        if (is_adjacent(v, w, g))
-                            continue;
-                        //get correct distance for boundary
-                        vector<double> point = pos[w];
-                        if (is_boundary(point, r, ranges))
-                            d = min(d, get_dist(pos[v], point));
-                    }
+                    double d = get_dist(pos[v], pos[w], ranges,
+                                        periodic_boundary);
 
-                    if (w > v && d <= r)
+                    if (w > v && d <= r &&
+                        (!periodic_boundary || !is_adjacent(v, w, g)))
                     {
                         #pragma omp critical
                         add_edge(v, w, g);
                     }
                 }
-
-                for (size_t j = 0; j < box.size(); ++j)
-                    box[j] -= (int(k / pow(3, j)) % 3) - 1;
+                get_box(points[i], w, box, ranges, periodic_boundary);
             }
         }
     }