Implement periodic delaunay triangulation

src/graph/generation/graph_generation.cc

@@ -105,7 +105,7 @@ void edge_property_union(GraphInterface& ugi, GraphInterface& gi,
                          boost::any p_vprop, boost::any p_eprop,
                          boost::any uprop, boost::any prop);
 void triangulation(GraphInterface& gi, python::object points, boost::any pos,
-                   string type);
+                   string type, bool periodic);
 
 using namespace boost::python;
 

src/graph/generation/graph_triangulation.cc

@@ -20,6 +20,7 @@
 //  follow the requirements of the GNU GPL in regard to all of the
 //  software in the executable aside from CGAL.
 
+#include <CGAL/version.h>
 #include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
 #include <CGAL/Triangulation_3.h>
 #include <CGAL/Delaunay_triangulation_3.h>
@@ -27,6 +28,7 @@ typedef CGAL::Exact_predicates_inexact_constructions_kernel Kernel;
 typedef CGAL::Triangulation_3<Kernel> SimpleTriangulation;
 typedef CGAL::Delaunay_triangulation_3<Kernel> DelaunayTriangulation;
 
+
 namespace std
 {
 bool operator==(const SimpleTriangulation::Vertex& a,
@@ -36,19 +38,35 @@ bool operator==(const SimpleTriangulation::Vertex& a,
 }
 }
 
+// periodic triangulation is only available in more recent versions of CGAL
+#if (CGAL_VERSION_NR >= 1030500000)
+#include <CGAL/Periodic_3_triangulation_traits_3.h>
+#include <CGAL/Periodic_3_Delaunay_triangulation_3.h>
+typedef CGAL::Periodic_3_triangulation_traits_3<Kernel> GT;
+typedef CGAL::Periodic_3_Delaunay_triangulation_3<GT>
+    PeriodicDelaunayTriangulation;
+namespace std
+{
+bool operator==(const PeriodicDelaunayTriangulation::Vertex& a,
+                const PeriodicDelaunayTriangulation::Vertex& b)
+{
+    return a.point() == b.point();
+}
+}
+#endif
+
 #include "graph.hh"
 #include "graph_util.hh"
 #include "graph_filtering.hh"
 #include "graph_triangulation.hh"
 #include "numpy_bind.hh"
 
-
 using namespace std;
 using namespace boost;
 using namespace graph_tool;
 
 void triangulation(GraphInterface& gi, python::object points, boost::any pos,
-                   string type)
+                   string type, bool periodic)
 {
     UndirectedAdaptor<GraphInterface::multigraph_t> g(gi.GetGraph());
     multi_array_ref<double,2> points_array = get_array<double,2>(points);
@@ -57,9 +75,27 @@ void triangulation(GraphInterface& gi, python::object points, boost::any pos,
     pos_type_t pos_map = any_cast<pos_type_t>(pos);
 
     if (type == "simple")
+    {
         get_triangulation<SimpleTriangulation>()(g, points_array, pos_map);
+    }
     else if (type == "delaunay")
-        get_triangulation<DelaunayTriangulation>()(g, points_array,
-                                                   pos_map);
+    {
+        if (!periodic)
+        {
+            get_triangulation<DelaunayTriangulation>()(g, points_array,
+                                                       pos_map);
+        }
+        else
+        {
+#if (CGAL_VERSION_NR >= 1030500000)
+            get_triangulation<PeriodicDelaunayTriangulation>()(g, points_array,
+                                                               pos_map);
+#else
+            throw ValueException("Periodic Delaunay triangulation is only "
+                                 "available with versions of CGAL newer than "
+                                 "3.5.0.");
+#endif
+        }
+    }
     gi.ReIndexEdges();
 }

src/graph/generation/graph_triangulation.hh

@@ -89,7 +89,7 @@ struct get_triangulation
             {
                 typename graph_traits<Graph>::vertex_descriptor target
                     = iter->second;
-                if (!is_adjacent(_source, target, _g))
+                if (!is_adjacent(_source, target, _g) && _source != target)
                     add_edge(_source, target, _g);
             }
             return *this;

src/graph_tool/generation/__init__.py

@@ -527,7 +527,7 @@ def graph_union(g1, g2, props=[], include=False):
         return g1
 
 @_limit_args({"type":["simple", "delaunay"]})
-def triangulation(points, type="simple"):
+def triangulation(points, type="simple", periodic=False):
     r"""
     Generate a 2D or 3D triangulation graph from a given point set.
 
@@ -538,6 +538,9 @@ def triangulation(points, type="simple"):
         is the number of points, and d is the space dimension (either 2 or 3).
     type : string (optional, default: 'simple')
         Type of triangulation. May be either 'simple' or 'delaunay'.
+    periodic : bool (optional, default: False)
+        If True, periodic boundary conditions will be used. This is parameter is
+        valid only for type="delaunay", and is otherwise ignored.
 
     Returns
     -------
@@ -623,6 +626,6 @@ def triangulation(points, type="simple"):
     g = Graph(directed=False)
     pos = g.new_vertex_property("vector<double>")
     libgraph_tool_generation.triangulation(g._Graph__graph, points,
-                                           _prop("v", g, pos), type)
+                                           _prop("v", g, pos), type, periodic)
     return g, pos