WIP: Making a Python package using Cython to call the CPP

.gitignore

@@ -1,3 +1,7 @@
+.DS_Store
+concaveman.egg-info
+build/
+
 # Prerequisites
 *.d
 

demos/demos.py

@@ -0,0 +1,59 @@
+#!/usr/bin/env python
+
+"""
+demos / checks for concaveman
+"""
+
+import numpy as np
+
+import matplotlib.pyplot as plt
+
+from concaveman import ConcaveHull
+
+def circular_points(N):
+    """
+    create a bunch of random points in a circle
+    """
+    # create a bunch of random points
+    points_in_circle = []
+
+    while len(points_in_circle) < N:
+        n = max(100, N // 2)
+#        print(f"picking {n} random points")
+        points = (np.random.random((n, 2)) - 0.5) * 100
+
+        # check which ones are in the circle
+        points_in_circle.extend(points[np.hypot(points[:,0], points[:,1]) <= 50.0])
+        # print(f"there are now {len(points_in_circle)} in the circle")
+    return np.array(points_in_circle[:N])
+
+# get some points
+points = circular_points(100)
+
+# compute the hulls
+hull = ConcaveHull(points, concavity=1.0)
+hull_1 = hull.cc_hull
+
+hull.concavity = 3.0
+hull.compute()
+hull_2 = hull.cc_hull
+
+
+hull.concavity = np.inf
+hull.compute()
+hull_inf = hull.cc_hull
+
+
+fig, ax = plt.subplots()
+ax.plot(points[:, 0], points[:, 1], '.')
+ax.plot(hull_1[:, 0], hull_1[:, 1], label="alpha=1.0")
+ax.plot(hull_2[:, 0], hull_2[:, 1], label="alpha=3.0")
+ax.plot(hull_inf[:, 0], hull_inf[:, 1], label="alpha=Inf")
+ax.legend()
+
+ax.set_xlim(-80, 80)
+ax.set_ylim(-80, 80)
+ax.set_aspect('equal', 'box')
+
+plt.show()
+

setup.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python
+
+from setuptools import setup, Extension
+
+from Cython.Build import cythonize
+
+import os
+import numpy as np  # For the include directory.
+
+from pathlib import Path
+
+rootpath = Path(__file__).parent.resolve()
+
+long_description = open(rootpath / 'README.md').read()
+
+include_dirs = [np.get_include(),
+                str(rootpath / 'src' / 'cpp')]
+
+# # Need the stdint header for Windows (VS2008).
+# if sys.platform.startswith('win') and sys.version_info.major <= 2:
+#     include_dirs.append(os.path.join('rootpath', 'src', 'win_headers'))
+
+
+ext_modules = cythonize([Extension("concaveman.build_hull",
+                                   ["src/concaveman/build_hull.pyx",
+                                    "src/cpp/concaveman.cpp"],
+                                   include_dirs=include_dirs,
+                                   language="c++",
+                                   # this work on the mac, should work in Linux
+                                   # for Windows: who knows?
+                                   extra_compile_args=['-std=c++11'],
+                                   )])
+
+
+def extract_version():
+    fname = os.path.join(rootpath, 'src', 'concaveman', '__init__.py')
+    with open(fname) as f:
+        for line in f:
+            if (line.startswith('__version__')):
+                version = line.split('=')[1].strip().strip('"')
+                break
+        else:
+            raise ValueError("Couldn't find __version__ in %s" % fname)
+    return version
+
+setup(
+    name="concaveman",
+    version=extract_version(),
+    description="Python wrappers around a C++ concave hull Implementation",
+    long_description=long_description,
+    author="s.adaszewski, Christopher H. Barker",
+    author_email="Chris.Barker@noaa.gov",
+    url="https://github.com/NOAA-ORR-ERD",
+    license="BSD",
+    # keywords = "",
+    ext_modules=ext_modules,
+    packages=["concaveman", "concaveman/tests"],
+    package_dir={'': 'src'},
+    install_requires=['numpy', 'scipy'],
+    setup_requires=['cython>0.29'],
+    tests_require=['pytest'],
+    classifiers=[
+        "Development Status :: 4 - Beta",
+        "License :: Public Domain",
+        "Intended Audience :: Developers",
+        "Intended Audience :: Science/Research",
+        "Operating System :: OS Independent",
+        "Programming Language :: C++",
+        "Programming Language :: Cython",
+        "Programming Language :: Python :: 3 :: Only",
+        "Programming Language :: Python :: Implementation :: CPython",
+        "Topic :: Utilities",
+        "Topic :: Scientific/Engineering",
+        "Topic :: Scientific/Engineering :: Visualization",
+    ],
+)

src/concaveman/__init__.py

@@ -0,0 +1,16 @@
+
+__version__ = "0.1.0"
+
+import numpy as np
+
+from .concaveman import ConcaveHull
+
+
+# an example square with a couple interior points
+EXAMPLE_SQUARE = np.array(((5.0, 10.0),
+                           (5.0, 20.0),
+                           (10.0, 20.0),
+                           (10.0, 10.0),
+                           (6.0, 13.0),
+                           (8.0, 18.0),
+                           ))

src/concaveman/build_hull.pyx

@@ -0,0 +1,106 @@
+# cython: language_level=3
+
+"""
+wrapper of concaveman
+
+NOTE: after working on this for a while, I realized that concaveman.cpp was
+mostly a wrapper that was set up to call the templates from C, with
+regular C arrays.
+
+But we should be able to call the Templates directly from Cython, so save that
+translation step.
+
+but other than some wasted data copying, this works.
+
+NOTE: I removed the function pointer for the free() call, as this code is now using that
+      malloc'ed pointer to build the returning numpy array -- so should be freed by
+      python when the numpy array is deleted.
+"""
+
+import numpy as np
+cimport numpy as cnp
+
+from libc.stdlib cimport free
+from libc.string cimport memcpy
+
+
+cdef extern from "pyconcaveman.h":
+
+    void pyconcaveman2d(double *points_c,
+                        size_t num_points,
+                        int *hull_points_c,
+                        size_t num_hull_points,
+                        double concavity,
+                        double lengthThreshold,
+                        double **concave_points_c,
+                        size_t *num_concave_points,
+                        )
+
+cpdef concave_hull(cnp.float64_t[:,:] points,
+                   int[:] hull,
+                   double concavity=2.0,
+                   double length_threshold=0.0):
+    """
+    compute the concave hull of a bunch of points
+
+    :param points: The points to make the hull out of.
+    :type points: Nx2 array of float64 type
+
+    :param hull: The convex hull of the points
+    :type points: 1D array of int32
+
+    :param concavity=2.0: concavity parameter: large means smoother
+    :type concavity: python float
+
+    :param length_threshold:
+    :type length_threshold: python float
+    """
+
+    if points.shape[1] != 2:
+        raise ValueError('points must be an Nx2 array')
+
+    # now a memoryview, so any isn't helpful here.
+    # should this check be at a higher level?
+    # if np.any(hull >= len(points)) or np.any(hull < 0):
+    #     raise ValueError('hull indices out of bounds')
+
+    cdef cnp.float64_t* p_concave_points = NULL
+    cdef size_t num_concave_points = 0
+
+    # print("in cython: about to call pyconcaveman2d")
+
+    pyconcaveman2d(&points[0, 0],
+                   len(points),
+                   &hull[0],
+                   len(hull),
+                   concavity,
+                   length_threshold,
+                   &p_concave_points,
+                   &num_concave_points,
+                   )
+
+    # print("cpp concave hull returned")
+    # print("num concave points:", num_concave_points)
+
+    # create an array to return
+    cdef cnp.ndarray[cnp.float64_t, ndim=2, mode="c"] arr_concave_points
+    arr_concave_points = np.zeros((num_concave_points, 2), dtype=np.float64)
+
+    # loping through element by element
+    # cdef cnp.float64_t* temp = &arr_concave_points[0, 0]
+    # for i in range(num_concave_points * 2):
+    #     temp[i] = p_concave_points[i]
+
+    # copy the data to the ndarray's memory block directly
+    memcpy(&arr_concave_points[0, 0],
+           p_concave_points,
+           sizeof(cnp.float64_t) * num_concave_points * 2)
+
+    # print('in cython again: concave_points:\n', arr_concave_points)
+
+    # As we are compiling the C++ code all together, the
+    # plain free() should work.
+    # on OS-X, it doesn't crash nor leak memory.
+    free(p_concave_points)
+
+    return arr_concave_points

src/concaveman/concaveman.py

@@ -0,0 +1,100 @@
+import numpy as np
+import scipy.spatial
+
+from .build_hull import concave_hull
+
+
+class ConcaveHull():
+    """
+    class modeled after the scipy.spatial.ConvexHull implementation
+
+    One of these days, we could add area and other computations.
+    """
+
+    def __init__(self, points, hull=None, concavity=2.0, length_threshold=0.0):
+
+        self.points = points
+
+        if hull is not None:
+            if len(hull.shape) != 1:
+                raise ValueError('hull must be a 1-D array')
+
+            if np.any(hull >= len(points)) or np.any(hull < 0):
+                raise ValueError('hull indices out of bounds')
+            self.cvx_hull = hull
+            self.cvx_hull = scipy.spatial.ConvexHull(points).vertices
+
+        self.concavity = concavity
+        self.length_threshold = length_threshold
+
+        self.compute()
+
+    @property
+    def concavity(self):
+        return self._concavity
+    @concavity.setter
+    def concavity(self, value):
+        self._concavity = value
+        self._cc_hull = None
+
+    @property
+    def length_threshold(self):
+        return self._length_threshold
+    @length_threshold.setter
+    def length_threshold(self, value):
+        self._length_threshold = value
+        self._cc_hull = None
+
+    @property
+    def points(self):
+        return self._points
+    @points.setter
+    def points(self, points):
+        points = np.asarray(points, dtype=np.float64)
+        if len(points.shape) != 2 or points.shape[1] != 2:
+            raise ValueError('points must be an Nx2 array')
+        self._points = points
+        self.cvx_hull = None
+        self._cc_hull = None
+
+    @property
+    def cc_hull(self):
+        if self._cc_hull is None:
+            self.compute()
+        return self._cc_hull
+
+
+    def compute(self):
+        """
+        compute the concave hull -- should be called after changing any parameters
+        """
+        if self.cvx_hull is None:
+            self.cvx_hull = scipy.spatial.ConvexHull(self.points).vertices
+
+        self._cc_hull = concave_hull(self.points,
+                                     self.cvx_hull,
+                                     self.concavity,
+                                     self.length_threshold)
+
+
+
+
+
+# def concaveman2d(points, hull=None, concavity=2.0, length_threshold=0.0):
+#     points = np.asrray(points, dytpe=np.float64)
+
+#     if len(points.shape) != 2 or points.shape[0] != 2:
+#         raise ValueError('points must be an Nx2 array')
+
+#     if hull is not None:
+#         if len(hull.shape) != 1:
+#             raise ValueError('hull must be a 1-D array')
+
+#         if np.any(hull >= len(points)) or np.any(hull < 0):
+#             raise ValueError('hull indices out of bounds')
+#     else:
+#         # compute the convex hull
+#         # compute the ConvexHull
+#         hull = scipy.spatial.ConvexHull(points).vertices
+
+#     return concave_hull(points, hull, concavity, length_threshold)