error is now compute outside of fitcircle

Author: AntSimi

src/py_eddy_tracker/dataset/grid.py

@@ -50,7 +50,6 @@
 from ..generic import (
     distance,
     interp2d_geo,
-    fit_circle,
     uniform_resample,
     coordinates_to_local,
     local_to_coordinates,
@@ -63,6 +62,7 @@
     winding_number_poly,
     create_vertice,
     poly_area,
+    fit_circle,
 )
 
 logger = logging.getLogger("pet")
@@ -153,7 +153,7 @@ def _circle_from_equal_area(vertice):
         # logger.debug('%d coordinates %s,%s', len(lons),lons,
         # lats)
         return 0, -90, nan, nan
-    return lon0, lat0, (poly_area(create_vertice(c_x, c_y)) / pi) ** 0.5, nan
+    return lon0, lat0, (poly_area(c_x, c_y) / pi) ** 0.5, nan
 
 
 @njit(cache=True, fastmath=True)

src/py_eddy_tracker/generic.py

@@ -41,7 +41,6 @@
     concatenate,
 )
 from numba import njit, prange, types as numba_types
-from numpy.linalg import lstsq
 from .poly import winding_number_grid_in_poly
 
 
@@ -195,79 +194,6 @@ def custom_convolution(data, mask, kernel):
     return out
 
 
-@njit(cache=True)
-def fit_circle(x_vec, y_vec):
-    nb_elt = x_vec.shape[0]
-    p_inon_x = empty(nb_elt)
-    p_inon_y = empty(nb_elt)
-
-    # last coordinates == first
-    x_mean = x_vec[1:].mean()
-    y_mean = y_vec[1:].mean()
-
-    norme = (x_vec[1:] - x_mean) ** 2 + (y_vec[1:] - y_mean) ** 2
-    norme_max = norme.max()
-    scale = norme_max ** 0.5
-
-    # Form matrix equation and solve it
-    # Maybe put f4
-    datas = ones((nb_elt - 1, 3))
-    datas[:, 0] = 2.0 * (x_vec[1:] - x_mean) / scale
-    datas[:, 1] = 2.0 * (y_vec[1:] - y_mean) / scale
-
-    (center_x, center_y, radius), _, _, _ = lstsq(datas, norme / norme_max)
-
-    # Unscale data and get circle variables
-    radius += center_x ** 2 + center_y ** 2
-    radius **= 0.5
-    center_x *= scale
-    center_y *= scale
-    # radius of fitted circle
-    radius *= scale
-    # center X-position of fitted circle
-    center_x += x_mean
-    # center Y-position of fitted circle
-    center_y += y_mean
-
-    # area of fitted circle
-    c_area = (radius ** 2) * pi
-    # Find distance between circle center and contour points_inside_poly
-    for i_elt in range(nb_elt):
-        # Find distance between circle center and contour points_inside_poly
-        dist_poly = (
-            (x_vec[i_elt] - center_x) ** 2 + (y_vec[i_elt] - center_y) ** 2
-        ) ** 0.5
-        # Indices of polygon points outside circle
-        # p_inon_? : polygon x or y points inside & on the circle
-        if dist_poly > radius:
-            p_inon_y[i_elt] = center_y + radius * (y_vec[i_elt] - center_y) / dist_poly
-            p_inon_x[i_elt] = center_x - (center_x - x_vec[i_elt]) * (
-                center_y - p_inon_y[i_elt]
-            ) / (center_y - y_vec[i_elt])
-        else:
-            p_inon_x[i_elt] = x_vec[i_elt]
-            p_inon_y[i_elt] = y_vec[i_elt]
-
-    # Area of closed contour/polygon enclosed by the circle
-    p_area_incirc = 0
-    p_area = 0
-    for i_elt in range(nb_elt - 1):
-        # Indices of polygon points outside circle
-        # p_inon_? : polygon x or y points inside & on the circle
-        p_area_incirc += (
-            p_inon_x[i_elt] * p_inon_y[1 + i_elt]
-            - p_inon_x[i_elt + 1] * p_inon_y[i_elt]
-        )
-        # Shape test
-        # Area and centroid of closed contour/polygon
-        p_area += x_vec[i_elt] * y_vec[1 + i_elt] - x_vec[1 + i_elt] * y_vec[i_elt]
-    p_area = abs(p_area) * 0.5
-    p_area_incirc = abs(p_area_incirc) * 0.5
-
-    a_err = (c_area - 2 * p_area_incirc + p_area) * 100.0 / c_area
-    return center_x, center_y, radius, a_err
-
-
 @njit(cache=True, fastmath=True)
 def uniform_resample(x_val, y_val, num_fac=2, fixed_size=None):
     """

src/py_eddy_tracker/poly.py

@@ -21,7 +21,8 @@
 ===========================================================================
 """
 
-from numpy import empty, where, array
+from numpy import empty, where, array, ones, pi
+from numpy.linalg import lstsq
 from numba import njit, prange, types as numba_types
 from Polygon import Polygon
 
@@ -74,16 +75,27 @@ def poly_contain_poly(xy_poly_out, xy_poly_in):
 
 
 @njit(cache=True)
-def poly_area(vertice):
+def poly_area_vertice(v):
     """
-    :param vertice vertice: polygon vertice
+    :param vertice v: polygon vertice
     :return: area of polygon in coordinates unit
     :rtype: float
     """
-    p_area = vertice[0, 0] * (vertice[1, 1] - vertice[1, -2])
-    nb_elt = vertice.shape[0]
-    for i_elt in range(1, nb_elt - 1):
-        p_area += vertice[i_elt, 0] * (vertice[1 + i_elt, 1] - vertice[i_elt - 1, 1])
+    return poly_area(v[:, 0], v[:, 1])
+
+
+@njit(cache=True)
+def poly_area(x, y):
+    """
+    :param array x:
+    :param array y:
+    :return: area of polygon in coordinates unit
+    :rtype: float
+    """
+    p_area = x[0] * (y[1] - y[-2])
+    nb = x.shape[0]
+    for i in range(1, nb - 1):
+        p_area += x[i] * (y[1 + i] - y[i - 1])
     return abs(p_area) * 0.5
 
 
@@ -295,3 +307,81 @@ def polygon_overlap(p0, p1, minimal_area=False):
         else:
             cost[i] = intersection / (p0 + p_).area()
     return cost
+
+
+@njit(cache=True)
+def fit_circle(x, y):
+    """
+    From a polygon, function will fit a circle
+    Must be call with local coordinates (in m, to get a radius in m)
+
+    :param array x: x of polygon
+    :param array y: y of polygon
+    :return: x0, y0, radius, shape_error
+    :rtype: (float,float,float,float)
+    """
+    nb_elt = x.shape[0]
+
+    # last coordinates == first
+    x_mean = x[1:].mean()
+    y_mean = y[1:].mean()
+
+    norme = (x[1:] - x_mean) ** 2 + (y[1:] - y_mean) ** 2
+    norme_max = norme.max()
+    scale = norme_max ** 0.5
+
+    # Form matrix equation and solve it
+    # Maybe put f4
+    datas = ones((nb_elt - 1, 3))
+    datas[:, 0] = 2.0 * (x[1:] - x_mean) / scale
+    datas[:, 1] = 2.0 * (y[1:] - y_mean) / scale
+
+    (x0, y0, radius), _, _, _ = lstsq(datas, norme / norme_max)
+
+    # Unscale data and get circle variables
+    radius += x0 ** 2 + y0 ** 2
+    radius **= 0.5
+    x0 *= scale
+    y0 *= scale
+    # radius of fitted circle
+    radius *= scale
+    # center X-position of fitted circle
+    x0 += x_mean
+    # center Y-position of fitted circle
+    y0 += y_mean
+
+    err = shape_error(x, y, x0, y0, radius)
+    return x0, y0, radius, err
+
+
+@njit(cache=True, fastmath=True)
+def shape_error(x, y, x0, y0, r):
+    """
+    With a polygon(x,y) in local coordinates
+    and circle properties(x0, y0, r), function compute a shape error:
+
+    .. math:: ShapeError = \\frac{Polygon_{area} + Circle_{area} - 2 * Intersection_{area}}{Circle_{area}} * 100
+
+    When error > 100, area of difference is bigger than circle area
+
+    :param array x: x of polygon
+    :param array y: y of polygon
+    :param float x0: x center of circle
+    :param float y0: y center of circle
+    :param float r: radius of circle
+    :return: shape error
+    :rtype: float
+    """
+    # circle area
+    c_area = (r ** 2) * pi
+    p_area = poly_area(x, y)
+    nb = x.shape[0]
+    x, y = x.copy(), y.copy()
+    # Find distance between circle center and polygon
+    for i in range(nb):
+        dx, dy = x[i] - x0, y[i] - y0
+        rd = r / (dx ** 2 + dy ** 2) ** 0.5
+        if rd < 1:
+            x[i] = x0 + dx * rd
+            y[i] = y0 + dy * rd
+    return 100 + (p_area - 2 * poly_area(x, y)) / c_area * 100

tests/test_poly.py

@@ -0,0 +1,18 @@
+from py_eddy_tracker.poly import poly_area_vertice, fit_circle
+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)))
+
+
+def test_poly_area():
+    assert 1 == poly_area_vertice(V.T)
+
+
+def test_fit_circle():
+    x0, y0, r, err = fit_circle(*V)
+    assert x0 == approx(2.5, rel=1e-10)
+    assert y0 == approx(-9.5, rel=1e-10)
+    assert r == approx(2 ** 0.5 / 2, rel=1e-10)
+    assert err == approx((r ** 2 * pi - 1) / (r ** 2 * pi) * 100, rel=1e-10)