diff --git a/doc/source/pysteps_reference/utils.rst b/doc/source/pysteps_reference/utils.rst index f43454bfe..cd13c95ca 100644 --- a/doc/source/pysteps_reference/utils.rst +++ b/doc/source/pysteps_reference/utils.rst @@ -7,8 +7,11 @@ Implementation of miscellaneous utility functions. .. automodule:: pysteps.utils.interface .. automodule:: pysteps.utils.arrays +.. automodule:: pysteps.utils.cleansing .. automodule:: pysteps.utils.conversion .. automodule:: pysteps.utils.dimension .. automodule:: pysteps.utils.fft +.. automodule:: pysteps.utils.images +.. automodule:: pysteps.utils.interpolate .. automodule:: pysteps.utils.spectral .. automodule:: pysteps.utils.transformation diff --git a/examples/LK_buffer_mask.py b/examples/LK_buffer_mask.py index a6b38b35d..2a3c90a72 100644 --- a/examples/LK_buffer_mask.py +++ b/examples/LK_buffer_mask.py @@ -1,15 +1,16 @@ +# -*- coding: utf-8 -*- """ Handling of no-data in Lucas-Kanade =================================== Areas of missing data in radar images are typically caused by visibility limits -such as beam blockage and the radar coverage itself. These artifacts can mislead +such as beam blockage and the radar coverage itself. These artifacts can mislead the echo tracking algorithms. For instance, precipitation leaving the domain might be erroneously detected as having nearly stationary velocity. -This example shows how the Lucas-Kanade algorithm can be tuned to avoid the -erroneous interpretation of velocities near the maximum range of the radars by -buffering the no-data mask in the radar image in order to exclude all vectors +This example shows how the Lucas-Kanade algorithm can be tuned to avoid the +erroneous interpretation of velocities near the maximum range of the radars by +buffering the no-data mask in the radar image in order to exclude all vectors detected nearby no-data areas. """ @@ -71,7 +72,9 @@ mask[~np.isnan(ref_mm)] = np.nan # Log-transform the data [dBR] -R, metadata = transformation.dB_transform(R, metadata, threshold=0.1, zerovalue=-15.0) +R, metadata = transformation.dB_transform( + R, metadata, threshold=0.1, zerovalue=-15.0 +) # Keep the reference frame in dBR (for plotting purposes) ref_dbr = R[0].copy() @@ -109,10 +112,20 @@ R.data[R.mask] = np.nan # Use default settings (i.e., no buffering of the radar mask) -x, y, u, v = LK_optflow(R, dense=False, buffer_mask=0, quality_level_ST=0.1) +fd_kwargs1 = {"buffer_mask":0} +xy, uv = LK_optflow(R, dense=False, fd_kwargs=fd_kwargs1) plt.imshow(ref_dbr, cmap=plt.get_cmap("Greys")) plt.imshow(mask, cmap=colors.ListedColormap(["black"]), alpha=0.5) -plt.quiver(x, y, u, v, color="red", angles="xy", scale_units="xy", scale=0.2) +plt.quiver( + xy[:, 0], + xy[:, 1], + uv[:, 0], + uv[:, 1], + color="red", + angles="xy", + scale_units="xy", + scale=0.2, +) circle = plt.Circle((620, 245), 100, color="b", clip_on=False, fill=False) plt.gca().add_artist(circle) plt.title("buffer_mask = 0 (default)") @@ -129,13 +142,24 @@ # 'x,y,u,v = LK_optflow(.....)'. # with buffer -x, y, u, v = LK_optflow(R, dense=False, buffer_mask=20, quality_level_ST=0.2) +buffer = 10 +fd_kwargs2 = {"buffer_mask":buffer} +xy, uv = LK_optflow(R, dense=False, fd_kwargs=fd_kwargs2) plt.imshow(ref_dbr, cmap=plt.get_cmap("Greys")) plt.imshow(mask, cmap=colors.ListedColormap(["black"]), alpha=0.5) -plt.quiver(x, y, u, v, color="red", angles="xy", scale_units="xy", scale=0.2) +plt.quiver( + xy[:, 0], + xy[:, 1], + uv[:, 0], + uv[:, 1], + color="red", + angles="xy", + scale_units="xy", + scale=0.2, +) circle = plt.Circle((620, 245), 100, color="b", clip_on=False, fill=False) plt.gca().add_artist(circle) -plt.title("buffer_mask = 20") +plt.title("buffer_mask = %i" % buffer) plt.show() ################################################################################ @@ -148,13 +172,13 @@ # the negative bias that is introduced by the the erroneous interpretation of # velocities near the maximum range of the radars. -UV1 = LK_optflow(R, dense=True, buffer_mask=0, quality_level_ST=0.1) -UV2 = LK_optflow(R, dense=True, buffer_mask=20, quality_level_ST=0.2) +UV1 = LK_optflow(R, dense=True, fd_kwargs=fd_kwargs1) +UV2 = LK_optflow(R, dense=True, fd_kwargs=fd_kwargs2) V1 = np.sqrt(UV1[0] ** 2 + UV1[1] ** 2) V2 = np.sqrt(UV2[0] ** 2 + UV2[1] ** 2) -plt.imshow((V1 - V2) / V2, cmap=cm.RdBu_r, vmin=-0.1, vmax=0.1) +plt.imshow((V1 - V2) / V2, cmap=cm.RdBu_r, vmin=-0.5, vmax=0.5) plt.colorbar(fraction=0.04, pad=0.04) plt.title("Relative difference in motion speed") plt.show() @@ -184,7 +208,13 @@ # Find the veriyfing observations in the archive fns = io.archive.find_by_date( - date, root_path, path_fmt, fn_pattern, fn_ext, timestep=5, num_next_files=12 + date, + root_path, + path_fmt, + fn_pattern, + fn_ext, + timestep=5, + num_next_files=12, ) # Read and convert the radar composites @@ -200,8 +230,8 @@ score_2.append(skill(R_f2[i, :, :], R_o[i + 1, :, :])["corr_s"]) x = (np.arange(12) + 1) * 5 # [min] -plt.plot(x, score_1, label="no mask buffer") -plt.plot(x, score_2, label="with mask buffer") +plt.plot(x, score_1, label="buffer_mask = 0") +plt.plot(x, score_2, label="buffer_mask = %i" % buffer) plt.legend() plt.xlabel("Lead time [min]") plt.ylabel("Corr. coeff. []") diff --git a/pysteps/decorators.py b/pysteps/decorators.py new file mode 100644 index 000000000..166eaca45 --- /dev/null +++ b/pysteps/decorators.py @@ -0,0 +1,56 @@ +""" +pysteps.decorators +================== + +Decorators used to define reusable building blocks that can change or extend +the behavior of some functions in pysteps. + +.. autosummary:: + :toctree: ../generated/ + + check_motion_input_image +""" +from functools import wraps + +import numpy as np + + +def check_input_frames(minimum_input_frames=2, + maximum_input_frames=np.inf, + just_ndim=False): + """ + Check that the input_images used as inputs in the optical-flow + methods has the correct shape (t, x, y ). + """ + + def _check_input_frames(motion_method_func): + @wraps(motion_method_func) + def new_function(*args, **kwargs): + """ + Return new function with the checks prepended to the + target motion_method_func function. + """ + + input_images = args[0] + if input_images.ndim != 3: + raise ValueError( + "input_images dimension mismatch.\n" + f"input_images.shape: {str(input_images.shape)}\n" + "(t, x, y ) dimensions expected" + ) + + if not just_ndim: + num_of_frames = input_images.shape[0] + + if minimum_input_frames < num_of_frames > maximum_input_frames: + raise ValueError( + f"input_images frames {num_of_frames} mismatch.\n" + f"Minimum frames: {minimum_input_frames}\n" + f"Maximum frames: {maximum_input_frames}\n" + ) + + return motion_method_func(*args, **kwargs) + + return new_function + + return _check_input_frames diff --git a/pysteps/motion/darts.py b/pysteps/motion/darts.py index 085cef25c..1391a739e 100644 --- a/pysteps/motion/darts.py +++ b/pysteps/motion/darts.py @@ -11,18 +11,23 @@ """ import sys -import time + import numpy as np +import time from numpy.linalg import lstsq, svd -from .. import utils -def DARTS(R, **kwargs): +from pysteps import utils +from pysteps.decorators import check_input_frames + + +@check_input_frames(just_ndim=True) +def DARTS(input_images, **kwargs): """Compute the advection field from a sequence of input images by using the DARTS method. :cite:`RCW2011` Parameters ---------- - R : array-like + input_images : array-like Array of shape (T,m,n) containing a sequence of T two-dimensional input images of shape (m,n). @@ -51,7 +56,7 @@ def DARTS(R, **kwargs): n_threads : int Number of threads to use for the FFT computation. Applicable if fft_method is 'pyfftw'. - print_info : bool + verbose : bool If True, print information messages. lsq_method : {1, 2} The method to use for solving the linear equations in the least squares @@ -76,12 +81,11 @@ def DARTS(R, **kwargs): M_y = kwargs.get("M_y", 2) fft_method = kwargs.get("fft_method", "numpy") output_type = kwargs.get("output_type", "spatial") - print_info = kwargs.get("print_info", False) lsq_method = kwargs.get("lsq_method", 2) verbose = kwargs.get("verbose", True) - if N_t >= R.shape[0]: - raise ValueError("N_t = %d >= %d = T, but N_t < T required" % (N_t, R.shape[0])) + if N_t >= input_images.shape[0]: + raise ValueError("N_t = %d >= %d = T, but N_t < T required" % (N_t, input_images.shape[0])) if output_type not in ["spatial", "spectral"]: raise ValueError("invalid output_type=%s, must be 'spatial' or 'spectral'" % output_type) @@ -90,16 +94,16 @@ def DARTS(R, **kwargs): print("Computing the motion field with the DARTS method.") t0 = time.time() - R = np.moveaxis(R, (0, 1, 2), (2, 0, 1)) + input_images = np.moveaxis(input_images, (0, 1, 2), (2, 0, 1)) - fft = utils.get_method(fft_method, shape=R.shape[:2], fftn_shape=R.shape, + fft = utils.get_method(fft_method, shape=input_images.shape[:2], fftn_shape=input_images.shape, **kwargs) - T_x = R.shape[1] - T_y = R.shape[0] - T_t = R.shape[2] + T_x = input_images.shape[1] + T_y = input_images.shape[0] + T_t = input_images.shape[2] - if print_info: + if verbose: print("-----") print("DARTS") print("-----") @@ -108,45 +112,43 @@ def DARTS(R, **kwargs): sys.stdout.flush() starttime = time.time() - R = fft.fftn(R) + input_images = fft.fftn(input_images) - if print_info: + if verbose: print("Done in %.2f seconds." % (time.time() - starttime)) print(" Constructing the y-vector..."), sys.stdout.flush() starttime = time.time() - m = (2*N_x+1)*(2*N_y+1)*(2*N_t+1) - n = (2*M_x+1)*(2*M_y+1) + m = (2 * N_x + 1) * (2 * N_y + 1) * (2 * N_t + 1) + n = (2 * M_x + 1) * (2 * M_y + 1) y = np.zeros(m, dtype=complex) - k_t, k_y, k_x = np.unravel_index(np.arange(m), (2*N_t+1, 2*N_y+1, 2*N_x+1)) + k_t, k_y, k_x = np.unravel_index(np.arange(m), (2 * N_t + 1, 2 * N_y + 1, 2 * N_x + 1)) for i in range(m): k_x_ = k_x[i] - N_x k_y_ = k_y[i] - N_y k_t_ = k_t[i] - N_t - R_ = R[k_y_, k_x_, k_t_] - - y[i] = k_t_ * R_ + y[i] = k_t_ * input_images[k_y_, k_x_, k_t_] - if print_info: + if verbose: print("Done in %.2f seconds." % (time.time() - starttime)) A = np.zeros((m, n), dtype=complex) B = np.zeros((m, n), dtype=complex) - if print_info: + if verbose: print(" Constructing the H-matrix..."), sys.stdout.flush() starttime = time.time() - c1 = -1.0*T_t / (T_x * T_y) + c1 = -1.0 * T_t / (T_x * T_y) - kp_y, kp_x = np.unravel_index(np.arange(n), (2*M_y+1, 2*M_x+1)) + kp_y, kp_x = np.unravel_index(np.arange(n), (2 * M_y + 1, 2 * M_x + 1)) for i in range(m): k_x_ = k_x[i] - N_x @@ -159,7 +161,7 @@ def DARTS(R, **kwargs): i_ = k_y_ - kp_y_ j_ = k_x_ - kp_x_ - R_ = R[i_, j_, k_t_] + R_ = input_images[i_, j_, k_t_] c2 = c1 / T_y * i_ A[i, :] = c2 * R_ @@ -167,7 +169,7 @@ def DARTS(R, **kwargs): c2 = c1 / T_x * j_ B[i, :] = c2 * R_ - if print_info: + if verbose: print("Done in %.2f seconds." % (time.time() - starttime)) print(" Solving the linear systems..."), @@ -179,40 +181,39 @@ def DARTS(R, **kwargs): else: x = _leastsq(A, B, y) - if print_info: + if verbose: print("Done in %.2f seconds." % (time.time() - starttime)) - h, w = 2*M_y+1, 2*M_x+1 + h, w = 2 * M_y + 1, 2 * M_x + 1 U = np.zeros((h, w), dtype=complex) V = np.zeros((h, w), dtype=complex) - i, j = np.unravel_index(np.arange(h*w), (h, w)) + i, j = np.unravel_index(np.arange(h * w), (h, w)) - V[i, j] = x[0:h*w] - U[i, j] = x[h*w:2*h*w] + V[i, j] = x[0:h * w] + U[i, j] = x[h * w:2 * h * w] - k_x, k_y = np.meshgrid(np.arange(-M_x, M_x+1), np.arange(-M_y, M_y+1)) + k_x, k_y = np.meshgrid(np.arange(-M_x, M_x + 1), np.arange(-M_y, M_y + 1)) if output_type == "spatial": - U = np.real(fft.ifft2(_fill(U, R.shape[0], R.shape[1], k_x, k_y))) - V = np.real(fft.ifft2(_fill(V, R.shape[0], R.shape[1], k_x, k_y))) + U = np.real(fft.ifft2(_fill(U, input_images.shape[0], input_images.shape[1], k_x, k_y))) + V = np.real(fft.ifft2(_fill(V, input_images.shape[0], input_images.shape[1], k_x, k_y))) if verbose: print("--- %s seconds ---" % (time.time() - t0)) return np.stack([U, V]) + def _leastsq(A, B, y): M = np.hstack([A, B]) M_ct = M.conjugate().T MM = np.dot(M_ct, M) - M = None - U, s, V = svd(MM, full_matrices=False) - MM = None - mask = s > 0.01*s[0] + + mask = s > 0.01 * s[0] s = 1.0 / s[mask] MM_inv = np.dot(np.dot(V[:len(s), :].conjugate().T, np.diag(s)), @@ -220,6 +221,7 @@ def _leastsq(A, B, y): return np.dot(MM_inv, np.dot(M_ct, y)) + def _fill(X, h, w, k_x, k_y): X_f = np.zeros((h, w), dtype=complex) X_f[k_y, k_x] = X diff --git a/pysteps/motion/lucaskanade.py b/pysteps/motion/lucaskanade.py index fee0b89a7..db0eccf83 100644 --- a/pysteps/motion/lucaskanade.py +++ b/pysteps/motion/lucaskanade.py @@ -1,887 +1,440 @@ +# -*- coding: utf-8 -*- """ pysteps.motion.lucaskanade ========================== -OpenCV implementation of the Lucas-Kanade method with interpolated motion -vectors for areas with no precipitation. +The Lucas-Kanade (LK) local feature tracking module. + +This module implements the interface to the local `Lucas-Kanade`_ routine +available in OpenCV_. + +For its dense method, it additionally interpolates the sparse vectors over a +regular grid to return a motion field. + +.. _OpenCV: https://opencv.org/ + +.. _`Lucas-Kanade`:\ + https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga473e4b886d0bcc6b65831eb88ed93323 .. autosummary:: :toctree: ../generated/ dense_lucaskanade - + track_features """ import numpy as np from numpy.ma.core import MaskedArray + +from pysteps.decorators import check_input_frames from pysteps.exceptions import MissingOptionalDependency +from pysteps import utils +from pysteps.utils.cleansing import decluster, detect_outliers +from pysteps.utils.images import morph_opening + try: import cv2 - cv2_imported = True + CV2_IMPORTED = True except ImportError: - cv2_imported = False -import scipy.spatial + CV2_IMPORTED = False + import time import warnings -def dense_lucaskanade(R, **kwargs): - """ - .. _opencv: https://opencv.org/ +@check_input_frames(2) +def dense_lucaskanade(input_images, + lk_kwargs=None, + fd_method="ShiTomasi", + fd_kwargs=None, + interp_method="rbfinterp2d", + interp_kwargs=None, + dense=True, + nr_std_outlier=3, + k_outlier=30, + size_opening=3, + decl_scale=10, + verbose=False): + """Run the Lucas-Kanade optical flow routine and interpolate the motion + vectors. - .. _`Lucas-Kanade`: https://docs.opencv.org/3.4/dc/d6b/\ - group__video__track.html#ga473e4b886d0bcc6b65831eb88ed93323 + .. _OpenCV: https://opencv.org/ - OpenCV_ implementation of the local `Lucas-Kanade`_ method with - interpolation of the sparse motion vectors to fill the whole grid. + .. _`Lucas-Kanade`:\ + https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga473e4b886d0bcc6b65831eb88ed93323 - .. _MaskedArray: https://docs.scipy.org/doc/numpy/reference/\ - maskedarray.baseclass.html#numpy.ma.MaskedArray + .. _MaskedArray:\ + https://docs.scipy.org/doc/numpy/reference/maskedarray.baseclass.html#numpy.ma.MaskedArray + + Interface to the OpenCV_ implementation of the local `Lucas-Kanade`_ optical + flow method applied in combination to a feature detection routine. + + The sparse motion vectors are finally interpolated to return the whole + motion field. - .. _ndarray:\ - https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html - - .. _Shi-Tomasi: https://docs.opencv.org/3.4.1/dd/d1a/group__\ - imgproc__feature.html#ga1d6bb77486c8f92d79c8793ad995d541 - Parameters ---------- - R : ndarray_ or MaskedArray_ - Array of shape (T,m,n) containing a sequence of T two-dimensional input - images of shape (m,n). - In case of an ndarray_, invalid values (Nans or infs) are masked. - The mask in the MaskedArray_ defines a region where velocity vectors are - not computed. + input_images : array_like or MaskedArray_ + Array of shape (T, m, n) containing a sequence of *T* two-dimensional + input images of shape (m, n). The indexing order in **input_images** is + assumed to be (time, latitude, longitude). + + *T* = 2 is the minimum required number of images. + With *T* > 2, all the resulting sparse vectors are pooled together for + the final interpolation on a regular grid. + + In case of array_like, invalid values (Nans or infs) are masked, + otherwise the mask of the MaskedArray_ is used. Such mask defines a + region where features are not detected for the tracking algorithm. + + lk_kwargs : dict, optional + Optional dictionary containing keyword arguments for the `Lucas-Kanade`_ + features tracking algorithm. See the documentation of + :py:func:`pysteps.motion.lucaskanade.track_features`. + + fd_method : {"ShiTomasi"}, optional + Name of the feature detection routine. See feature detection methods in + :py:mod:`pysteps.utils.images`. + + fd_kwargs : dict, optional + Optional dictionary containing keyword arguments for the features + detection algorithm. + See the documentation of :py:mod:`pysteps.utils.images`. + + interp_method : {"rbfinterp2d"}, optional + Name of the interpolation method to use. See interpolation methods in + :py:mod:`pysteps.utils.interpolate`. + + interp_kwargs : dict, optional + Optional dictionary containing keyword arguments for the interpolation + algorithm. See the documentation of :py:mod:`pysteps.utils.interpolate`. - Other Parameters - ---------------- dense : bool, optional - If True (the default), it returns the three-dimensional array (2,m,n) - containing the dense x- and y-components of the motion field. If false, - it returns the sparse motion vectors as 1D arrays x, y, u, v, where - x, y define the vector locations, u, v define the x and y direction - components of the vectors. - - buffer_mask : int, optional - A mask buffer width in pixels. This extends the input mask (if any) - to help avoiding the erroneous interpretation of velocities near the - maximum range of the radars (0 by default). - - max_corners_ST : int, optional - The maxCorners parameter in the `Shi-Tomasi`_ corner detection method. - It represents the maximum number of points to be tracked (corners), - by default this is 500. If set to zero, all detected corners are used. - - quality_level_ST : float, optional - The qualityLevel parameter in the `Shi-Tomasi`_ corner detection method. - It represents the minimal accepted quality for the points to be tracked - (corners), by default this is set to 0.1. Higher quality thresholds can - lead to no detection at all. - - min_distance_ST : int, optional - The minDistance parameter in the `Shi-Tomasi`_ corner detection method. - It represents minimum possible Euclidean distance in pixels - between corners, by default this is set to 3 pixels. - - block_size_ST : int, optional - The blockSize parameter in the `Shi-Tomasi`_ corner detection method. - It represents the window size in pixels used for computing a derivative - covariation matrix over each pixel neighborhood, by default this is set - to 15 pixels. - - winsize_LK : tuple of int, optional - The winSize parameter in the `Lucas-Kanade`_ optical flow method. - It represents the size of the search window that it is used at each - pyramid level, by default this is set to (50, 50) pixels. + If True, return the three-dimensional array (2, m, n) containing + the dense x- and y-components of the motion field. - nr_levels_LK : int, optional - The maxLevel parameter in the `Lucas-Kanade`_ optical flow method. - It represents the 0-based maximal pyramid level number, by default this - is set to 3. + If False, return the sparse motion vectors as 2-D **xy** and **uv** + arrays, where **xy** defines the vector positions, **uv** defines the + x and y direction components of the vectors. nr_std_outlier : int, optional - Maximum acceptable deviation from the mean/median in terms of - number of standard deviations. Any anomaly larger than - this value is flagged as outlier and excluded from the interpolation. - By default this is set to 3. - - multivariate_outlier : bool, optional - If true (the default), the outlier detection is computed in terms of - the Mahalanobis distance. If false, the outlier detection is simply - computed in terms of velocity. - - k_outlier : int, optional + Maximum acceptable deviation from the mean in terms of number of + standard deviations. Any sparse vector with a deviation larger than + this threshold is flagged as outlier and excluded from the + interpolation. + See the documentation of + :py:func:`pysteps.utils.cleansing.detect_outliers`. + + k_outlier : int or None, optional The number of nearest neighbours used to localize the outlier detection. - If set equal to 0, it employs all the data points. - The default is 30. + If set to None, it employs all the data points (global detection). + See the documentation of + :py:func:`pysteps.utils.cleansing.detect_outliers`. size_opening : int, optional The size of the structuring element kernel in pixels. This is used to perform a binary morphological opening on the input fields in order to - filter isolated echoes due to clutter. By default this is set to 3. - If set to zero, the fitlering is not perfomed. - - decl_grid : int, optional - The cell size in pixels of the declustering grid that is used to filter - out outliers in a sparse motion field and get more representative data - points before the interpolation. This simply computes new sparse vectors - over a coarser grid by taking the median of all vectors within one cell. - By default this is set to 20 pixels. If set to less than 2 pixels, the - declustering is not perfomed. - - min_nr_samples : int, optional - The minimum number of samples necessary for computing the median vector - within given declustering cell, otherwise all sparse vectors in that - cell are discarded. By default this is set to 2. - - rbfunction : string, optional - The name of the radial basis function used for the interpolation of the - sparse vectors. This is based on the Euclidian norm d. By default this - is set to "inverse" and the available names are "nearest", "inverse", - "gaussian". - - k : int, optional - The number of nearest neighbours used for fast interpolation, by default - this is set to 20. If set equal to zero, it employs all the neighbours. - - epsilon : float, optional - The adjustable constant used in the gaussian and inverse radial basis - functions. by default this is computed as the median distance between - the sparse vectors. - - nchunks : int, optional - Split the grid points in n chunks to limit the memory usage during the - interpolation. By default this is set to 5, if set to 1 the interpolation - is computed with the whole grid. - - extra_vectors : ndarray_, optional - Additional sparse motion vectors as 2d array (columns: x,y,u,v; rows: - nbr. of vectors) to be integrated with the sparse vectors from the - Lucas-Kanade local tracking. - x and y must be in pixel coordinates, with (0,0) being the upper-left - corner of the field R. u and v must be in pixel units. By default this - is set to None. + filter isolated echoes due to clutter. If set to zero, the filtering + is not perfomed. + See the documentation of + :py:func:`pysteps.utils.images.morph_opening`. + + decl_scale : int, optional + The scale declustering parameter in pixels used to reduce the number of + redundant sparse vectors before the interpolation. + Sparse vectors within this declustering scale are averaged together. + If set to less than 2 pixels, the declustering is not perfomed. + See the documentation of + :py:func:`pysteps.cleansing.cleansing.decluster`. verbose : bool, optional - If set to True, it prints information about the program (True by default). + If set to True, print some information about the program. Returns ------- - out : ndarray_ - If dense=True (the default), it returns the three-dimensional array (2,m,n) - containing the dense x- and y-components of the motion field in units of - pixels / timestep as given by the input array R. - If dense=False, it returns a tuple containing the one-dimensional arrays - x, y, u, v, where x, y define the vector locations, u, v define the x - and y direction components of the vectors. - Return an empty array when no motion vectors are found. + + out : array_like or tuple + If **dense=True** (the default), it returns the three-dimensional array + (2, m, n) containing the dense x- and y-components of the motion field + in units of [pixels/timestep] as given by the input array input_images. + Return a zero motion field of shape (2, m, n) when no motion is + detected. + + If **dense=False**, it returns a tuple containing the 2-dimensional + arrays **xy** and **uv**, where x, y define the vector locations, + u, v define the x and y direction components of the vectors. + Return two empty arrays when no motion is detected. + + See also + -------- + + pysteps.motion.lucaskanade.track_features References ---------- - - Bouguet, J.-Y.: Pyramidal implementation of the affine Lucas Kanade - feature tracker description of the algorithm, Intel Corp., 5, 4, + + Bouguet, J.-Y.: Pyramidal implementation of the affine Lucas Kanade + feature tracker description of the algorithm, Intel Corp., 5, 4, https://doi.org/10.1109/HPDC.2004.1323531, 2001 - Lucas, B. D. and Kanade, T.: An iterative image registration technique with - an application to stereo vision, in: Proceedings of the 1981 DARPA Imaging + Lucas, B. D. and Kanade, T.: An iterative image registration technique with + an application to stereo vision, in: Proceedings of the 1981 DARPA Imaging Understanding Workshop, pp. 121–130, 1981. - """ - if len(R.shape) != 3: - raise ValueError( - "R has %i dimensions, but a three-dimensional array is expected" - % len(R.shape) - ) - - if R.shape[0] < 2: - raise ValueError( - "R has %i frame, but at least two frames are expected" % R.shape[0] - ) - - R = R.copy() - - # defaults - dense = kwargs.get("dense", True) - max_corners_ST = kwargs.get("max_corners_ST", 500) - quality_level_ST = kwargs.get("quality_level_ST", 0.1) - min_distance_ST = kwargs.get("min_distance_ST", 3) - block_size_ST = kwargs.get("block_size_ST", 15) - winsize_LK = kwargs.get("winsize_LK", (50, 50)) - nr_levels_LK = kwargs.get("nr_levels_LK", 3) - nr_std_outlier = kwargs.get("nr_std_outlier", 3) - nr_IQR_outlier = kwargs.get("nr_IQR_outlier", None) - if nr_IQR_outlier is not None: - nr_std_outlier = nr_IQR_outlier - warnings.warn( - "the 'nr_IQR_outlier' argument will be deprecated in the next release; " - + "use 'nr_std_outlier' instead.", - category=FutureWarning, - ) - multivariate_outlier = kwargs.get("multivariate_outlier", True) - k_outlier = kwargs.get("k_outlier", 30) - size_opening = kwargs.get("size_opening", 3) - decl_grid = kwargs.get("decl_grid", 20) - min_nr_samples = kwargs.get("min_nr_samples", 2) - rbfunction = kwargs.get("rbfunction", "inverse") - k = kwargs.get("k", 50) - epsilon = kwargs.get("epsilon", None) - nchunks = kwargs.get("nchunks", 5) - extra_vectors = kwargs.get("extra_vectors", None) - if extra_vectors is not None: - if len(extra_vectors.shape) != 2: - raise ValueError( - "extra_vectors has %i dimensions, but 2 dimensions are expected" - % len(extra_vectors.shape) - ) - if extra_vectors.shape[1] != 4: - raise ValueError( - "extra_vectors has %i columns, but 4 columns are expected" - % extra_vectors.shape[1] - ) - verbose = kwargs.get("verbose", True) - buffer_mask = kwargs.get("buffer_mask", 0) + input_images = input_images.copy() if verbose: print("Computing the motion field with the Lucas-Kanade method.") t0 = time.time() - # Get mask - if isinstance(R, MaskedArray): - mask = np.ma.getmaskarray(R).copy() - else: - R = np.ma.masked_invalid(R) - mask = np.ma.getmaskarray(R).copy() - R[mask] = np.nanmin(R) # Remove any Nan from the raw data - - nr_fields = R.shape[0] - domain_size = (R.shape[1], R.shape[2]) - y0Stack = [] - x0Stack = [] - uStack = [] - vStack = [] - for n in range(nr_fields - 1): + nr_fields = input_images.shape[0] + domain_size = (input_images.shape[1], input_images.shape[2]) - # extract consecutive images - prvs = R[n, :, :].copy() - next = R[n + 1, :, :].copy() + feature_detection_method = utils.get_method(fd_method) + interpolation_method = utils.get_method(interp_method) - # skip loop if no precip - if ~np.any(prvs > prvs.min()) or ~np.any(next > next.min()): - continue + if fd_kwargs is None: + fd_kwargs = dict() - # scale between 0 and 255 - prvs = (prvs - prvs.min()) / (prvs.max() - prvs.min()) * 255 - next = (next - next.min()) / (next.max() - next.min()) * 255 + if lk_kwargs is None: + lk_kwargs = dict() + + if interp_kwargs is None: + interp_kwargs = dict() + + xy = np.empty(shape=(0, 2)) + uv = np.empty(shape=(0, 2)) + for n in range(nr_fields - 1): + + # extract consecutive images + prvs_img = input_images[n, :, :].copy() + next_img = input_images[n + 1, :, :].copy() - # convert to 8-bit - prvs = np.ndarray.astype(prvs, "uint8") - next = np.ndarray.astype(next, "uint8") - mask_ = np.ndarray.astype(mask[n, :, :], "uint8") + if ~isinstance(prvs_img, MaskedArray): + prvs_img = np.ma.masked_invalid(prvs_img) + np.ma.set_fill_value(prvs_img, prvs_img.min()) - # buffer the quality mask to ensure that no vectors are computed nearby - # the edges of the radar mask - if buffer_mask > 0: - mask_ = cv2.dilate( - mask_, np.ones((int(buffer_mask), int(buffer_mask)), np.uint8), 1 - ) + if ~isinstance(next_img, MaskedArray): + next_img = np.ma.masked_invalid(next_img) + np.ma.set_fill_value(next_img, next_img.min()) # remove small noise with a morphological operator (opening) if size_opening > 0: - prvs = _clean_image(prvs, n=size_opening) - next = _clean_image(next, n=size_opening) - - # Shi-Tomasi good features to track - # TODO: implement different feature detection algorithms (e.g. Harris) - mask_ = (-1 * mask_ + 1).astype("uint8") - p0 = _ShiTomasi_features_to_track( - prvs, - max_corners_ST, - quality_level_ST, - min_distance_ST, - block_size_ST, - mask_, - ) + prvs_img = morph_opening(prvs_img, prvs_img.min(), size_opening) + next_img = morph_opening(next_img, next_img.min(), size_opening) + + # features detection + points = feature_detection_method(prvs_img, **fd_kwargs) # skip loop if no features to track - if p0 is None: + if points.shape[0] == 0: continue # get sparse u, v vectors with Lucas-Kanade tracking - x0, y0, u, v = _LucasKanade_features_tracking( - prvs, next, p0, winsize_LK, nr_levels_LK - ) + xy_, uv_ = track_features(prvs_img, next_img, points, **lk_kwargs) + # skip loop if no vectors - if x0 is None: + if xy_.shape[0] == 0: continue - # stack vectors within time window as column vectors - x0Stack.append(x0.flatten()[:, None]) - y0Stack.append(y0.flatten()[:, None]) - uStack.append(u.flatten()[:, None]) - vStack.append(v.flatten()[:, None]) + # stack vectors + xy = np.append(xy, xy_, axis=0) + uv = np.append(uv, uv_, axis=0) # return zero motion field is no sparse vectors are found - if len(x0Stack) == 0: + if xy.shape[0] == 0: if dense: return np.zeros((2, domain_size[0], domain_size[1])) else: - rzero = np.array([0]) - return rzero, rzero, rzero, rzero - - # convert lists of arrays into single arrays - x = np.vstack(x0Stack) - y = np.vstack(y0Stack) - u = np.vstack(uStack) - v = np.vstack(vStack) - - # exclude outlier vectors - x, y, u, v = _outlier_removal( - x, y, u, v, nr_std_outlier, multivariate_outlier, k_outlier, verbose - ) + return xy, uv + + # detect and remove outliers + outliers = detect_outliers(uv, nr_std_outlier, xy, k_outlier, verbose) + xy = xy[~outliers, :] + uv = uv[~outliers, :] if verbose: - print("--- LK found %i sparse vectors ---" % x.size) + print("--- LK found %i sparse vectors ---" % xy.shape[0]) # return sparse vectors if required if not dense: - return x, y, u, v + return xy, uv # decluster sparse motion vectors - if decl_grid > 1: - x, y, u, v = _declustering(x, y, u, v, decl_grid, min_nr_samples) - - # append extra vectors if provided - if extra_vectors is not None: - x = np.concatenate((x, extra_vectors[:, 0])) - y = np.concatenate((y, extra_vectors[:, 1])) - u = np.concatenate((u, extra_vectors[:, 2])) - v = np.concatenate((v, extra_vectors[:, 3])) + if decl_scale > 1: + xy, uv = decluster(xy, uv, decl_scale, 1, verbose) # return zero motion field if no sparse vectors are left for interpolation - if x.size == 0: + if xy.shape[0] == 0: return np.zeros((2, domain_size[0], domain_size[1])) - if verbose: - print("--- %i sparse vectors left after declustering ---" % x.size) - - # kernel interpolation - _, _, UV = _interpolate_sparse_vectors( - x, - y, - u, - v, - domain_size, - rbfunction=rbfunction, - k=k, - epsilon=epsilon, - nchunks=nchunks, - ) + # interpolation + xgrid = np.arange(domain_size[1]) + ygrid = np.arange(domain_size[0]) + UV = interpolation_method(xy, uv, xgrid, ygrid, **interp_kwargs) if verbose: - print("--- %.2f seconds ---" % (time.time() - t0)) + print("--- total time: %.2f seconds ---" % (time.time() - t0)) return UV -def _ShiTomasi_features_to_track( - R, max_corners_ST, quality_level_ST, min_distance_ST, block_size_ST, mask +def track_features( + prvs_image, + next_image, + points, + winsize=(50, 50), + nr_levels=3, + criteria=(3, 10, 0), + flags=0, + min_eig_thr=1e-4, + verbose=False, ): - """Call the Shi-Tomasi corner detection algorithm. - - Parameters - ---------- - R : array-like - Array of shape (m,n) containing the input precipitation field passed as - 8-bit image. - - max_corners_ST : int - Maximum number of corners to return. If there are more corners than are - found, the strongest of them is returned. - - quality_level_ST : float - Parameter characterizing the minimal accepted quality of image corners. - See original documentation for more details (https://docs.opencv.org). - - min_distance_ST : int - Minimum possible Euclidean distance between the returned corners [px]. - - block_size_ST : int - Size of an average block for computing a derivative covariation matrix - over each pixel neighborhood. - - mask : ndarray_ - Array of shape (m,n). It specifies the region in which the corners are - detected. - - Returns - ------- - p0 : list - Output vector of detected corners. - """ - if not cv2_imported: - raise MissingOptionalDependency( - "opencv package is required for the Lucas-Kanade " - "optical flow method but it is not installed" - ) - - if len(R.shape) != 2: - raise ValueError("R must be a two-dimensional array") - if R.dtype != "uint8": - raise ValueError("R must be passed as 8-bit image") - - # ShiTomasi corner detection parameters - ShiTomasi_params = dict( - maxCorners=max_corners_ST, - qualityLevel=quality_level_ST, - minDistance=min_distance_ST, - blockSize=block_size_ST, - ) + Interface to the OpenCV `Lucas-Kanade`_ features tracking algorithm + (cv.calcOpticalFlowPyrLK). - # detect corners - p0 = cv2.goodFeaturesToTrack(R, mask=mask, **ShiTomasi_params) + .. _`Lucas-Kanade`:\ + https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga473e4b886d0bcc6b65831eb88ed93323 - return p0 + .. _calcOpticalFlowPyrLK:\ + https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga473e4b886d0bcc6b65831eb88ed93323 -def _LucasKanade_features_tracking(prvs, next, p0, winsize_LK, nr_levels_LK): - """Call the Lucas-Kanade features tracking algorithm. + .. _MaskedArray:\ + https://docs.scipy.org/doc/numpy/reference/maskedarray.baseclass.html#numpy.ma.MaskedArray Parameters ---------- - prvs : array-like - Array of shape (m,n) containing the first 8-bit input image. - next : array-like - Array of shape (m,n) containing the successive 8-bit input image. - p0 : list - Vector of 2D points for which the flow needs to be found. - Point coordinates must be single-precision floating-point numbers. - winsize_LK : tuple - Size of the search window at each pyramid level. - Small windows (e.g. 10) lead to unrealistic motion. - nr_levels_LK : int - 0-based maximal pyramid level number. - Not very sensitive parameter. - - Returns - ------- - x0 : array-like - Output vector of x-coordinates of detected point motions. - y0 : array-like - Output vector of y-coordinates of detected point motions. - u : array-like - Output vector of u-components of detected point motions. - v : array-like - Output vector of v-components of detected point motions. - - """ - if not cv2_imported: - raise MissingOptionalDependency( - "opencv package is required for the Lucas-Kanade method " - "optical flow method but it is not installed" - ) - # LK parameters - lk_params = dict( - winSize=winsize_LK, - maxLevel=nr_levels_LK, - criteria=(cv2.TERM_CRITERIA_EPS | cv2.TERM_CRITERIA_COUNT, 10, 0), - ) + prvs_image : array_like or MaskedArray_ + Array of shape (m, n) containing the first image. + Invalid values (Nans or infs) are filled using the min value. - # Lucas-Kande - p1, st, err = cv2.calcOpticalFlowPyrLK(prvs, next, p0, None, **lk_params) + next_image : array_like or MaskedArray_ + Array of shape (m, n) containing the successive image. + Invalid values (Nans or infs) are filled using the min value. - # keep only features that have been found - st = st[:, 0] == 1 - if np.any(st): - p1 = p1[st, :, :] - p0 = p0[st, :, :] - err = err[st, :] + points : array_like + Array of shape (p, 2) indicating the pixel coordinates of the + tracking points (corners). - # extract vectors - x0 = p0[:, :, 0] - y0 = p0[:, :, 1] - u = np.array((p1 - p0)[:, :, 0]) - v = np.array((p1 - p0)[:, :, 1]) - else: - x0 = y0 = u = v = None - - return x0, y0, u, v - - -def _clean_image(R, n=3, thr=0): - """Apply a binary morphological opening to filter small isolated echoes. - - Parameters - ---------- - R : array-like - Array of shape (m,n) containing the input precipitation field. - n : int - The structuring element size [px]. - thr : float - The rain/no-rain threshold to convert the image into a binary image. - - Returns - ------- - R : array - Array of shape (m,n) containing the cleaned precipitation field. - - """ - if not cv2_imported: - raise MissingOptionalDependency( - "opencv package is required for the Lucas-Kanade method " - "optical flow method but it is not installed" - ) - - # convert to binary image (rain/no rain) - field_bin = np.ndarray.astype(R > thr, "uint8") - - # build a structuring element of size (nx) - kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (n, n)) - - # apply morphological opening (i.e. erosion then dilation) - field_bin_out = cv2.morphologyEx(field_bin, cv2.MORPH_OPEN, kernel) - - # build mask to be applied on the original image - mask = (field_bin - field_bin_out) > 0 + winsize : tuple of int, optional + The **winSize** parameter in calcOpticalFlowPyrLK_. + It represents the size of the search window that it is used at each + pyramid level. - # filter out small isolated echoes based on mask - R[mask] = np.nanmin(R) + nr_levels : int, optional + The **maxLevel** parameter in calcOpticalFlowPyrLK_. + It represents the 0-based maximal pyramid level number. - return R + criteria : tuple of int, optional + The **TermCriteria** parameter in calcOpticalFlowPyrLK_ , + which specifies the termination criteria of the iterative search + algorithm. + flags : int, optional + Operation flags, see documentation calcOpticalFlowPyrLK_. -def _outlier_removal(x, y, u, v, thr, multivariate=True, k=30, verbose=False): + min_eig_thr : float, optional + The **minEigThreshold** parameter in calcOpticalFlowPyrLK_. - """Outlier removal. - - Parameters - ---------- - x : array_like - X-coordinates of the origins of the velocity vectors. - y : array_like - Y-coordinates of the origins of the velocity vectors. - u : array_like - X-components of the velocities. - v : array_like - Y-components of the velocities. - thr : float - Threshold for outlier detection defined as measure of deviation from - the mean/median in terms of standard deviations. - multivariate : bool, optional - If true (the default), the outlier detection is computed in terms of - the Mahalanobis distance. If false, the outlier detection is simply - computed in terms of velocity. - k : int, optinal - The number of nearest neighbours used to localize the outlier detection. - If set equal to 0, it employs all the data points. - The default is 30. + verbose : bool, optional + Print the number of vectors that have been found. Returns ------- - A four-element tuple (x,y,u,v) containing the x- and y-coordinates and - velocity components of the motion vectors. - """ - - if multivariate: - data = np.concatenate((u, v), axis=1) - - # globally - if k <= 0: - - if not multivariate: - - # in terms of velocity - - vel = np.sqrt(u ** 2 + v ** 2) # [px/timesteps] - q1, q2, q3 = np.percentile(vel, [16, 50, 84]) - min_speed_thr = np.max((0, q2 - thr * (q3 - q1) / 2)) - max_speed_thr = q2 + thr * (q3 - q1) / 2 - keep = np.logical_and(vel < max_speed_thr, vel >= min_speed_thr) - - else: - # mahalanobis distance + xy : array_like + Array of shape (d, 2) with the x- and y-coordinates of *d* <= *p* + detected sparse motion vectors. - data = data - np.mean(data, axis=0) - V = np.cov(data.T) - VI = np.linalg.inv(V) - MD = np.sqrt(np.dot(np.dot(data, VI), data.T).diagonal()) - keep = MD < thr + uv : array_like + Array of shape (d, 2) with the u- and v-components of *d* <= *p* + detected sparse motion vectors. - # locally - else: - - points = np.concatenate((x, y), axis=1) - tree = scipy.spatial.cKDTree(points) - _, inds = tree.query(points, k=np.min((k + 1, points.shape[0]))) - keep = [] - for i in range(inds.shape[0]): - - if not multivariate: - - # in terms of velocity - - thisvel = np.sqrt(u[i] ** 2 + v[i] ** 2) # [px/timesteps] - neighboursvel = np.sqrt(u[inds[i, 1:]] ** 2 + v[inds[i, 1:]] ** 2) - q1, q2, q3 = np.percentile(neighboursvel, [16, 50, 84]) - min_speed_thr = np.max((0, q2 - thr * (q3 - q1) / 2)) - max_speed_thr = q2 + thr * (q3 - q1) / 2 - keep.append(thisvel < max_speed_thr and thisvel > min_speed_thr) + Notes + ----- - else: + The tracking points can be obtained with the + :py:func:`pysteps.utils.images.ShiTomasi_detection` routine. - # mahalanobis distance + See also + -------- - thisdata = data[i, :] - neighbours = data[inds[i, 1:], :].copy() - thisdata = thisdata - np.mean(neighbours, axis=0) - neighbours = neighbours - np.mean(neighbours, axis=0) - V = np.cov(neighbours.T) - try: - VI = np.linalg.inv(V) - MD = np.sqrt(np.dot(np.dot(thisdata, VI), thisdata.T)) - - except np.linalg.LinAlgError: - MD = 0 - - keep.append(MD < thr) - - keep = np.array(keep) - - if verbose: - print("--- %i outliers removed ---" % np.sum(~keep)) + pysteps.motion.lucaskanade.dense_lucaskanade - x = x[keep] - y = y[keep] - u = u[keep] - v = v[keep] - - return x, y, u, v - - -def _declustering(x, y, u, v, decl_grid, min_nr_samples): - """Filter out outliers in a sparse motion field and get more representative - data points. The method assigns data points to a (RxR) declustering grid - and then take the median of all values within one cell. - - Parameters + References ---------- - x : array_like - X-coordinates of the origins of the velocity vectors. - y : array_like - Y-coordinates of the origins of the velocity vectors. - u : array_like - X-components of the velocities. - v : array_like - Y-components of the velocities. - decl_grid : int - Size of the declustering grid [px]. - min_nr_samples : int - The minimum number of samples for computing the median within given - declustering cell. - Returns - ------- - A four-element tuple (x,y,u,v) containing the x- and y-coordinates and - velocity components of the declustered motion vectors. + Bouguet, J.-Y.: Pyramidal implementation of the affine Lucas Kanade + feature tracker description of the algorithm, Intel Corp., 5, 4, + https://doi.org/10.1109/HPDC.2004.1323531, 2001 + Lucas, B. D. and Kanade, T.: An iterative image registration technique with + an application to stereo vision, in: Proceedings of the 1981 DARPA Imaging + Understanding Workshop, pp. 121–130, 1981. """ - # make sure these are all numpy vertical arrays - x = np.array(x).flatten()[:, None] - y = np.array(y).flatten()[:, None] - u = np.array(u).flatten()[:, None] - v = np.array(v).flatten()[:, None] - - # return empty arrays if the number of sparse vectors is < min_nr_samples - if x.size < min_nr_samples: - return np.array([]), np.array([]), np.array([]), np.array([]) - - # discretize coordinates into declustering grid - xT = x / float(decl_grid) - yT = y / float(decl_grid) - - # round coordinates to low integer - xT = np.floor(xT) - yT = np.floor(yT) + if not CV2_IMPORTED: + raise MissingOptionalDependency( + "opencv package is required for the calcOpticalFlowPyrLK() " + "routine but it is not installed" + ) - # keep only unique combinations of coordinates - xy = np.concatenate((xT, yT), axis=1) - xyb = np.ascontiguousarray(xy).view( - np.dtype((np.void, xy.dtype.itemsize * xy.shape[1])) + prvs_img = np.copy(prvs_image) + next_img = np.copy(next_image) + p0 = np.copy(points) + + if ~isinstance(prvs_img, MaskedArray): + prvs_img = np.ma.masked_invalid(prvs_img) + np.ma.set_fill_value(prvs_img, prvs_img.min()) + + if ~isinstance(next_img, MaskedArray): + next_img = np.ma.masked_invalid(next_img) + np.ma.set_fill_value(next_img, next_img.min()) + + # scale between 0 and 255 + prvs_img = ((prvs_img.filled() - prvs_img.min()) / + (prvs_img.max() - prvs_img.min()) * 255) + + next_img = ((next_img.filled() - next_img.min()) / + (next_img.max() - next_img.min()) * 255) + + # convert to 8-bit + prvs_img = np.ndarray.astype(prvs_img, "uint8") + next_img = np.ndarray.astype(next_img, "uint8") + + # Lucas-Kanade + # TODO: use the error returned by the OpenCV routine + params = dict( + winSize=winsize, + maxLevel=nr_levels, + criteria=criteria, + flags=flags, + minEigThreshold=min_eig_thr, ) - _, idx = np.unique(xyb, return_index=True) - uxy = xy[idx] - - # now loop through these unique values and average vectors which belong to - # the same declustering grid cell - xN = [] - yN = [] - uN = [] - vN = [] - for i in range(uxy.shape[0]): - idx = np.logical_and(xT == uxy[i, 0], yT == uxy[i, 1]) - npoints = np.sum(idx) - if npoints >= min_nr_samples: - xN.append(np.median(x[idx])) - yN.append(np.median(y[idx])) - uN.append(np.median(u[idx])) - vN.append(np.median(v[idx])) - - # convert to numpy arrays - x = np.array(xN) - y = np.array(yN) - u = np.array(uN) - v = np.array(vN) - - return x, y, u, v - - -def _interpolate_sparse_vectors( - x, y, u, v, domain_size, rbfunction="inverse", k=20, epsilon=None, nchunks=5 -): - - """Interpolation of sparse motion vectors to produce a dense field of motion - vectors. + p1, st, __ = cv2.calcOpticalFlowPyrLK(prvs_img, next_img, + p0, None, **params) - Parameters - ---------- - x : array-like - x coordinates of the sparse motion vectors - y : array-like - y coordinates of the sparse motion vectors - u : array_like - u components of the sparse motion vectors - v : array_like - v components of the sparse motion vectors - domain_size : tuple - size of the domain of the dense motion field [px] - rbfunction : string - the radial basis rbfunction, based on the Euclidian norm, d. - default : inverse - available : nearest, inverse, gaussian - k : int or "all" - the number of nearest neighbours used to speed-up the interpolation - If set equal to "all", it employs all the sparse vectors - default : 20 - epsilon : float - adjustable constant for gaussian or inverse functions - default : median distance between sparse vectors - nchunks : int - split the grid points in n chunks to limit the memory usage during the - interpolation - default : 5 - - Returns - ------- - X : array-like - grid - Y : array-like - grid - UV : array-like - Three-dimensional array (2,domain_size[0],domain_size[1]) - containing the dense U, V motion fields. - - """ - - # make sure these are vertical arrays - x = np.array(x).flatten()[:, None] - y = np.array(y).flatten()[:, None] - u = np.array(u).flatten()[:, None] - v = np.array(v).flatten()[:, None] - points = np.concatenate((x, y), axis=1) - npoints = points.shape[0] - - if len(domain_size) == 1: - domain_size = (domain_size, domain_size) - - # generate the grid - xgrid = np.arange(domain_size[1]) - ygrid = np.arange(domain_size[0]) - X, Y = np.meshgrid(xgrid, ygrid) - grid = np.column_stack((X.ravel(), Y.ravel())) - - U = np.zeros(grid.shape[0]) - V = np.zeros(grid.shape[0]) + # keep only features that have been found + st = st.squeeze() == 1 + if np.any(st): + p1 = p1[st, :] + p0 = p0[st, :] - # create cKDTree object to represent source grid - if k > 0: - k = np.min((k, npoints)) - tree = scipy.spatial.cKDTree(points) + # extract vectors + xy = p0 + uv = p1 - p0 - # split grid points in n chunks - if nchunks > 1: - subgrids = np.array_split(grid, nchunks, 0) - subgrids = [x for x in subgrids if x.size > 0] else: - subgrids = [grid] - - # loop subgrids - i0 = 0 - for i, subgrid in enumerate(subgrids): + xy = uv = np.empty(shape=(0, 2)) - idelta = subgrid.shape[0] - - if rbfunction.lower() == "nearest": - # find indices of the nearest neighbours - _, inds = tree.query(subgrid, k=1) - - U[i0 : (i0 + idelta)] = u.ravel()[inds] - V[i0 : (i0 + idelta)] = v.ravel()[inds] + if verbose: + print("--- %i sparse vectors found ---" % xy.shape[0]) - else: - if k <= 0: - d = scipy.spatial.distance.cdist( - points, subgrid, "euclidean" - ).transpose() - inds = np.arange(u.size)[None, :] * np.ones( - (subgrid.shape[0], u.size) - ).astype(int) - - else: - # find indices of the k-nearest neighbours - d, inds = tree.query(subgrid, k=k) - - if inds.ndim == 1: - inds = inds[:, None] - d = d[:, None] - - # the bandwidth - if epsilon is None: - epsilon = 1 - if npoints > 1: - dpoints = scipy.spatial.distance.pdist(points, "euclidean") - epsilon = np.median(dpoints) - - # the interpolation weights - if rbfunction.lower() == "inverse": - w = 1.0 / np.sqrt((d / epsilon) ** 2 + 1) - elif rbfunction.lower() == "gaussian": - w = np.exp(-0.5 * (d / epsilon) ** 2) - else: - raise ValueError("unknown radial fucntion %s" % rbfunction) - - if not np.all(np.sum(w, axis=1)): - w[np.sum(w, axis=1) == 0, :] = 1.0 - - U[i0 : (i0 + idelta)] = np.sum(w * u.ravel()[inds], axis=1) / np.sum( - w, axis=1 - ) - V[i0 : (i0 + idelta)] = np.sum(w * v.ravel()[inds], axis=1) / np.sum( - w, axis=1 - ) - - i0 += idelta - - # reshape back to original size - U = U.reshape(domain_size[0], domain_size[1]) - V = V.reshape(domain_size[0], domain_size[1]) - UV = np.stack([U, V]) - - return X, Y, UV + return xy, uv diff --git a/pysteps/motion/proesmans.py b/pysteps/motion/proesmans.py index feb46d35f..8e987b286 100644 --- a/pysteps/motion/proesmans.py +++ b/pysteps/motion/proesmans.py @@ -12,9 +12,14 @@ import numpy as np from scipy.ndimage import gaussian_filter + +from pysteps.decorators import check_input_frames from pysteps.motion._proesmans import _compute_advection_field -def proesmans(input_images, lam=50.0, num_iter=100, num_levels=6, filter_std=0.0): + +@check_input_frames(2, 2) +def proesmans(input_images, lam=50.0, num_iter=100, + num_levels=6, filter_std=0.0, verbose=True, ): """Implementation of the anisotropic diffusion method of Proesmans et al. (1994). @@ -32,6 +37,8 @@ def proesmans(input_images, lam=50.0, num_iter=100, num_levels=6, filter_std=0.0 filter_std : float Standard deviation of an optional Gaussian filter that is applied before computing the optical flow. + verbose : bool, optional + Verbosity enabled if True (default). Returns ------- @@ -46,18 +53,16 @@ def proesmans(input_images, lam=50.0, num_iter=100, num_levels=6, filter_std=0.0 :cite:`PGPO1994` """ - if (input_images.ndim != 3) or input_images.shape[0] != 2: - raise ValueError("input_images dimension mismatch.\n" + - "input_images.shape: " + str(input_images.shape) + - "\n(2, m, n) expected") - + del verbose # Not used + im1 = input_images[-2, :, :].copy() im2 = input_images[-1, :, :].copy() im = np.stack([im1, im2]) im_min = np.min(im) im_max = np.max(im) - im = (im - im_min) / (im_max - im_min) * 255.0 + if im_max - im_min > 1e-8: + im = (im - im_min) / (im_max - im_min) * 255.0 if filter_std > 0.0: im[0, :, :] = gaussian_filter(im[0, :, :], filter_std) diff --git a/pysteps/motion/vet.py b/pysteps/motion/vet.py index 6474af3db..a3363540a 100644 --- a/pysteps/motion/vet.py +++ b/pysteps/motion/vet.py @@ -39,6 +39,7 @@ from scipy.ndimage.interpolation import zoom from scipy.optimize import minimize +from pysteps.decorators import check_input_frames from pysteps.motion._vet import _warp, _cost_function @@ -314,6 +315,7 @@ def vet_cost_function(sector_displacement_1d, return residuals + smoothness_penalty +@check_input_frames(2, 3) def vet(input_images, sectors=((32, 16, 4, 2), (32, 16, 4, 2)), smooth_gain=1e6, @@ -488,15 +490,10 @@ def debug_print(*args, **kwargs): debug_print("Running VET algorithm") - if (input_images.ndim != 3) or (1 < input_images.shape[0] > 3): - raise ValueError("input_images dimension mismatch.\n" + - "input_images.shape: " + str(input_images.shape) + - "\n(2, x, y ) or (2, x, y ) dimensions expected") - valid_indexing = ['yx', 'xy', 'ij'] if indexing not in valid_indexing: - raise ValueError("Invalid indexing valus: {0}\n".format(indexing) + raise ValueError("Invalid indexing values: {0}\n".format(indexing) + "Supported values: {0}".format(str(valid_indexing))) # Get mask @@ -570,7 +567,8 @@ def debug_print(*args, **kwargs): if (pad_i != (0, 0)) or (pad_j != (0, 0)): - _input_images = numpy.pad(input_images, ((0, 0), pad_i, pad_j), 'edge') + _input_images = numpy.pad(input_images, ((0, 0), pad_i, pad_j), + 'edge') _mask = numpy.pad(mask, (pad_i, pad_j), 'constant', diff --git a/pysteps/tests/test_interfaces.py b/pysteps/tests/test_interfaces.py index 76177c80c..583fd7cd6 100644 --- a/pysteps/tests/test_interfaces.py +++ b/pysteps/tests/test_interfaces.py @@ -206,19 +206,35 @@ def test_nowcasts_interface(): def test_utils_interface(): """Test utils module interface.""" + from pysteps.utils import arrays + from pysteps.utils import cleansing from pysteps.utils import conversion - from pysteps.utils import transformation from pysteps.utils import dimension + from pysteps.utils import images + from pysteps.utils import interpolate + from pysteps.utils import spectral + from pysteps.utils import transformation method_getter = pysteps.utils.interface.get_method - valid_names_func_pair = [('mm/h', conversion.to_rainrate), + valid_names_func_pair = [('centred_coord', arrays.compute_centred_coord_array), + ('decluster', cleansing.decluster), + ('detect_outliers', cleansing.detect_outliers), + ('mm/h', conversion.to_rainrate), ('rainrate', conversion.to_rainrate), ('mm', conversion.to_raindepth), ('raindepth', conversion.to_raindepth), ('dbz', conversion.to_reflectivity), ('reflectivity', conversion.to_reflectivity), - ('rainrate', conversion.to_rainrate), + ('accumulate', dimension.aggregate_fields_time), + ('clip', dimension.clip_domain), + ('square', dimension.square_domain), + ('upscale', dimension.aggregate_fields_space), + ('shitomasi', images.ShiTomasi_detection), + ('morph_opening', images.morph_opening), + ('rbfinterp2d', interpolate.rbfinterp2d), + ('rapsd', spectral.rapsd), + ('rm_rdisc', spectral.remove_rain_norain_discontinuity), ('boxcox', transformation.boxcox_transform), ('box-cox', transformation.boxcox_transform), ('db', transformation.dB_transform), @@ -226,10 +242,6 @@ def test_utils_interface(): ('log', transformation.boxcox_transform), ('nqt', transformation.NQ_transform), ('sqrt', transformation.sqrt_transform), - ('accumulate', dimension.aggregate_fields_time), - ('clip', dimension.clip_domain), - ('square', dimension.square_domain), - ('upscale', dimension.aggregate_fields_space), ] invalid_names = ['random', 'invalid'] diff --git a/pysteps/tests/test_motion.py b/pysteps/tests/test_motion.py index 84b8f0d09..8d3d366e4 100644 --- a/pysteps/tests/test_motion.py +++ b/pysteps/tests/test_motion.py @@ -17,6 +17,8 @@ the retrieval. """ +from contextlib import contextmanager + import numpy as np import pytest from scipy.ndimage import uniform_filter @@ -26,6 +28,15 @@ from pysteps.motion.vet import morph from pysteps.tests.helpers import get_precipitation_fields + +@contextmanager +def not_raises(_exception): + try: + yield + except _exception: + raise pytest.fail("DID RAISE {0}".format(_exception)) + + reference_field = get_precipitation_fields(num_prev_files=0) @@ -226,9 +237,13 @@ def test_optflow_method_convergence(input_precip, optflow_method_name, no_precip_args_names = ("optflow_method_name, num_times") -no_precip_args_values = [('lk', 2), ('lk', 3), - ('vet', 2), ('vet', 3), - ('darts', 9)] +no_precip_args_values = [('lk', 2), + ('lk', 3), + ('vet', 2), + ('vet', 3), + ('darts', 9), + ('proesmans', 2) + ] @pytest.mark.parametrize(no_precip_args_names, no_precip_args_values) @@ -256,9 +271,46 @@ def test_no_precipitation(optflow_method_name, num_times): assert np.abs(uv_motion).max() < 0.01 +input_tests_args_names = ("optflow_method_name", + "minimum_input_frames", + "maximum_input_frames") +input_tests_args_values = [ + ('lk', 2, np.inf), + ('vet', 2, 3), + ('darts', 9, 9), + ('proesmans', 2, 2), +] + + +@pytest.mark.parametrize(input_tests_args_names, input_tests_args_values) +def test_input_shape_checks(optflow_method_name, + minimum_input_frames, + maximum_input_frames): + image_size = 100 + motion_method = motion.get_method(optflow_method_name) + + if maximum_input_frames == np.inf: + maximum_input_frames = minimum_input_frames + 10 + + with not_raises(Exception): + for frames in range(minimum_input_frames, maximum_input_frames + 1): + motion_method(np.zeros((frames, image_size, image_size)), verbose=False) + + with pytest.raises(ValueError): + motion_method(np.zeros((2,))) + motion_method(np.zeros((2, 2))) + for frames in range(minimum_input_frames): + motion_method(np.zeros((frames, image_size, image_size)), + verbose=False) + for frames in range(maximum_input_frames + 1, maximum_input_frames + 4): + motion_method(np.zeros((frames, image_size, image_size)), + verbose=False) + + def test_vet_cost_function(): """ - Test that the vet cost_function computation gives always the same result with the same input. + Test that the vet cost_function computation gives always the same result + with the same input. Useful to test if the parallelization in VET produce undesired results. """ @@ -281,7 +333,6 @@ def test_vet_cost_function(): 1e6, # smooth_gain debug=False) - print(returned_values) tolerance = 1e-12 errors = np.abs(returned_values - returned_values[0]) # errors should contain all zeros diff --git a/pysteps/tests/test_nowcasts_steps.py b/pysteps/tests/test_nowcasts_steps.py index bf7f46716..bf4de1ffd 100644 --- a/pysteps/tests/test_nowcasts_steps.py +++ b/pysteps/tests/test_nowcasts_steps.py @@ -64,11 +64,11 @@ def _import_mch_gif(prv, nxt): ) steps_arg_values = [ - (5, 6, 2, None, None, 1.51), - (5, 6, 2, "incremental", None, 6.38), - (5, 6, 2, "sprog", None, 7.35), - (5, 6, 2, "obs", None, 7.36), - (5, 6, 2, None, "cdf", 0.66), + (5, 6, 2, None, None, 1.55), + (5, 6, 2, "incremental", None, 6.65), + (5, 6, 2, "sprog", None, 7.65), + (5, 6, 2, "obs", None, 7.65), + (5, 6, 2, None, "cdf", 0.70), (5, 6, 2, None, "mean", 1.55), ] diff --git a/pysteps/utils/__init__.py b/pysteps/utils/__init__.py index 07752892d..283ad3883 100644 --- a/pysteps/utils/__init__.py +++ b/pysteps/utils/__init__.py @@ -1,9 +1,12 @@ """Miscellaneous utility functions.""" from .arrays import * +from .cleansing import * from .conversion import * from .dimension import * +from .images import * from .interface import get_method +from .interpolate import * from .fft import * from .spectral import * from .transformation import * \ No newline at end of file diff --git a/pysteps/utils/cleansing.py b/pysteps/utils/cleansing.py new file mode 100644 index 000000000..4c38cdfc5 --- /dev/null +++ b/pysteps/utils/cleansing.py @@ -0,0 +1,283 @@ +# -*- coding: utf-8 -*- +""" +pysteps.utils.cleansing +======================= + +Data cleansing routines for pysteps. + +.. autosummary:: + :toctree: ../generated/ + + decluster + detect_outliers +""" + +import numpy as np +import scipy.spatial + + +def decluster(coord, input_array, scale, min_samples=1, verbose=False): + """Decluster a set of sparse data points by aggregating, that is, taking + the median value of all values lying within a certain distance (i.e., a + cluster). + + Parameters + ---------- + + coord : array_like + Array of shape (n, d) containing the coordinates of the input data into + a space of *d* dimensions. + + input_array : array_like + Array of shape (n) or (n, m), where *n* is the number of samples and + *m* the number of variables. + All values in **input_array** are required to have finite values. + + scale : float or array_like + The **scale** parameter in the same units of **coord**. + It can be a scalar or an array_like of shape (d). + Data points within the declustering **scale** are aggregated. + + min_samples : int, optional + The minimum number of samples for computing the median within a given + cluster. + + verbose : bool, optional + Print out information. + + Returns + ------- + + out : tuple of ndarrays + A two-element tuple (**out_coord**, **output_array**) containing the + declustered coordinates (l, d) and **input_array** (l, m), where *l* is + the new number of samples with *l* <= *n*. + """ + + coord = np.copy(coord) + input_array = np.copy(input_array) + + # check inputs + if np.any(~np.isfinite(input_array)): + raise ValueError("input_array contains non-finite values") + + if input_array.ndim == 1: + nvar = 1 + input_array = input_array[:, None] + elif input_array.ndim == 2: + nvar = input_array.shape[1] + else: + raise ValueError( + "input_array must have 1 (n) or 2 dimensions (n, m), but it has %i" + % input_array.ndim + ) + + if coord.ndim != 2: + raise ValueError( + "coord must have 2 dimensions (n, d), but it has %i" % coord.ndim + ) + if coord.shape[0] != input_array.shape[0]: + raise ValueError( + "the number of samples in the input_array does not match the " + + "number of coordinates %i!=%i" + % (input_array.shape[0], coord.shape[0]) + ) + + if np.isscalar(scale): + scale = np.float(scale) + else: + scale = np.copy(scale) + if scale.ndim != 1: + raise ValueError( + "scale must have 1 dimension (d), but it has %i" % scale.ndim + ) + if scale.shape[0] != coord.shape[1]: + raise ValueError( + "scale must have %i elements, but it has %i" + % (coord.shape[1], scale.shape[0]) + ) + scale = scale[None, :] + + # reduce original coordinates + coord_ = np.floor(coord / scale) + + # keep only unique pairs of the reduced coordinates + coordb_ = np.ascontiguousarray(coord_).view( + np.dtype((np.void, coord_.dtype.itemsize * coord_.shape[1])) + ) + __, idx = np.unique(coordb_, return_index=True) + ucoord_ = coord_[idx] + # TODO: why not simply using np.unique(coord_, axis=0) ? + + # loop through these unique values and average data points which belong to + # the same cluster + dinput = np.empty(shape=(0, nvar)) + dcoord = np.empty(shape=(0, coord.shape[1])) + for i in range(ucoord_.shape[0]): + idx = np.all(coord_ == ucoord_[i, :], axis=1) + npoints = np.sum(idx) + if npoints >= min_samples: + dinput = np.append( + dinput, np.median(input_array[idx, :], axis=0)[None, :], axis=0 + ) + dcoord = np.append( + dcoord, np.median(coord[idx, :], axis=0)[None, :], axis=0 + ) + + if verbose: + print("--- %i samples left after declustering ---" % dinput.shape[0]) + + return dcoord.squeeze(), dinput + + +def detect_outliers(input_array, thr, coord=None, k=None, verbose=False): + """Detect outliers in a (multivariate and georeferenced) dataset. + + Assume a (multivariate) Gaussian distribution and detect outliers based on + the number of standard deviations from the mean. + + If spatial information is provided through coordinates, the outlier + detection can be localized by considering only the k-nearest neighbours + when computing the local mean and standard deviation. + + Parameters + ---------- + + input_array : array_like + Array of shape (n) or (n, m), where *n* is the number of samples and + *m* the number of variables. If *m* > 1, the Mahalanobis distance + is used. + All values in **input_array** are required to have finite values. + + thr : float + The number of standard deviations from the mean that defines an outlier. + + coord : array_like, optional + Array of shape (n, d) containing the coordinates of the input data into + a space of *d* dimensions. + Passing **coord** requires that **k** is not None. + + k : int or None, optional + The number of nearest neighbours used to localize the outlier detection. + If set to None (the default), it employs all the data points (global + detection). Setting **k** requires that **coord** is not None. + + verbose : bool, optional + Print out information. + + Returns + ------- + + out : array_like + A boolean array of the same shape as **input_array**, with True values + indicating the outliers detected in **input_array**. + """ + + input_array = np.copy(input_array) + + if np.any(~np.isfinite(input_array)): + raise ValueError("input_array contains non-finite values") + + if input_array.ndim == 1: + nvar = 1 + elif input_array.ndim == 2: + nvar = input_array.shape[1] + else: + raise ValueError( + "input_array must have 1 (n) or 2 dimensions (n, m), but it has %i" + % coord.ndim + ) + + if coord is not None: + + coord = np.copy(coord) + if coord.ndim == 1: + coord = coord[:, None] + + elif coord.ndim > 2: + raise ValueError( + "coord must have 2 dimensions (n, d), but it has %i" + % coord.ndim + ) + + if coord.shape[0] != input_array.shape[0]: + raise ValueError( + "the number of samples in input_array does not match the " + + "number of coordinates %i!=%i" + % (input_array.shape[0], coord.shape[0]) + ) + + if k is None: + raise ValueError("coord is set but k is None") + + k = np.min((coord.shape[0], k + 1)) + + else: + if k is not None: + raise ValueError("k is set but coord=None") + + # global + + if k is None: + + if nvar == 1: + + # univariate + + zdata = (input_array - np.mean(input_array)) / np.std(input_array) + outliers = zdata > thr + + else: + + # multivariate (mahalanobis distance) + + zdata = input_array - np.mean(input_array, axis=0) + V = np.cov(zdata.T) + VI = np.linalg.inv(V) + try: + VI = np.linalg.inv(V) + MD = np.sqrt(np.dot(np.dot(zdata, VI), zdata.T).diagonal()) + except np.linalg.LinAlgError: + MD = np.zeros(input_array.shape) + outliers = MD > thr + + # local + + else: + + tree = scipy.spatial.cKDTree(coord) + __, inds = tree.query(coord, k=k) + outliers = np.empty(shape=0, dtype=bool) + for i in range(inds.shape[0]): + + if nvar == 1: + + # univariate + + thisdata = input_array[i] + neighbours = input_array[inds[i, 1:]] + thiszdata = (thisdata - np.mean(neighbours)) / np.std( + neighbours + ) + outliers = np.append(outliers, thiszdata > thr) + + else: + + # multivariate (mahalanobis distance) + + thisdata = input_array[i, :] + neighbours = input_array[inds[i, 1:], :].copy() + thiszdata = thisdata - np.mean(neighbours, axis=0) + neighbours = neighbours - np.mean(neighbours, axis=0) + V = np.cov(neighbours.T) + try: + VI = np.linalg.inv(V) + MD = np.sqrt(np.dot(np.dot(thiszdata, VI), thiszdata.T)) + except np.linalg.LinAlgError: + MD = 0 + outliers = np.append(outliers, MD > thr) + + if verbose: + print("--- %i outliers detected ---" % np.sum(outliers)) + + return outliers diff --git a/pysteps/utils/images.py b/pysteps/utils/images.py new file mode 100644 index 000000000..b0972808a --- /dev/null +++ b/pysteps/utils/images.py @@ -0,0 +1,214 @@ +# -*- coding: utf-8 -*- +""" +pysteps.utils.images +==================== + +Image processing routines for pysteps. + +.. _`Shi-Tomasi`:\ + https://docs.opencv.org/3.4.1/dd/d1a/group__imgproc__feature.html#ga1d6bb77486c8f92d79c8793ad995d541 + + +.. autosummary:: + :toctree: ../generated/ + + ShiTomasi_detection + morph_opening +""" + +import numpy as np +from numpy.ma.core import MaskedArray + +try: + import cv2 + + CV2_IMPORTED = True +except ImportError: + CV2_IMPORTED = False + + +def ShiTomasi_detection(input_image, max_corners=500, quality_level=0.1, + min_distance=3, block_size=15, buffer_mask=0, + use_harris = False, k = 0.04, + verbose=False, + **kwargs): + """ + Interface to the OpenCV `Shi-Tomasi`_ features detection method to detect + corners in an image. + + Corners are used for local tracking methods. + + .. _`Shi-Tomasi`:\ + https://docs.opencv.org/3.4.1/dd/d1a/group__imgproc__feature.html#ga1d6bb77486c8f92d79c8793ad995d541 + + .. _MaskedArray:\ + https://docs.scipy.org/doc/numpy/reference/maskedarray.baseclass.html#numpy.ma.MaskedArray + + .. _`Harris detector`:\ + https://docs.opencv.org/3.4.1/dd/d1a/group__imgproc__feature.html#gac1fc3598018010880e370e2f709b4345 + + .. _cornerMinEigenVal:\ + https://docs.opencv.org/3.4.1/dd/d1a/group__imgproc__feature.html#ga3dbce297c1feb859ee36707e1003e0a8 + + + Parameters + ---------- + + input_image : array_like or MaskedArray_ + Array of shape (m, n) containing the input image. + + In case of array_like, invalid values (Nans or infs) are masked, + otherwise the mask of the MaskedArray_ is used. Such mask defines a + region where features are not detected. + + The fill value for the masked pixels is taken as the minimum of all + valid pixels. + + max_corners : int, optional + The **maxCorners** parameter in the `Shi-Tomasi`_ corner detection + method. + It represents the maximum number of points to be tracked (corners). + If set to zero, all detected corners are used. + + quality_level : float, optional + The **qualityLevel** parameter in the `Shi-Tomasi`_ corner detection + method. + It represents the minimal accepted quality for the image corners. + + min_distance : int, optional + The **minDistance** parameter in the `Shi-Tomasi`_ corner detection + method. + It represents minimum possible Euclidean distance in pixels between + corners. + + block_size : int, optional + The **blockSize** parameter in the `Shi-Tomasi`_ corner detection method. + It represents the window size in pixels used for computing a derivative + covariation matrix over each pixel neighborhood. + + use_harris : bool, optional + Whether to use a `Harris detector`_ or cornerMinEigenVal_. + + k : float, optional + Free parameter of the Harris detector. + + buffer_mask : int, optional + A mask buffer width in pixels. This extends the input mask (if any) + to limit edge effects. + + verbose : bool, optional + Print the number of features detected. + + Returns + ------- + + points : array_like + Array of shape (p, 2) indicating the pixel coordinates of *p* detected + corners. + + References + ---------- + + Jianbo Shi and Carlo Tomasi. Good features to track. In Computer Vision and + Pattern Recognition, 1994. Proceedings CVPR'94., 1994 IEEE Computer Society + Conference on, pages 593–600. IEEE, 1994. + """ + if not CV2_IMPORTED: + raise MissingOptionalDependency( + "opencv package is required for the goodFeaturesToTrack() " + "routine but it is not installed" + ) + + input_image = np.copy(input_image) + + if input_image.ndim != 2: + raise ValueError("input_image must be a two-dimensional array") + + # masked array + if ~isinstance(input_image, MaskedArray): + input_image = np.ma.masked_invalid(input_image) + np.ma.set_fill_value(input_image, input_image.min()) + + # buffer the quality mask to ensure that no vectors are computed nearby + # the edges of the radar mask + mask = np.ma.getmaskarray(input_image).astype("uint8") + if buffer_mask > 0: + mask = cv2.dilate( + mask, np.ones((int(buffer_mask), int(buffer_mask)), np.uint8), 1 + ) + input_image[mask] = np.ma.masked + + # scale image between 0 and 255 + input_image = ( + (input_image.filled() - input_image.min()) + / (input_image.max() - input_image.min()) + * 255 + ) + + # convert to 8-bit + input_image = np.ndarray.astype(input_image, "uint8") + mask = (-1 * mask + 1).astype("uint8") + + params = dict( + maxCorners=max_corners, + qualityLevel=quality_level, + minDistance=min_distance, + useHarrisDetector=use_harris, + k=k, + ) + points = cv2.goodFeaturesToTrack(input_image, mask=mask, **params) + if points is None: + points = np.empty(shape=(0, 2)) + else: + points = points.squeeze() + + if verbose: + print("--- %i good features to track detected ---" % points.shape[0]) + + return points + + +def morph_opening(input_image, thr, n): + """Filter out small scale noise on the image by applying a binary + morphological opening, that is, erosion followed by dilation. + + Parameters + ---------- + + input_image : array_like + Array of shape (m, n) containing the input image. + + thr : float + The threshold used to convert the image into a binary image. + + n : int + The structuring element size [pixels]. + + Returns + ------- + + input_image : array_like + Array of shape (m,n) containing the filtered image. + """ + if not CV2_IMPORTED: + raise MissingOptionalDependency( + "opencv package is required for the morphologyEx " + "routine but it is not installed" + ) + + # Convert to binary image + field_bin = np.ndarray.astype(input_image > thr, "uint8") + + # Build a structuring element of size n + kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (n, n)) + + # Apply morphological opening (i.e. erosion then dilation) + field_bin_out = cv2.morphologyEx(field_bin, cv2.MORPH_OPEN, kernel) + + # Build mask to be applied on the original image + mask = (field_bin - field_bin_out) > 0 + + # Filter out small isolated pixels based on mask + input_image[mask] = np.nanmin(input_image) + + return input_image diff --git a/pysteps/utils/interface.py b/pysteps/utils/interface.py index 10e0ac726..fcbe00096 100644 --- a/pysteps/utils/interface.py +++ b/pysteps/utils/interface.py @@ -11,11 +11,15 @@ """ from . import arrays +from . import cleansing from . import conversion -from . import transformation from . import dimension from . import fft +from . import images +from . import interpolate from . import spectral +from . import transformation + def get_method(name, **kwargs): """Return a callable function for the utility method corresponding to the @@ -23,66 +27,79 @@ def get_method(name, **kwargs): Arrays methods: - +-------------------+--------------------------------------------------------+ - | Name | Description | - +===================+========================================================+ - | centred_coord | compute a 2D coordinate array | - +-------------------+--------------------------------------------------------+ + +-------------------+-----------------------------------------------------+ + | Name | Description | + +===================+=====================================================+ + | centred_coord | compute a 2D coordinate array | + +-------------------+-----------------------------------------------------+ - Conversion methods: + Cleansing methods: - +-------------------+--------------------------------------------------------+ - | Name | Description | - +===================+========================================================+ - | mm/h or rainrate | convert to rain rate [mm/h] | - +-------------------+--------------------------------------------------------+ - | mm or raindepth | convert to rain depth [mm] | - +-------------------+--------------------------------------------------------+ - | dbz or | convert to reflectivity [dBZ] | - | reflectivity | | - +-------------------+--------------------------------------------------------+ + +-------------------+-----------------------------------------------------+ + | Name | Description | + +===================+=====================================================+ + | decluster | decluster a set of sparse data points | + +-------------------+-----------------------------------------------------+ + | detect_outliers | detect outliers in a dataset | + +-------------------+-----------------------------------------------------+ - Transformation methods: + Conversion methods: - +-------------------+--------------------------------------------------------+ - | Name | Description | - +===================+========================================================+ - | boxcox or box-cox | one-parameter Box-Cox transform | - +-------------------+--------------------------------------------------------+ - | db or decibel | transform to units of decibel | - +-------------------+--------------------------------------------------------+ - | log | log transform | - +-------------------+--------------------------------------------------------+ - | nqt | Normal Quantile Transform | - +-------------------+--------------------------------------------------------+ - | sqrt | square-root transform | - +-------------------+--------------------------------------------------------+ + +-------------------+-----------------------------------------------------+ + | Name | Description | + +===================+=====================================================+ + | mm/h or rainrate | convert to rain rate [mm/h] | + +-------------------+-----------------------------------------------------+ + | mm or raindepth | convert to rain depth [mm] | + +-------------------+-----------------------------------------------------+ + | dbz or | convert to reflectivity [dBZ] | + | reflectivity | | + +-------------------+-----------------------------------------------------+ Dimension methods: - +-------------------+--------------------------------------------------------+ - | Name | Description | - +===================+========================================================+ - | accumulate | aggregate fields in time | - +-------------------+--------------------------------------------------------+ - | clip | resize the field domain by geographical coordinates | - +-------------------+--------------------------------------------------------+ - | square | either pad or crop the data to get a square domain | - +-------------------+--------------------------------------------------------+ - | upscale | upscale the field | - +-------------------+--------------------------------------------------------+ + +-------------------+-----------------------------------------------------+ + | Name | Description | + +===================+=====================================================+ + | accumulate | aggregate fields in time | + +-------------------+-----------------------------------------------------+ + | clip | resize the field domain by geographical coordinates | + +-------------------+-----------------------------------------------------+ + | square | either pad or crop the data to get a square domain | + +-------------------+-----------------------------------------------------+ + | upscale | upscale the field | + +-------------------+-----------------------------------------------------+ FFT methods (wrappers to different implementations): - +-------------------+--------------------------------------------------------+ - | Name | Description | - +===================+========================================================+ - | numpy | numpy.fft | - +-------------------+--------------------------------------------------------+ - | scipy | scipy.fftpack | - +-------------------+--------------------------------------------------------+ - | pyfftw | pyfftw.interfaces.numpy_fft | - +-------------------+--------------------------------------------------------+ + +-------------------+-----------------------------------------------------+ + | Name | Description | + +===================+=====================================================+ + | numpy | numpy.fft | + +-------------------+-----------------------------------------------------+ + | scipy | scipy.fftpack | + +-------------------+-----------------------------------------------------+ + | pyfftw | pyfftw.interfaces.numpy_fft | + +-------------------+-----------------------------------------------------+ + + Image processing methods: + + +-------------------+-----------------------------------------------------+ + | Name | Description | + +===================+=====================================================+ + | ShiTomasi | Shi-Tomasi corner detection on an image | + +-------------------+-----------------------------------------------------+ + | morph_opening | filter small scale noise on an image | + +-------------------+-----------------------------------------------------+ + + Interpolation methods: + + +-------------------+-----------------------------------------------------+ + | Name | Description | + +===================+=====================================================+ + | rbfinterp2d | fast kernel interpolation of a (multivariate) array | + | | over a 2D grid using a radial basis function | + +-------------------+-----------------------------------------------------+ Additional keyword arguments are passed to the initializer of the FFT methods, see utils.fft. @@ -97,6 +114,22 @@ def get_method(name, **kwargs): | rm_rdisc | remove the rain / no-rain discontinuity | +-------------------+-----------------------------------------------------+ + Transformation methods: + + +-------------------+-----------------------------------------------------+ + | Name | Description | + +===================+=====================================================+ + | boxcox or box-cox | one-parameter Box-Cox transform | + +-------------------+-----------------------------------------------------+ + | db or decibel | transform to units of decibel | + +-------------------+-----------------------------------------------------+ + | log | log transform | + +-------------------+-----------------------------------------------------+ + | nqt | Normal Quantile Transform | + +-------------------+-----------------------------------------------------+ + | sqrt | square-root transform | + +-------------------+-----------------------------------------------------+ + """ if name is None: @@ -107,30 +140,50 @@ def get_method(name, **kwargs): def donothing(R, metadata=None, *args, **kwargs): return R.copy(), {} if metadata is None else metadata.copy() - methods_objects = dict() - methods_objects["none"] = donothing + methods_objects = dict() + methods_objects["none"] = donothing + # arrays methods methods_objects["centred_coord"] = arrays.compute_centred_coord_array + + # cleansing methods + methods_objects["decluster"] = cleansing.decluster + methods_objects["detect_outliers"] = cleansing.detect_outliers + # conversion methods - methods_objects["mm/h"] = conversion.to_rainrate - methods_objects["rainrate"] = conversion.to_rainrate - methods_objects["mm"] = conversion.to_raindepth - methods_objects["raindepth"] = conversion.to_raindepth - methods_objects["dbz"] = conversion.to_reflectivity - methods_objects["reflectivity"] = conversion.to_reflectivity - # transformation methods - methods_objects["boxcox"] = transformation.boxcox_transform - methods_objects["box-cox"] = transformation.boxcox_transform - methods_objects["db"] = transformation.dB_transform - methods_objects["decibel"] = transformation.dB_transform - methods_objects["log"] = transformation.boxcox_transform - methods_objects["nqt"] = transformation.NQ_transform - methods_objects["sqrt"] = transformation.sqrt_transform + methods_objects["mm/h"] = conversion.to_rainrate + methods_objects["rainrate"] = conversion.to_rainrate + methods_objects["mm"] = conversion.to_raindepth + methods_objects["raindepth"] = conversion.to_raindepth + methods_objects["dbz"] = conversion.to_reflectivity + methods_objects["reflectivity"] = conversion.to_reflectivity + # dimension methods - methods_objects["accumulate"] = dimension.aggregate_fields_time - methods_objects["clip"] = dimension.clip_domain - methods_objects["square"] = dimension.square_domain - methods_objects["upscale"] = dimension.aggregate_fields_space + methods_objects["accumulate"] = dimension.aggregate_fields_time + methods_objects["clip"] = dimension.clip_domain + methods_objects["square"] = dimension.square_domain + methods_objects["upscale"] = dimension.aggregate_fields_space + + # image processing methods + methods_objects["shitomasi"] = images.ShiTomasi_detection + methods_objects["morph_opening"] = images.morph_opening + + # interpolation methods + methods_objects["rbfinterp2d"] = interpolate.rbfinterp2d + + # spectral methods + methods_objects["rapsd"] = spectral.rapsd + methods_objects["rm_rdisc"] = spectral.remove_rain_norain_discontinuity + + # transformation methods + methods_objects["boxcox"] = transformation.boxcox_transform + methods_objects["box-cox"] = transformation.boxcox_transform + methods_objects["db"] = transformation.dB_transform + methods_objects["decibel"] = transformation.dB_transform + methods_objects["log"] = transformation.boxcox_transform + methods_objects["nqt"] = transformation.NQ_transform + methods_objects["sqrt"] = transformation.sqrt_transform + # FFT methods if name in ["numpy", "pyfftw", "scipy"]: if "shape" not in kwargs.keys(): @@ -140,11 +193,11 @@ def donothing(R, metadata=None, *args, **kwargs): try: return methods_objects[name] except KeyError as e: - raise ValueError("Unknown method %s\n" % e + - "Supported methods:%s" % str(methods_objects.keys())) - # spectral methods - methods_objects["rapsd"] = spectral.rapsd - methods_objects["rm_rdisc"] = spectral.remove_rain_norain_discontinuity + raise ValueError( + "Unknown method %s\n" % e + + "Supported methods:%s" % str(methods_objects.keys()) + ) + def _get_fft_method(name, **kwargs): kwargs = kwargs.copy() @@ -158,6 +211,8 @@ def _get_fft_method(name, **kwargs): elif name == "pyfftw": return fft.get_pyfftw(shape, **kwargs) else: - raise ValueError("Unknown method {}\n".format(name) - + "The available methods are:" - + str(["numpy", "pyfftw", "scipy"])) from None + raise ValueError( + "Unknown method {}\n".format(name) + + "The available methods are:" + + str(["numpy", "pyfftw", "scipy"]) + ) from None diff --git a/pysteps/utils/interpolate.py b/pysteps/utils/interpolate.py new file mode 100644 index 000000000..c137952ea --- /dev/null +++ b/pysteps/utils/interpolate.py @@ -0,0 +1,223 @@ +# -*- coding: utf-8 -*- +""" +pysteps.utils.interpolate +========================= + +Interpolation routines for pysteps. + +.. autosummary:: + :toctree: ../generated/ + + rbfinterp2d + +""" + +import numpy as np +import scipy.spatial + + +def rbfinterp2d( + coord, + input_array, + xgrid, + ygrid, + rbfunction="gaussian", + epsilon=5, + k=50, + nchunks=5, +): + """Fast 2-D grid interpolation of a sparse (multivariate) array using a + radial basis function. + + Parameters + ---------- + + coord : array_like + Array of shape (n, 2) containing the coordinates of the data points into + a 2-dimensional space. + + input_array : array_like + Array of shape (n) or (n, m) containing the values of the data points, + where *n* is the number of data points and *m* the number of co-located + variables. + All values in **input_array** are required to have finite values. + + xgrid, ygrid : array_like + 1D arrays representing the coordinates of the 2-D output grid. + + rbfunction : {"gaussian", "multiquadric", "inverse quadratic", "inverse multiquadric", "bump"}, optional + The name of one of the available radial basis function based on a + normalized Euclidian norm. + + See also the Notes section below. + + epsilon : float, optional + The shape parameter used to scale the input to the radial kernel. + + A smaller value for **epsilon** produces a smoother interpolation. More + details provided in the wikipedia reference page. + + k : int or None, optional + The number of nearest neighbours used to speed-up the interpolation. + If set to None, it interpolates based on all the data points. + + nchunks : int, optional + The number of chunks in which the grid points are split to limit the + memory usage during the interpolation. + + Returns + ------- + + output_array : array_like + The interpolated field(s) having shape (m, ygrid.size, xgrid.size). + + Notes + ----- + + The coordinates are normalized before computing the Euclidean norms: + + x = (x - min(x)) / max[max(x) - min(x), max(y) - min(y)], + y = (y - min(y)) / max[max(x) - min(x), max(y) - min(y)], + + where the min and max values are taken as the 2nd and 98th percentiles. + + References + ---------- + + Wikipedia contributors, "Radial basis function," Wikipedia, The Free Encyclopedia, + https://en.wikipedia.org/w/index.php?title=Radial_basis_function&oldid=906155047 + (accessed August 19, 2019). + """ + + _rbfunctions = [ + "nearest", + "gaussian", + "inverse quadratic", + "inverse multiquadric", + "bump", + ] + + input_array = np.copy(input_array) + + if np.any(~np.isfinite(input_array)): + raise ValueError("input_array contains non-finite values") + + if input_array.ndim == 1: + nvar = 1 + input_array = input_array[:, None] + + elif input_array.ndim == 2: + nvar = input_array.shape[1] + + else: + raise ValueError( + "input_array must have 1 (n) or 2 dimensions (n, m), but it has %i" + % input_array.ndim + ) + + npoints = input_array.shape[0] + + coord = np.copy(coord) + + if coord.ndim != 2: + raise ValueError( + "coord must have 2 dimensions (n, 2), but it has %i" % coord.ndim + ) + + if npoints != coord.shape[0]: + raise ValueError( + "the number of samples in the input_array does not match the " + + "number of coordinates %i!=%i" % (npoints, coord.shape[0]) + ) + + # normalize coordinates + qcoord = np.percentile(coord, [2, 98], axis=0) + dextent = np.max(np.diff(qcoord, axis=0)) + coord = ( coord - qcoord[0, :] ) / dextent + + rbfunction = rbfunction.lower() + if rbfunction not in _rbfunctions: + raise ValueError( + "Unknown rbfunction '{}'\n".format(rbfunction) + + "The available rbfunctions are: " + + str(_rbfunctions) + ) from None + + # generate the target grid + X, Y = np.meshgrid(xgrid, ygrid) + grid = np.column_stack((X.ravel(), Y.ravel())) + # normalize the grid coordinates + grid = (grid - qcoord[0, :] ) / dextent + + # k-nearest interpolation + if k is not None and k > 0: + k = int(np.min((k, npoints))) + + # create cKDTree object to represent source grid + tree = scipy.spatial.cKDTree(coord) + + else: + k = 0 + + # split grid points in n chunks + if nchunks > 1: + subgrids = np.array_split(grid, nchunks, 0) + subgrids = [x for x in subgrids if x.size > 0] + + else: + subgrids = [grid] + + # loop subgrids + i0 = 0 + output_array = np.zeros((grid.shape[0], nvar)) + for i, subgrid in enumerate(subgrids): + idelta = subgrid.shape[0] + + if k == 0: + # use all points + d = scipy.spatial.distance.cdist( + coord, subgrid, "euclidean" + ).transpose() + inds = np.arange(npoints)[None, :] * np.ones( + (subgrid.shape[0], npoints) + ).astype(int) + + else: + # use k-nearest neighbours + d, inds = tree.query(subgrid, k=k) + + if k == 1: + # nearest neighbour + output_array[i0 : (i0 + idelta), :] = input_array[inds, :] + + else: + + # the interpolation weights + if rbfunction == "gaussian": + w = np.exp(-(d * epsilon) ** 2) + + elif rbfunction == "inverse quadratic": + w = 1.0 / (1 + (epsilon * d) ** 2) + + elif rbfunction == "inverse multiquadric": + w = 1.0 / np.sqrt(1 + (epsilon * d) ** 2) + + elif rbfunction == "bump": + w = np.exp(-1.0 / (1 - (epsilon * d) ** 2)) + w[d >= 1 / epsilon] = 0.0 + + if not np.all(np.sum(w, axis=1)): + w[np.sum(w, axis=1) == 0, :] = 1.0 + + # interpolate + for j in range(nvar): + output_array[i0 : (i0 + idelta), j] = np.sum( + w * input_array[inds, j], axis=1 + ) / np.sum(w, axis=1) + + i0 += idelta + + # reshape to final grid size + output_array = output_array.reshape(ygrid.size, xgrid.size, nvar) + + return np.moveaxis(output_array, -1, 0).squeeze()