Add support for triangulation() in generation module

This also includes the library CGAL as a dependency.

configure.ac

@@ -175,6 +175,11 @@ if test "$BOOST_GRAPH_LIB" = ""; then
 fi
 [CPPFLAGS="${CPPFLAGS} ${BOOST_CPPFLAGS}"]
 
+dnl CGAL
+AC_CHECK_LIB(CGAL,main)
+[CGAL_LIBADD="-lCGAL"]
+AC_SUBST(CGAL_LIBADD)
+
 dnl Checks for header files.
 
 dnl numpy
@@ -207,6 +212,11 @@ AC_CHECK_HEADER([expat.h],
                 [],
                 [AC_MSG_ERROR([Expat header not found])])
 
+dnl cgal
+AC_CHECK_HEADER([CGAL/version.h],
+                [],
+                [AC_MSG_ERROR([CGAL headers not found])])
+
 dnl Checks for typedefs, structures, and compiler characteristics.
 dnl Checks for library functions.
 

src/graph/generation/Makefile.am

@@ -10,7 +10,7 @@ libgraph_tool_generation_LTLIBRARIES = libgraph_tool_generation.la
 
 libgraph_tool_generation_la_includedir = $(pythondir)/graph_tool/include
 
-libgraph_tool_generation_la_LIBADD = $(MOD_LIBADD)
+libgraph_tool_generation_la_LIBADD = $(MOD_LIBADD) $(CGAL_LIBADD)
 
 libgraph_tool_generation_la_LDFLAGS = $(MOD_LDFLAGS)
 
@@ -21,14 +21,12 @@ libgraph_tool_generation_la_SOURCES = \
     graph_line_graph.cc \
     graph_union.cc \
     graph_union_vprop.cc \
-    graph_union_eprop.cc
-
-
-
-
+    graph_union_eprop.cc \
+    graph_triangulation.cc
 
 libgraph_tool_generation_la_include_HEADERS = \
     graph_generation.hh \
     graph_rewiring.hh \
     graph_predecessor.hh \
-    graph_union.hh
\ No newline at end of file
+    graph_union.hh \
+    graph_triangulation.hh

src/graph/generation/graph_generation.cc

@@ -104,6 +104,8 @@ void vertex_property_union(GraphInterface& ugi, GraphInterface& gi,
 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);
 
 using namespace boost::python;
 
@@ -116,4 +118,5 @@ BOOST_PYTHON_MODULE(libgraph_tool_generation)
     def("graph_union", &graph_union);
     def("vertex_property_union", &vertex_property_union);
     def("edge_property_union", &edge_property_union);
+    def("triangulation", &triangulation);
 }

src/graph/generation/graph_triangulation.cc

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007  Tiago de Paula Peixoto <tiago@forked.de>
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 3
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+//  As a special exception, you have permission to link this program
+//  with the CGAL library and distribute executables, as long as you
+//  follow the requirements of the GNU GPL in regard to all of the
+//  software in the executable aside from CGAL.
+
+#include <CGAL/Exact_predicates_inexact_constructions_kernel.h>
+#include <CGAL/Triangulation_3.h>
+#include <CGAL/Delaunay_triangulation_3.h>
+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,
+                const SimpleTriangulation::Vertex& b)
+{
+    return a.point() == b.point();
+}
+}
+
+#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)
+{
+    UndirectedAdaptor<GraphInterface::multigraph_t> g(gi.GetGraph());
+    multi_array_ref<double,2> points_array = get_array<double,2>(points);
+    typedef property_map_type::apply
+        <vector<double>, GraphInterface::vertex_index_map_t>::type pos_type_t;
+    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);
+    gi.ReIndexEdges();
+}

src/graph/generation/graph_triangulation.hh

+// graph-tool -- a general graph modification and manipulation thingy
+//
+// Copyright (C) 2007  Tiago de Paula Peixoto <tiago@forked.de>
+//
+// This program is free software; you can redistribute it and/or
+// modify it under the terms of the GNU General Public License
+// as published by the Free Software Foundation; either version 3
+// of the License, or (at your option) any later version.
+//
+// This program is distributed in the hope that it will be useful,
+// but WITHOUT ANY WARRANTY; without even the implied warranty of
+// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+// GNU General Public License for more details.
+//
+// You should have received a copy of the GNU General Public License
+// along with this program. If not, see <http://www.gnu.org/licenses/>.
+
+//  As a special exception, you have permission to link this program
+//  with the CGAL library and distribute executables, as long as you
+//  follow the requirements of the GNU GPL in regard to all of the
+//  software in the executable aside from CGAL.
+
+#ifndef GRAPH_TRIANGULATION_HH
+#define GRAPH_TRIANGULATION_HH
+
+#include <tr1/unordered_map>
+#include <boost/functional/hash.hpp>
+
+#include "graph_util.hh"
+
+namespace graph_tool
+{
+using namespace std;
+using namespace boost;
+
+
+template <class Graph>
+bool is_adjacent(typename graph_traits<Graph>::vertex_descriptor v1,
+                 typename graph_traits<Graph>::vertex_descriptor v2,
+                 Graph& g)
+{
+    typename graph_traits<Graph>::out_edge_iterator e, e_end;
+    for (tie(e, e_end) = out_edges(v1, g); e != e_end; ++e)
+        if (target(*e, g) == v2)
+            return true;
+    return false;
+}
+
+struct hash_point
+{
+    template <class Vertex>
+    std::size_t operator()(const Vertex& v) const
+    {
+        size_t seed = 42;
+        hash_combine(seed, v.point().x());
+        hash_combine(seed, v.point().y());
+        hash_combine(seed, v.point().z());
+        return seed;
+    }
+};
+
+template <class Triang>
+struct get_triangulation
+{
+    // this will insert edges in the graph
+    template <class Graph, class VertexMap>
+    class edge_inserter
+    {
+    public:
+        typedef output_iterator_tag iterator_category;
+        typedef typename graph_traits<Graph>::vertex_descriptor value_type;
+        typedef size_t difference_type;
+        typedef typename graph_traits<Graph>::vertex_descriptor* pointer;
+        typedef typename graph_traits<Graph>::vertex_descriptor& reference;
+
+        edge_inserter(Graph& g, const typename VertexMap::key_type& v,
+                      VertexMap& vertex_map): _g(g), _vertex_map(vertex_map),
+                                              _source(vertex_map[v]) {}
+
+        edge_inserter& operator++() { return *this; }
+        edge_inserter& operator++(int) { return *this; }
+        edge_inserter& operator*() { return *this; }
+
+        template <class Vertex>
+        edge_inserter& operator=(const Vertex& v)
+        {
+            typename VertexMap::iterator iter = _vertex_map.find(*v);
+            if (iter != _vertex_map.end())
+            {
+                typename graph_traits<Graph>::vertex_descriptor target
+                    = iter->second;
+                if (!is_adjacent(_source, target, _g))
+                    add_edge(_source, target, _g);
+            }
+            return *this;
+        }
+
+    private:
+        Graph& _g;
+        VertexMap _vertex_map;
+        typename graph_traits<Graph>::vertex_descriptor _source;
+    };
+
+    template <class Graph, class Points, class PosMap>
+    void operator()(Graph& g, Points& points, PosMap pos) const
+    {
+        typedef tr1::unordered_map
+            <typename Triang::Vertex,
+             typename graph_traits<Graph>::vertex_descriptor,
+             hash_point> vertex_map_t;
+        vertex_map_t vertex_map;
+
+        Triang T;
+        for (size_t i = 0; i < points.shape()[0]; ++i)
+        {
+            typename Triang::Point p(points[i][0], points[i][1], points[i][2]);
+            typename graph_traits<Graph>::vertex_descriptor v = add_vertex(g);
+            vertex_map[*T.insert(p)] = v;
+            pos[v].resize(3);
+            for (size_t j = 0; j < 3; ++j)
+                pos[v][j] = points[i][j];
+        }
+
+        typename Triang::Vertex_iterator v_begin, v_end;
+        v_begin = T.vertices_begin();
+        v_end = T.vertices_end();
+        while (v_begin != v_end)
+        {
+            typename Triang::Vertex_handle v = v_begin++;
+            if (vertex_map.find(*v) == vertex_map.end())
+                continue;
+            edge_inserter<Graph, vertex_map_t> insert(g, *v, vertex_map);
+            T.incident_vertices(v, insert);
+        }
+    }
+
+};
+
+} // namespace graph_tool
+
+#endif // GRAPH_TRIANGULATION_HH

src/graph_tool/generation/__init__.py

@@ -31,6 +31,7 @@ Summary
    predecessor_tree
    line_graph
    graph_union
+   triangulation
 
 Contents
 ++++++++
@@ -39,11 +40,11 @@ Contents
 from .. dl_import import dl_import
 dl_import("import libgraph_tool_generation")
 
-from .. core import Graph, _check_prop_scalar, _prop
+from .. core import Graph, _check_prop_scalar, _prop, _limit_args
 import sys, numpy, numpy.random
 
 __all__ = ["random_graph", "random_rewire", "predecessor_tree", "line_graph",
-           "graph_union"]
+           "graph_union", "triangulation"]
 
 def _corr_wrap(i, j, corr):
     return corr(i[1], j[1])
@@ -485,3 +486,97 @@ def graph_union(g1, g2, props=[], include=False):
         return g1, n_props
     else:
         return g1
+
+@_limit_args({"type":["simple", "delaunay"]})
+def triangulation(points, type="simple"):
+    r"""
+    Generate a 2D or 3D triangulation graph from a given point set.
+
+    Parameters
+    ----------
+    points : :class:`~numpy.ndarray`
+        Point set for the triangulation. It may be either a N x d array, where N
+        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'.
+
+    Returns
+    -------
+    triangulation_graph : :class:`~graph_tool.Graph`
+        The generated graph.
+    pos : :class:`~graph_tool.PropertyMap`
+        Vertex property map with the Cartesian coordinates.
+
+    See Also
+    --------
+    random_graph: random graph generation
+
+    Notes
+    -----
+
+    A triangulation [cgal_triang]_ is division of the convex hull of a point set
+    into triangles, using only that set as triangle vertices.
+
+    In simple triangulations (`type="simple"`), the insertion of a point is done
+    by locating a face that contains the point, and splitting this face into
+    three new faces (the order of insertion is therefore important). If the
+    point falls outside the convex hull, the triangulation is restored by
+    flips. Apart from the location, insertion takes a time O(1). This bound is
+    only an amortized bound for points located outside the convex hull.
+
+    Delaunay triangulations (`type="delaunay"`) have the specific empty sphere
+    property, that is, the circumscribing sphere of each cell of such a
+    triangulation does not contain any other vertex of the triangulation in its
+    interior. These triangulations are uniquely defined except in degenerate
+    cases where five points are co-spherical. Note however that the CGAL
+    implementation computes a unique triangulation even in these cases.
+
+    Examples
+    --------
+    >>> from numpy.random import seed, random
+    >>> seed(42)
+    >>> points = random((500,2))
+    >>> for i in xrange(points.shape[0]):
+    ...     x,y = 1,1
+    ...     while sqrt(x**2 + y**2) > 0.5:
+    ...          x, y = random()-0.5, random()-0.5
+    ...     points[i,:] = [x,y]
+    >>> g, pos = gt.triangulation(points)
+    >>> b = gt.betweenness(g)
+    >>> gt.graph_draw(g, pos=pos, pin=True, size=(8,8), vsize=0.08, vcolor=b[0],
+    ...               ecolor=b[1], output="triang.png")
+    <...>
+    >>> g, pos = gt.triangulation(points, type="delaunay")
+    >>> b = gt.betweenness(g)
+    >>> gt.graph_draw(g, pos=pos, pin=True, size=(8,8), vsize=0.08, vcolor=b[0],
+    ...               ecolor=b[1], output="triang-delaunay.png")
+    <...>
+
+    2D triangulation of random points:
+
+    .. image:: triang.png
+    .. image:: triang-delaunay.png
+
+    *Left:* Simple triangulation. *Right:* Delaunay triangulation. The colors
+    correspond to the betweeness centrality.
+
+    References
+    ----------
+    .. [cgal_triang] http://www.cgal.org/Manual/last/doc_html/cgal_manual/Triangulation_3/Chapter_main.html
+
+    """
+
+    if points.shape[1] not in [2,3]:
+        raise ValueError("points array must have shape N x d, with d either 2 or 3.")
+    # copy points to ensure continuity and correct data type
+    points = numpy.array(points, dtype='float64')
+    if points.shape[1] == 2:
+        npoints = numpy.zeros((points.shape[0], 3))
+        npoints[:,:2] = points
+        points = npoints
+    g = Graph(directed=False)
+    pos = g.new_vertex_property("vector<double>")
+    libgraph_tool_generation.triangulation(g._Graph__graph, points,
+                                           _prop("v", g, pos), type)
+    return g, pos
+