Add method to get convex hull from a concave polygon

Author: AntSimi

setup.py

@@ -9,6 +9,7 @@
 
 setup(
     name="pyEddyTracker",
+    python_requires=">=3.7",
     version=versioneer.get_version(),
     cmdclass=versioneer.get_cmdclass(),
     description="Py-Eddy-Tracker libraries",

src/py_eddy_tracker/observations/observation.py

@@ -58,6 +58,7 @@
     close_center,
     get_pixel_in_regular,
     winding_number_poly,
+    convexs,
 )
 
 logger = logging.getLogger("pet")
@@ -101,9 +102,7 @@ def shifted_ellipsoid_degrees_mask2(lon0, lat0, lon1, lat1, minor=1.5, major=1.5
 
 class EddiesObservations(object):
     """
-    Class to hold eddy properties *amplitude* and counts of
-    *local maxima/minima* within a closed region of a sea level anomaly field.
-
+    Class to store eddy observations.
     """
 
     __slots__ = (
@@ -1606,6 +1605,17 @@ def last_obs(self):
         m[:-1][self["n"][1:] == 0] = True
         return self.extract_with_mask(m)
 
+    def is_convex(self, intern=False):
+        """
+        Get flag of eddy convexity
+
+        :param bool intern: If true use speed contour instead of effective contour
+        :return: true if contour is convex
+        :rtype: array[bool]
+        """
+        xname, yname = self.intern(intern)
+        return convexs(self[xname], self[yname])
+
     def inside(self, x, y, intern=False):
         """
         True for each postion inside an eddy
@@ -1769,9 +1779,9 @@ def insidepoly(x_p, y_p, x_c, y_c):
         x_c_max, y_c_max = x_c[i].max(), y_c[i].max()
         v = create_vertice(x_c[i], y_c[i])
         for j in range(nb_p):
-            x, y = x_p[j], y_p[j]
             if flag[j]:
                 continue
+            x, y = x_p[j], y_p[j]
             if x > x_c_min and x < x_c_max and y > y_c_min and y < y_c_max:
                 if winding_number_poly(x, y, v) != 0:
                     flag[j] = True

src/py_eddy_tracker/poly.py

@@ -88,6 +88,101 @@ def poly_area(x, y):
     return abs(p_area) * 0.5
 
 
+@njit(cache=True)
+def convexs(x, y):
+    """
+    Check if polygons are convex
+
+    :param array[float] x:
+    :param array[float] y:
+    :return: True if convex
+    :rtype: array[bool]
+    """
+    nb_poly = x.shape[0]
+    flag = empty(nb_poly, dtype=numba_types.bool_)
+    for i in range(nb_poly):
+        flag[i] = convex(x[i], y[i])
+    return flag
+
+
+@njit(cache=True)
+def convex(x, y):
+    """
+    Check if polygon is convex
+
+    :param array[float] x:
+    :param array[float] y:
+    :return: True if convex
+    :rtype: bool
+    """
+    nb = x.shape[0]
+    x0, y0, x1, y1, x2, y2 = x[-2], y[-2], x[-1], y[-1], x[1], y[1]
+    # if first is left it must be always left if it's right it must be always right
+    ref = is_left(x0, y0, x1, y1, x2, y2)
+    # We skip 0 because it's same than -1
+    # We skip 1 because we tested previously
+    for i in range(2, nb):
+        # shift position
+        x0, y0, x1, y1 = x1, y1, x2, y2
+        x2, y2 = x[i], y[i]
+        # test
+        if ref != is_left(x0, y0, x1, y1, x2, y2):
+            return False
+    return True
+
+
+@njit(cache=True)
+def get_convex_hull(x, y):
+    """
+    Get convex polygon which enclosed current polygon
+
+    Work only if contour is describe anti-clockwise
+
+    :param array[float] x:
+    :param array[float] y:
+    :return: a convex polygon
+    :rtype: array,array
+    """
+    nb = x.shape[0] - 1
+    indices = list()
+    # leftmost point
+    i_first = x[:-1].argmin()
+    indices.append(i_first)
+    i_next = (i_first + 1) % nb
+    # Will define bounds line
+    x0, y0, x1, y1 = x[i_first], y[i_first], x[i_next], y[i_next]
+    xf, yf = x0, y0
+    # we will check if no point are right
+    while True:
+        i_test = (i_next + 1) % nb
+        # value to test
+        xt, yt = x[i_test], y[i_test]
+        # We will test all the position until we touch first one,
+        # If all next position are on the left we keep x1, y1
+        # if not we will replace by xt,yt which are more outter
+        while is_left(x0, y0, x1, y1, xt, yt):
+            i_test += 1
+            i_test %= nb
+            if i_test == i_first:
+                x0, y0 = x1, y1
+                indices.append(i_next)
+                i_next += 1
+                i_next %= nb
+                x1, y1 = x[i_next], y[i_next]
+                break
+            xt, yt = x[i_test], y[i_test]
+        if i_test != i_first:
+            i_next = i_test
+            x1, y1 = x[i_next], y[i_next]
+        if i_next == (i_first - 1) % nb:
+            if is_left(x0, y0, x1, y1, xf, yf):
+                indices.append(i_next)
+            break
+    indices.append(i_first)
+    indices = array(indices)
+    return x[indices], y[indices]
+
+
 @njit(cache=True)
 def winding_number_poly(x, y, xy_poly):
     """
@@ -155,7 +250,7 @@ def winding_number_grid_in_poly(x_1d, y_1d, i_x0, i_x1, x_size, i_y0, xy_poly):
 
 
 @njit(cache=True, fastmath=True)
-def close_center(x0, y0, x1, y1, delta=.1):
+def close_center(x0, y0, x1, y1, delta=0.1):
     """
     Compute an overlap with circle parameter and return a percentage
 

tests/test_poly.py

@@ -1,9 +1,10 @@
-from py_eddy_tracker.poly import poly_area_vertice, fit_circle
+from py_eddy_tracker.poly import poly_area_vertice, fit_circle, convex, get_convex_hull
 from numpy import array, pi
 from pytest import approx
 
 # Vertices for next test
 V = array(((2, 2, 3, 3, 2), (-10, -9, -9, -10, -10)))
+V_concave = array(((2, 2, 2.5, 3, 3, 2), (-10, -9, -9.5, -9, -10, -10)))
 
 
 def test_poly_area():
@@ -16,3 +17,14 @@ def test_fit_circle():
     assert y0 == approx(-9.5, rel=1e-10)
     assert r == approx(2 ** 0.5 / 2, rel=1e-10)
     assert err == approx((1 - 2 / pi) * 100, rel=1e-10)
+
+
+def test_convex():
+    assert convex(*V) is True
+    assert convex(*V[::-1]) is True
+    assert convex(*V_concave) is False
+    assert convex(*V_concave[::-1]) is False
+
+
+def test_convex_hull():
+    assert convex(*get_convex_hull(*V_concave)) is True