Add Farneback motion estimation method

examples/optical_flow_methods_convergence.py

@@ -374,4 +374,22 @@ def plot_optflow_method_convergence(input_precip, optflow_method_name, motion_ty
 # ~~~~~~~~~~~~~~~~~
 plot_optflow_method_convergence(reference_field, "DARTS", "rotor")
 
+################################################################################
+# Farneback
+# ---------
+#
+# Constant motion x-direction
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+plot_optflow_method_convergence(reference_field, "farneback", "linear_x")
+
+################################################################################
+# Constant motion y-direction
+# ~~~~~~~~~~~~~~~~~~~~~~~~~~~
+plot_optflow_method_convergence(reference_field, "farneback", "linear_y")
+
+################################################################################
+# Rotational motion
+# ~~~~~~~~~~~~~~~~~
+plot_optflow_method_convergence(reference_field, "farneback", "rotor")
+
 # sphinx_gallery_thumbnail_number = 5

examples/plot_optical_flow.py

@@ -130,7 +130,7 @@
 #
 # This module implements the anisotropic diffusion method presented in Proesmans
 # et al. (1994), a robust optical flow technique which employs the notion of
-# inconsitency during the solution of the optical flow equations.
+# inconsistency during the solution of the optical flow equations.
 
 oflow_method = motion.get_method("proesmans")
 R[~np.isfinite(R)] = metadata["zerovalue"]
@@ -141,4 +141,22 @@
 quiver(V4, geodata=metadata, step=25)
 plt.show()
 
+################################################################################
+# Farnebäck smoothed method
+# -------------------------
+#
+# This module implements the pyramidal decomposition method for motion estimation
+# of Farnebäck as implemented in OpenCV, with an option for smoothing and
+# renormalization of the motion fields proposed by Driedger et al.:
+# https://cmosarchives.ca/Congress_P_A/program_abstracts2022.pdf (p. 392).
+
+oflow_method = motion.get_method("farneback")
+R[~np.isfinite(R)] = metadata["zerovalue"]
+V5 = oflow_method(R[-2:, :, :], verbose=True)
+
+# Plot the motion field
+plot_precip_field(R_, geodata=metadata, title="Farneback")
+quiver(V5, geodata=metadata, step=25)
+plt.show()
+
 # sphinx_gallery_thumbnail_number = 1

pysteps/motion/farneback.py

@@ -0,0 +1,269 @@
+# -*- coding: utf-8 -*-
+"""
+pysteps.motion.farneback
+========================
+
+The Farneback dense optical flow module.
+
+This module implements the interface to the local `Farneback`_ routine
+available in OpenCV_.
+
+.. _OpenCV: https://opencv.org/
+
+.. _`Farneback`:\
+    https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga5d10ebbd59fe09c5f650289ec0ece5af
+
+.. autosummary::
+    :toctree: ../generated/
+
+    farneback
+"""
+
+import numpy as np
+from numpy.ma.core import MaskedArray
+import scipy.ndimage as sndi
+import time
+
+from pysteps.decorators import check_input_frames
+from pysteps.exceptions import MissingOptionalDependency
+from pysteps.utils.images import morph_opening
+
+try:
+    import cv2
+
+    CV2_IMPORTED = True
+except ImportError:
+    CV2_IMPORTED = False
+
+
+@check_input_frames(2)
+def farneback(
+    input_images,
+    pyr_scale=0.5,
+    levels=3,
+    winsize=15,
+    iterations=3,
+    poly_n=5,
+    poly_sigma=1.1,
+    flags=0,
+    size_opening=3,
+    sigma=60.0,
+    verbose=False,
+):
+    """Estimate a dense motion field from a sequence of 2D images using
+    the `Farneback`_ optical flow algorithm.
+
+    This function computes dense optical flow between each pair of consecutive
+    input frames using OpenCV's Farneback method. If more than two frames are
+    provided, the motion fields estimated from all consecutive pairs are
+    averaged to obtain a single representative advection field.
+
+    After the pairwise motion fields are averaged, the resulting motion field
+    can optionally be smoothed with a Gaussian filter. In that case, its
+    amplitude is rescaled so that the mean motion magnitude is preserved.
+
+    .. _OpenCV: https://opencv.org/
+
+    .. _`Farneback`:\
+        https://docs.opencv.org/3.4/dc/d6b/group__video__track.html#ga5d10ebbd59fe09c5f650289ec0ece5af
+
+    .. _MaskedArray:\
+        https://docs.scipy.org/doc/numpy/reference/maskedarray.baseclass.html#numpy.ma.MaskedArray
+
+    .. _ndarray:\
+    https://docs.scipy.org/doc/numpy/reference/generated/numpy.ndarray.html
+
+    Parameters
+    ----------
+    input_images: ndarray_ 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 motion vectors are averaged together.
+
+        In case of ndarray_, 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.
+
+    pyr_scale : float, optional
+        Parameter specifying the image scale (<1) used to build pyramids for
+        each image; pyr_scale=0.5 means a classical pyramid, where each next
+        layer is twice smaller than the previous one. This and the following
+        parameter descriptions are adapted from the original OpenCV
+        documentation (see https://docs.opencv.org).
+
+    levels : int, optional
+        Number of pyramid layers including the initial image; levels=1 means
+        that no extra layers are created and only the original images are used.
+
+    winsize : int, optional
+        Averaging window size; larger values increase the algorithm robustness
+        to image noise and give more stable motion estimates. Small windows
+        (e.g. 10) lead to unrealistic motion.
+    iterations : int, optional
+        Number of iterations the algorithm does at each pyramid level.
+
+    poly_n : int
+        Size of the pixel neighborhood used to find polynomial expansion in
+        each pixel; larger values mean that the image will be approximated with
+        smoother surfaces, yielding more robust algorithm and more blurred
+        motion field, typically poly_n = 5 or 7.
+
+    poly_sigma : float
+        Standard deviation of the Gaussian that is used to smooth derivatives
+        used as a basis for the polynomial expansion; for poly_n=5, you can set
+        poly_sigma=1.1, for poly_n=7, a good value would be poly_sigma=1.5.
+
+    flags : int, optional
+        Operation flags that can be a combination of the following:
+
+        OPTFLOW_USE_INITIAL_FLOW uses the input 'flow' as an initial flow
+        approximation.
+
+        OPTFLOW_FARNEBACK_GAUSSIAN uses the Gaussian winsize x winsize filter
+        instead of a box filter of the same size for optical flow estimation;
+        usually, this option gives a more accurate flow than with a box filter,
+        at the cost of lower speed; normally, winsize for a Gaussian window
+        should be set to a larger value to achieve the same level of robustness.
+
+    size_opening : int, optional
+        Non-OpenCV parameter:
+        The structuring element size for the filtering of isolated pixels [px].
+
+    sigma : float, optional
+        Non-OpenCV parameter:
+        The smoothing bandwidth of the motion field. The motion field amplitude
+        is adjusted by multiplying by the ratio of average magnitude before and
+        after smoothing to avoid damping of the motion field.
+
+    verbose: bool, optional
+        If set to True, print some information about the program.
+
+    Returns
+    -------
+    out : ndarray_, shape (2,m,n)
+        Return the advection field having shape
+        (2, m, n), where out[0, :, :] contains the x-components of the motion
+        vectors and out[1, :, :] contains the y-components.
+        The velocities are in units of pixels / timestep, where timestep is the
+        time difference between the two input images.
+        Return a zero motion field of shape (2, m, n) when no motion is
+        detected.
+
+    References
+    ----------
+    Farnebäck, G.: Two-frame motion estimation based on polynomial expansion, 
+    In Image Analysis, pages 363–370. Springer, 2003.
+    Driedger, N., Mahidjiba, A. and Hortal, A.P. (2022, June 1-8): Evaluation of optical flow
+    methods for radar precipitation extrapolation.
+    Canadian Meteorological and Oceanographic Society Congress, contributed abstract 11801.
+    """
+
+    if len(input_images.shape) != 3:
+        raise ValueError(
+            "input_images has %i dimensions, but a "
+            "three-dimensional array is expected" % len(input_images.shape)
+        )
+
+    input_images = input_images.copy()
+
+    if verbose:
+        print("Computing the motion field with the Farneback method.")
+        t0 = time.time()
+
+    if not CV2_IMPORTED:
+        raise MissingOptionalDependency(
+            "OpenCV (cv2) is required for the Farneback optical flow method, but it is not installed"
+        )
+
+    nr_pairs = input_images.shape[0] - 1
+    domain_size = (input_images.shape[1], input_images.shape[2])
+    u_sum = np.zeros(domain_size)
+    v_sum = np.zeros(domain_size)
+    for n in range(nr_pairs):
+        # extract consecutive images
+        prvs_img = input_images[n, :, :].copy()
+        next_img = input_images[n + 1, :, :].copy()
+
+        # Check if a MaskedArray is used. If not, mask the ndarray
+        if not isinstance(prvs_img, MaskedArray):
+            prvs_img = np.ma.masked_invalid(prvs_img)
+        np.ma.set_fill_value(prvs_img, prvs_img.min())
+
+        if not 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
+        im_min = prvs_img.min()
+        im_max = prvs_img.max()
+        if (im_max - im_min) > 1e-8:
+            prvs_img = (prvs_img.filled() - im_min) / (im_max - im_min) * 255
+        else:
+            prvs_img = prvs_img.filled() - im_min
+
+        im_min = next_img.min()
+        im_max = next_img.max()
+        if (im_max - im_min) > 1e-8:
+            next_img = (next_img.filled() - im_min) / (im_max - im_min) * 255
+        else:
+            next_img = next_img.filled() - im_min
+
+        # convert to 8-bit
+        prvs_img = np.ndarray.astype(prvs_img, "uint8")
+        next_img = np.ndarray.astype(next_img, "uint8")
+
+        # remove small noise with a morphological operator (opening)
+        if size_opening > 0:
+            prvs_img = morph_opening(prvs_img, prvs_img.min(), size_opening)
+            next_img = morph_opening(next_img, next_img.min(), size_opening)
+
+        flow = cv2.calcOpticalFlowFarneback(
+            prvs_img,
+            next_img,
+            None,
+            pyr_scale,
+            levels,
+            winsize,
+            iterations,
+            poly_n,
+            poly_sigma,
+            flags,
+        )
+
+        fa, fb = np.dsplit(flow, 2)
+        u_sum += fa.reshape(domain_size)
+        v_sum += fb.reshape(domain_size)
+
+    # Compute the average motion field
+    u = u_sum / nr_pairs
+    v = v_sum / nr_pairs
+
+    # Smoothing
+    if sigma > 0:
+        uv2 = u * u + v * v  # squared magnitude of motion field
+        us = sndi.gaussian_filter(u, sigma, mode="nearest")
+        vs = sndi.gaussian_filter(v, sigma, mode="nearest")
+        uvs2 = us * us + vs * vs  # squared magnitude of smoothed motion field
+
+        mean_uv2 = np.nanmean(uv2)
+        mean_uvs2 = np.nanmean(uvs2)
+        if mean_uvs2 > 0:
+            mult = np.sqrt(mean_uv2 / mean_uvs2)
+        else:
+            mult = 1.0
+    else:
+        mult = 1.0
+        us = u
+        vs = v
+    if verbose:
+        print("mult factor of smoothed motion field=", mult)
+
+    UV = np.stack([us * mult, vs * mult])
+
+    if verbose:
+        print("--- %s seconds ---" % (time.time() - t0))
+
+    return UV

pysteps/motion/interface.py

@@ -31,6 +31,7 @@
 from pysteps.motion.lucaskanade import dense_lucaskanade
 from pysteps.motion.proesmans import proesmans
 from pysteps.motion.vet import vet
+from pysteps.motion.farneback import farneback
 
 _methods = dict()
 _methods["constant"] = constant
@@ -39,6 +40,7 @@
 _methods["darts"] = DARTS
 _methods["proesmans"] = proesmans
 _methods["vet"] = vet
+_methods["farneback"] = farneback
 _methods[None] = lambda precip, *args, **kw: np.zeros(
     (2, precip.shape[1], precip.shape[2])
 )
@@ -73,6 +75,8 @@ def get_method(name):
     |                   | Laroche and Zawadzki (1995) and                      |
     |                   | Germann and Zawadzki (2002)                          |
     +-------------------+------------------------------------------------------+
+    |  farneback        | OpenCV implementation of the Farneback (2003) method.|
+    +-------------------+------------------------------------------------------+
 
     +--------------------------------------------------------------------------+
     | Methods implemented in C (these require separate compilation and linkage)|