diff --git a/floater/rclv.py b/floater/rclv.py index 158cd82..8dd5762 100644 --- a/floater/rclv.py +++ b/floater/rclv.py @@ -172,6 +172,57 @@ def contour_area(con): return region_area, hull_area, cd +def contour_CI(lx, ly, con, ji, region_area, dx, dy, border_i, border_j): + """Calculate the coherency index of a polygon contour. + + Parameters + ---------- + con : arraylike + A 2D array of vertices with shape (N,2) that follows the scikit + image conventions (con[:,0] are j indices) + ji : tuple + The index of the maximum in (j, i) order + lx : array_like + Array with shape (N,n,n) including the x position of particles at N time steps(position array is n*n, + and trajectories were unwrapped to correct the jump of displacement at the periodic boundaries) + ly : array_like + Array with shape (N,n,n) including the y position of particles at N time steps(position array is n*n, + and trajectories were unwrapped to correct the jump of displacement at the periodic boundaries) + region_area : float + Area of a polygon contour + dx, dy : float + Particle spacings in x and y directions, must have the same units as lx and ly, respectively. + + Returns + ------- + coherency index(CI) : float + """ + # the maximum + j, i = ji + + # get the position of the contour + con1 = con.copy() + con1[:, 0] += (j-border_j[0]) + con1[:, 1] += (i-border_i[0]) + + # label the points in eddy + mask = label_points_in_contours(lx[0].shape, [con1]) + + # get the positions of particles inside the eddy at each time + lx_p, ly_p = lx[:,mask==1], ly[:,mask==1] + + # calculate variance of particle positons at each time + var_p = np.var(lx_p,axis=-1) + np.var(ly_p,axis=-1) + + # calculate the minimum variance of particle postions + area1 = region_area*dx*dy # the real area of a polygon contour + var_min = area1/(2*np.pi) + + CI = (var_min - np.max(var_p))/var_min + + return CI + + def project_vertices(verts, lon0, lat0, dlon, dlat): """Project the logical coordinates of vertices into physical map coordiantes. @@ -207,7 +258,7 @@ def project_vertices(verts, lon0, lat0, dlon, dlat): def find_contour_around_maximum(data, ji, level, border_j=(5,5), - border_i=(5,5), max_footprint=None, proj_kwargs={}, + border_i=(5,5), max_width=100, max_footprint=None, proj_kwargs={}, periodic=(False, False)): j,i = ji max_val = data[j,i] @@ -269,15 +320,17 @@ def find_contour_around_maximum(data, ji, level, border_j=(5,5), if target_con is None and not ( grow_down or grow_up or grow_left or grow_right): raise ValueError("Couldn't find a contour") + if (np.array(border_i)>max_width).any() or (np.array(border_j)>max_width).any(): #set a limit on the width of the window + raise ValueError("Local region becomes too large.") return target_con, region_data, border_j, border_i -def convex_contour_around_maximum(data, ji, init_contour_step_frac=0.1, - border=5, - convex_def=0.01, convex_def_tol=0.001, +def convex_contour_around_maximum(data, lx, ly, ji, dx, dy, init_contour_step_frac=0.1, + border=5, max_width=100, + CI_th = -1.0, CI_tol = 0.1, convex_def=0.01, convex_def_tol=0.001, max_footprint=None, proj_kwargs=None, - periodic=(False, False), + periodic=(False, False), #False max_iters=1000, min_limit_diff=1e-10): """Find the largest convex contour around a maximum. @@ -287,15 +340,29 @@ def convex_contour_around_maximum(data, ji, init_contour_step_frac=0.1, The 2D data to contour ji : tuple The index of the maximum in (j, i) order + lx : array_like + Array with shape (N,n,n) including the x position of particles at N time steps(position array is n*n, + and trajectories were unwrapped to correct the jump of displacement at the periodic boundaries) + ly : array_like + Array with shape (N,n,n) including the y position of particles at N time steps(position array is n*n, + and trajectories were unwrapped to correct the jump of displacement at the periodic boundaries) + dx, dy : float + Particle spacings in x and y directions, must have the same units as lx and ly, respectively. init_contour_step_frac : float the value with which to increment the initial contour level (multiplied by the local maximum value) border: int the initial window around the maximum + max_width: int + The maximum width of the search window. Value depends on the resolution of particles convex_def : float, optional The target convexity deficiency allowed for the contour. convex_def_tol : float, optional The tolerance for which the convexity deficiency will be sought + CI_th : float, optional + The target coherency index allowed for the contour. Set it as -np.inf if don't need it. + CI_tol : float, optional + The tolerance for which the coherency index will be sought verbose : bool, optional Whether to print out diagnostic information proj_kwargs : dict, optional @@ -331,20 +398,31 @@ def convex_contour_around_maximum(data, ji, init_contour_step_frac=0.1, logger.debug("init_contour_step_frac %g" % init_contour_step_frac) logger.debug("init_contour_step_size %g" % init_contour_step_size) - lower_lim = 0 + lower_lim = 0 #max_value upper_lim = 2*init_contour_step_size # special logic needed to find the first contour greater than convex_def - exceeded = False + exceeded = False #False + overlap = False + CI = -np.inf cd = np.inf contour = None region_area = None for n_iter in range(max_iters): - logger.debug('iter %g, cd %g' % (n_iter, cd)) - if abs(cd - convex_def) <= convex_def_tol: + logger.debug('iter %g, CI %g, cd %g' % (n_iter, CI,cd)) +# logger.debug('iter %g, cd %g' % (n_iter, cd)) + if (abs(CI - CI_th) <= CI_tol) and (abs(cd - convex_def) <= convex_def_tol): + + logger.debug('CI %g is close to target %g within tolerance %g' % + (CI, CI_th, CI_tol)) logger.debug('cd %g is close to target %g within tolerance %g' % (cd, convex_def, convex_def_tol)) break + +# elif abs(cd - convex_def) <= convex_def_tol: +# logger.debug('cd %g is close to target %g within tolerance %g' % +# (cd, convex_def, convex_def_tol)) +# break # the new contour level to try logger.debug('current lims: (%10.20f, %10.20f)' % (lower_lim, upper_lim)) @@ -352,14 +430,19 @@ def convex_contour_around_maximum(data, ji, init_contour_step_frac=0.1, if (upper_lim - lower_lim) < min_limit_diff: logger.debug('limit diff below threshold; ending search') - break + level = lower_lim # switch to the lower_lim (cd or CI of lower_lim and upper_lim might have big difference + # although the contour values are very close) + overlap = True + if level==0: #this might be a weird situation + break + logger.debug(('contouring level: %10.20f border: ' % level) + repr(border_j) + repr(border_i)) try: # try to get a contour contour, region_data, border_j, border_i = find_contour_around_maximum( - data, (j, i), level, border_j, border_i, + data, (j, i), level, border_j, border_i, max_width=max_width, max_footprint=max_footprint, periodic=periodic) except ValueError as ve: # we will probably be here if the contour search ended up covering @@ -369,6 +452,9 @@ def convex_contour_around_maximum(data, ji, init_contour_step_frac=0.1, #if contour is None: upper_lim = level exceeded = True + + if overlap: # limit diff below threshold; ending search here + break continue #else: # break @@ -380,12 +466,30 @@ def convex_contour_around_maximum(data, ji, init_contour_step_frac=0.1, contour_proj = project_vertices(contour, **proj_kwargs) region_area, hull_area, cd = contour_area(contour_proj) - logger.debug('region_area: % 6.1f, hull_area: % 6.1f, convex_def: % 6.5e' - % (region_area, hull_area, cd)) + # get the coherency index + CI = contour_CI(lx, ly,contour_proj, ji, region_area, dx, dy, border_i, border_j) + + logger.debug('region_area: % 6.1f, hull_area: % 6.1f, convex_def: % 6.5e, CI: % 6.5e' + % (region_area, hull_area, cd, CI)) + + if overlap: + if (CI > CI_th and cd < convex_def): + break # limit diff below threshold; ending search here + else: + border_j = (border, border) + border_i = (border, border) + continue + + +# if exceeded and (region_area < 2): # contous area is too small; ending search here +# logger.debug('contour area is too small') +# break + # special logic needed to find the first contour greater than convex_def if not exceeded: - if cd < convex_def: + if (CI > CI_th and cd < convex_def): +# if cd < convex_def: # need to keep upper_lim until we exceed the convex def lower_lim = level upper_lim = level + 2*init_contour_step_size @@ -396,7 +500,8 @@ def convex_contour_around_maximum(data, ji, init_contour_step_frac=0.1, else: # from here we can assume that the target contour lies between # lower_lim and_upper_lim - if cd < convex_def: + if (CI > CI_th and cd < convex_def): +# if cd < convex_def: lower_lim = level else: upper_lim = level @@ -405,19 +510,29 @@ def convex_contour_around_maximum(data, ji, init_contour_step_frac=0.1, if contour is not None: contour[:, 0] += (j-border_j[0]) contour[:, 1] += (i-border_i[0]) - return contour, region_area, cd + return contour, region_area, cd, CI -def find_convex_contours(data, min_distance=5, min_area=100., - use_threadpool=False, lon=None, lat=None, +def find_convex_contours(data, lx, ly, dx, dy, min_distance=5, min_area=100., CI_th = -1.0, CI_tol = 0.1, + convex_def=0.01, convex_def_tol=0.001, + init_contour_step_frac=0.1, min_limit_diff=1e-10, max_width=100, + use_threadpool=False, lon=None, lat=None, #use_multi_process=False, progress=False, **contour_kwargs): """Find the outermost convex contours around the maxima of - data with specified convexity deficiency. + data with specified convexity deficiency. Parameters ---------- data : array_like The 2D data to contour + lx : array_like + Array with shape (N,n,n) including the x position of particles at N time steps(position array is n*n, + and trajectories were unwrapped to correct the jump of displacement at the periodic boundaries) + ly : array_like + Array with shape (N,n,n) including the y position of particles at N time steps(position array is n*n, + and trajectories were unwrapped to correct the jump of displacement at the periodic boundaries) + dx, dy : float + Particle spacings in x and y directions, must have the same units as lx and ly, respectively. min_distance : int, optional The minimum distance around maxima (pixel units) min_area : float, optional @@ -431,10 +546,16 @@ def find_convex_contours(data, min_distance=5, min_area=100., (multiplied by the local maximum value) border: int the initial window around the maximum + max_width: int + The maximum width of the search window. Value depends on the resolution of particles convex_def : float, optional The target convexity deficiency allowed for the contour. convex_def_tol : float, optional The tolerance for which the convexity deficiency will be sought + CI_th : float, optional + The target coherency index allowed for the contour. Set it as -np.inf if don't need it. + CI_tol : float, optional + The tolerance for which the coherency index will be sought verbose : bool, optional Whether to print out diagnostic information proj_kwargs : dict, optional @@ -474,6 +595,10 @@ def find_convex_contours(data, min_distance=5, min_area=100., from multiprocessing.pool import ThreadPool pool = ThreadPool() map_function = pool.imap_unordered +# if use_multi_process: +# from pathos.multiprocessing import ProcessingPool as Pool +# #map_function = Pool(processes = 4).imap_unordered +# map_function = Pool().uimap else: try: from itertools import imap @@ -496,10 +621,13 @@ def maybe_contour_maximum(ji): if 'proj_kwargs' in contour_kwargs: del contour_kwargs['proj_kwargs'] - contour, area, cd = convex_contour_around_maximum(data, ji, - **contour_kwargs) + contour, area, cd, CI = convex_contour_around_maximum(data, lx, ly, ji, dx=dx, dy=dy, + max_width=max_width,init_contour_step_frac=init_contour_step_frac, + min_limit_diff=min_limit_diff,CI_th = CI_th, CI_tol = CI_tol, + convex_def=convex_def, convex_def_tol=convex_def_tol, + **contour_kwargs) if area and (area >= min_area): - result = ji, contour, area, cd + result = ji, contour, area, cd, CI toc = time() logger.debug("point " + repr(tuple(ji)) + " took %g s" % (toc-tic)) return result