Merge pull request #7985 from andfoy/add_delaunay_2d

Add 2D Delaunay triangulation

cupyx/scipy/spatial/__init__.py

@@ -1 +1,2 @@
 from cupyx.scipy.spatial.distance import distance_matrix   # NOQA
+from cupyx.scipy.spatial._delaunay import Delaunay  # NOQA

cupyx/scipy/spatial/_delaunay.py

@@ -0,0 +1,158 @@
+
+import cupy
+from cupyx.scipy.spatial.delaunay_2d._tri import GDel2D
+
+
+class Delaunay:
+    """
+    Delaunay(points, furthest_site=False, incremental=False)
+
+    Delaunay tessellation in 2 dimensions.
+
+    Parameters
+    ----------
+    points : ndarray of floats, shape (npoints, ndim)
+        Coordinates of points to triangulate
+    furthest_site : bool, optional
+        Whether to compute a furthest-site Delaunay triangulation. This option
+        will be ignored, since it is not supported by CuPy
+        Default: False
+    incremental : bool, optional
+        Allow adding new points incrementally. This takes up some additional
+        resources. This option will be ignored, since it is not supported by
+        CuPy. Default: False
+
+    Attributes
+    ----------
+    points : ndarray of double, shape (npoints, ndim)
+        Coordinates of input points.
+    simplices : ndarray of ints, shape (nsimplex, ndim+1)
+        Indices of the points forming the simplices in the triangulation.
+        For 2-D, the points are oriented counterclockwise.
+    neighbors : ndarray of ints, shape (nsimplex, ndim+1)
+        Indices of neighbor simplices for each simplex.
+        The kth neighbor is opposite to the kth vertex.
+        For simplices at the boundary, -1 denotes no neighbor.0
+    vertex_to_simplex : ndarray of int, shape (npoints,)
+        Lookup array, from a vertex, to some simplex which it is a part of.
+        If qhull option "Qc" was not specified, the list will contain -1
+        for points that are not vertices of the tessellation.
+    convex_hull : ndarray of int, shape (nfaces, ndim)
+        Vertices of facets forming the convex hull of the point set.
+        The array contains the indices of the points belonging to
+        the (N-1)-dimensional facets that form the convex hull
+        of the triangulation.
+
+        .. note::
+
+           Computing convex hulls via the Delaunay triangulation is
+           inefficient and subject to increased numerical instability.
+           Use `ConvexHull` instead.
+    coplanar : ndarray of int, shape (ncoplanar, 3)
+        Indices of coplanar points and the corresponding indices of
+        the nearest facet and the nearest vertex.  Coplanar
+        points are input points which were *not* included in the
+        triangulation due to numerical precision issues.
+
+        If option "Qc" is not specified, this list is not computed.
+
+        .. versionadded:: 0.12.0
+    vertex_neighbor_vertices : tuple of two ndarrays of int; (indptr, indices)
+        Neighboring vertices of vertices. The indices of neighboring
+        vertices of vertex `k` are ``indices[indptr[k]:indptr[k+1]]``.
+    furthest_site
+        True if this was a furthest site triangulation and False if not.
+
+    Notes
+    -----
+    This implementation makes use of GDel2D to perform the triangulation in 2D.
+    See [1]_ for more information.
+
+    References
+    ----------
+    .. [1] A GPU accelerated algorithm for 3D Delaunay triangulation (2014).
+    Thanh-Tung Cao, Ashwin Nanjappa, Mingcen Gao, Tiow-Seng Tan.
+    Proc. 18th ACM SIGGRAPH Symp. Interactive 3D Graphics and Games, 47-55.
+    """
+
+    def __init__(self, points, furthest_site=False, incremental=False):
+        if points.shape[-1] != 2:
+            raise ValueError('Delaunay only supports 2D inputs at the moment.')
+
+        if furthest_site:
+            raise ValueError(
+                'furthest_site argument is not supported by CuPy.')
+
+        if incremental:
+            raise ValueError(
+                'incremental argument is not supported by CuPy.')
+
+        self.points = cupy.array(points, cupy.float64, copy=True)
+        # self._points = cupy.unique(self._points, axis=0)
+        self._triangulator = GDel2D(self.points)
+        self.simplices, self.neighbors = self._triangulator.compute()
+
+    def _find_simplex_coordinates(self, xi, eps, find_coords=False):
+        """
+        Find the simplices containing the given points.
+
+        Parameters
+        ----------
+        xi : ndarray of double, shape (..., ndim)
+            Points to locate
+        eps: float
+            Tolerance allowed in the inside-triangle check.
+        find_coords: bool, optional
+            Whether to return the barycentric coordinates of `xi`
+            with respect to the found simplices.
+
+        Returns
+        -------
+        i : ndarray of int, same shape as `xi`
+            Indices of simplices containing each point.
+            Points outside the triangulation get the value -1.
+        c : ndarray of float64, same shape as `xi`, optional
+            Barycentric coordinates of `xi` with respect to the enclosing
+            simplices. Returned only when `find_coords` is True.
+        """
+        out = cupy.empty((xi.shape[0],), dtype=cupy.int32)
+        c = cupy.empty(0, dtype=cupy.float64)
+        if find_coords:
+            c = cupy.empty((xi.shape[0], xi.shape[-1] + 1), dtype=cupy.float64)
+
+        out, c = self._triangulator.find_point_in_triangulation(
+            xi, eps, find_coords)
+
+        if find_coords:
+            return out, c
+        return out
+
+    def find_simplex(self, xi, bruteforce=False, tol=None):
+        """
+        Find the simplices containing the given points.
+
+        Parameters
+        ----------
+        xi : ndarray of double, shape (..., ndim)
+            Points to locate
+        bruteforce : bool, optional
+            Whether to only perform a brute-force search. Not used by CuPy
+        tol : float, optional
+            Tolerance allowed in the inside-triangle check.
+            Default is ``100*eps``.
+
+        Returns
+        -------
+        i : ndarray of int, same shape as `xi`
+            Indices of simplices containing each point.
+            Points outside the triangulation get the value -1.
+        """
+        if tol is None:
+            eps = 100 * cupy.finfo(cupy.double).eps
+        else:
+            eps = float(tol)
+
+        if xi.shape[-1] != 2:
+            raise ValueError('Delaunay only supports 2D inputs at the moment.')
+
+        return self._find_simplex_coordinates(xi, eps)