Replace nearest_path func by nearest_path which bbox which contain pt to speed up process

Author: adelepoulle

src/py_eddy_tracker/property_functions.py

@@ -237,7 +237,8 @@ def get_uavg(eddy, contours, centlon_e, centlat_e, poly_eff, grd,
             corrected_coll_index = nb_coll - coll_index - 1
 
         # Leave loop if no contours at level citer.index
-        theindex = eddy.swirl.get_index_nearest_path(
+        #~ theindex = eddy.swirl.get_index_nearest_path(
+        theindex = eddy.swirl.get_index_nearest_path_bbox_contain_pt(
             corrected_coll_index, centlon_e, centlat_e)
         if theindex is None:
             continue

src/py_eddy_tracker/property_objects.py

@@ -33,7 +33,8 @@
 from scipy.ndimage import binary_erosion
 from scipy.ndimage import minimum_filter
 from numpy import array, isfinite, ma, where, ones
-from .tools import index_from_nearest_path, distance_matrix
+from .tools import index_from_nearest_path, \
+    index_from_nearest_path_with_pt_in_bbox, distance_matrix
 import logging
 
 
@@ -224,28 +225,61 @@ class SwirlSpeed(object):
       grd:
         A grid object
     """
-    def __init__(self, contour):
+    def __init__(self, contour, nearest_contain_in_bbox=False):
         """
         c_i : index to contours
         l_i : index to levels
         """
         x_list, y_list, ci_list, li_list = [], [], [], []
+        x_min_list, y_min_list, x_max_list, y_max_list = [], [], [], []
 
         for cont in contour.collections:
             for coll in cont.get_paths():
-                x_list.append(coll.vertices[:, 0])
-                y_list.append(coll.vertices[:, 1])
-                ci_list.append(len(coll.vertices[:, 0]))
+                x_val, y_val = coll.vertices[:, 0], coll.vertices[:, 1]
+                x_min_list.append(x_val.min())
+                x_max_list.append(x_val.max())
+                y_min_list.append(y_val.min())
+                y_max_list.append(y_val.max())
+                x_list.append(x_val)
+                y_list.append(y_val)
+                ci_list.append(coll.vertices.shape[0])
             li_list.append(len(cont.get_paths()))
 
         self.x_value = array([val for sublist in x_list for val in sublist])
         self.y_value = array([val for sublist in y_list for val in sublist])
+
+        self.x_min_per_c = array(x_min_list)
+        self.y_min_per_c = array(y_min_list)
+        self.x_max_per_c = array(x_max_list)
+        self.y_max_per_c = array(y_max_list)
+
         self.nb_pt_per_c = array(ci_list, dtype='u4')
         self.c_i = array(self.nb_pt_per_c.cumsum() - self.nb_pt_per_c,
                          dtype='u4')
         self.nb_c_per_l = array(li_list, dtype='u4')
         self.l_i = array(self.nb_c_per_l.cumsum() - self.nb_c_per_l,
                          dtype='u4')
+                         
+        self.nearest_contain_in_bbox = nearest_contain_in_bbox
+
+    def get_index_nearest_path_bbox_contain_pt(self, level, xpt, ypt):
+        """
+        """
+        return index_from_nearest_path_with_pt_in_bbox(
+            level,
+            self.l_i,
+            self.nb_c_per_l,
+            self.nb_pt_per_c,
+            self.c_i,
+            self.x_value,
+            self.y_value,
+            self.x_min_per_c,
+            self.y_min_per_c,
+            self.x_max_per_c,
+            self.y_max_per_c,
+            xpt,
+            ypt
+            )
 
     def get_index_nearest_path(self, level, xpt, ypt):
         """

src/py_eddy_tracker/tools.pyx

@@ -15,6 +15,7 @@ cdef DTYPE_coord D2R = 0.017453292519943295
 cdef DTYPE_coord PI = 3.141592653589793
 cdef DTYPE_coord EARTH_DIAMETER = 6371.3150 * 2
 
+
 @wraparound(False)
 @boundscheck(False)
 def fit_circle_c(
@@ -200,6 +201,7 @@ def winding_number_poly(
                     wn -= 1
     return wn
 
+
 @wraparound(False)
 @boundscheck(False)
 def distance(
@@ -284,9 +286,8 @@ cdef dist_array_size(
         ):
     """Give slice to select data
     """
-    cdef DTYPE_ui i_elt, nb_pts
+    cdef DTYPE_ui i_elt
     i_end[0] = i_start + nb_c_per_l
-    nb_pts = 0
 
     c_start[0] = c_i[i_start]
     c_end[0] = c_start[0]
@@ -370,3 +371,53 @@ cdef nearest_contour_index(
             dist_ref = dist
             i_ref = i_elt
     nearesti[0] = i_ref - start
+
+
+@wraparound(False)
+@boundscheck(False)
+def index_from_nearest_path_with_pt_in_bbox(
+        DTYPE_ui level_index,
+        ndarray[DTYPE_ui] l_i,
+        ndarray[DTYPE_ui] nb_c_per_l,
+        ndarray[DTYPE_ui] nb_pt_per_c,
+        ndarray[DTYPE_ui] indices_of_first_pts,
+        ndarray[DTYPE_coord] x_value,
+        ndarray[DTYPE_coord] y_value,
+        ndarray[DTYPE_coord] x_min_per_c,
+        ndarray[DTYPE_coord] y_min_per_c,
+        ndarray[DTYPE_coord] x_max_per_c,
+        ndarray[DTYPE_coord] y_max_per_c,
+        DTYPE_coord xpt,
+        DTYPE_coord ypt,
+        ):
+    """Get index from nearest path in edge bbox contain pt
+    """
+    cdef DTYPE_ui i_start_c, i_end_c, i_elt_c, i_start_pt, i_end_pt, i_elt_pt
+    cdef DTYPE_ui i_ref, nb_contour, find_contour
+    cdef DTYPE_coord dist_ref, dist
+    find_contour = 0
+    i_start_c = l_i[level_index]
+    nb_contour = nb_c_per_l[level_index]
+    if nb_contour == 0:
+        return None
+    i_end_c = i_start_c + nb_c_per_l[level_index]
+
+    i_ref = i_start_c    
+    i_start_pt = indices_of_first_pts[i_start_c]
+    dist_ref = (x_value[i_start_pt] - xpt) ** 2 + (y_value[i_start_pt] - ypt) ** 2
+
+    for i_elt_c from i_start_c <= i_elt_c < i_end_c:
+        if x_min_per_c[i_elt_c] > xpt or x_max_per_c[i_elt_c] < xpt or y_min_per_c[i_elt_c] > ypt or y_max_per_c[i_elt_c] < ypt:
+            continue
+        i_start_pt = indices_of_first_pts[i_elt_c]
+        i_end_pt = i_start_pt + nb_pt_per_c[i_elt_c]
+        find_contour = 1
+        
+        for i_elt_pt from i_start_pt <= i_elt_pt < i_end_pt:
+            dist = (x_value[i_elt_pt] - xpt) ** 2 + (y_value[i_elt_pt] - ypt) ** 2
+            if dist < dist_ref:
+                dist_ref = dist
+                i_ref = i_elt_c            
+    if find_contour == 0:
+        return None
+    return i_ref - i_start_c