Merge pull request OpenPIV#327 from nepomnyi/master

Added normalized median filter to validation.py (i.e. the 2005 version of Westerweel's median filter).

Author: alexlib

openpiv/test/OpenPIV_results_16_test1/field_A0000.png

@@ -76,6 +76,62 @@ def test_local_median_validation(u_threshold=3, N=3):
     assert flag[N, N]
 
 
+def test_local_norm_median_validation():
+    """ test normalized local median
+
+    Args:
+        threshold (int, optional): _description_.
+        N (int, optional): _description_.
+    """
+
+    # The idea is the following. Let's create an array of the size of one filtering kernel.
+    # Calculate normalized median filter by hand and using local_norm_median_val(). See if
+    # the results are comparable. As outlined in the description of local_norm_median_val(),
+    # follow the paper J. Westerweel, F. Scarano, "Universal outlier detection for PIV data",
+    # Experiments in fluids, 39(6), p.1096-1100, 2005, equation 2, to perform calculations
+    # by hand.
+
+    # Case study:
+    u = np.asarray([[1,2,3],[4,5,6],[7,8,9]])
+    v = np.asarray([[2,3,4],[5,6,7],[8,9,10]])
+    ε = 0.1
+    thresh = 2
+
+    # Calculations by hand (according to equation 2 in the referenced article).
+    u0 = 5
+    v0 = 6
+    # Median formula for even number of observations (the center is excluded)
+    # that are placed in the ascending order (just like the given u and v arrays):
+    # Median = [(n/2)th term + ((n/2) + 1)th term]/2 (https://www.cuemath.com/data/median/).
+    um = 5 # = (4 + 6) / 2
+    vm = 6 # = (5 + 7) / 2 
+    # Arrays of residuals given by the formula rui = |ui-um| and rvi = |vi-vm|
+    rui = np.asarray([[4,3,2],[1,0,1],[2,3,4]])
+    rvi = np.asarray([[4,3,2],[1,0,1],[2,3,4]])
+    # Median residuals values (center excluded) must be calculated for ascending arrays:
+    ruiascend = [1, 1, 2, 2, 3, 3, 4, 4]
+    rviascend = [1, 1, 2, 2, 3, 3, 4, 4]
+    rum = 2.5 # = (2 + 3) / 2
+    rvm = 2.5 # = (2 + 3) / 2
+    # Now implement equation 2 from the referenced article:
+    r0ast_u = np.abs(u0 - um) / (rum + ε)
+    r0ast_v = np.abs(v0 - vm) / (rvm + ε)
+    # The method of comparison to the threshold is given at the very end of the referenced 
+    # article in the MATLAB code:
+    byHand = (r0ast_u**2 + r0ast_v**2)**0.5 > thresh
+    # Here, we calculated byHand for the velocity vector (u0,v0) coordinates of which are 
+    # located at u[1,1] and v[1,1].
+
+    # Calculations using local_norm_median_val() function.
+    byFunc = validation.local_norm_median_val(u, v, ε, thresh)
+    # Here by func returns a boolean matrix containing either True or False for each velocity
+    # vector (ui,vi).
+
+    # Compare the two results:
+    assert byHand == byFunc[1,1]
+
+
+
 def test_global_val():
     """ tests global validation
     """

openpiv/test/test_validation.py

@@ -202,7 +202,7 @@ def local_median_val(
     Returns
     -------
 
-    flag : boolean 2d np.ndarray
+    ind : boolean 2d np.ndarray
         a boolean array. True elements corresponds to outliers.
 
     """
@@ -228,6 +228,126 @@ def local_median_val(
     return ind
 
 
+def local_norm_median_val(
+        u: np.ndarray, 
+        v: np.ndarray, 
+        ε: float,
+        threshold: float,
+        size: int=1
+        )->np.ndarray:
+    """This function is adapted from OpenPIV's implementation of 
+    validation.local_median_val(). validation.local_median_val() is, 
+    basically, Westerweel's original median filter (with some changes). 
+    The current function builts upon validation.local_median_val() and implements
+    improved Westerweel's median filter (normalized filter) as described 
+    in 2007 edition of the German PIV book (paragraph 6.1.5) and in Westerweel's
+    article J. Westerweel, F. Scarano, "Universal outlier detection for PIV data",
+    Experiments in fluids, 39(6), p.1096-1100, 2005.
+    For the list of parameters, see the referenced article, equation 2 on p.1097.
+    The current function implements equation 2 from the referenced article in a
+    manner shown in the MATLAB script at the end of the article, on p.1100.
+
+    This validation method tests for the spatial consistency of the data.
+    Vectors are classified as outliers and replaced with Nan (Not a Number) if
+    the absolute difference with the local median is greater than a user
+    specified threshold. The median is computed for both velocity components.
+
+    The image masked areas (obstacles, reflections) are marked as masked array:
+       u = np.ma.masked(u, flag = image_mask)
+    and it should not be replaced by the local median, but remain masked. 
+
+
+    Parameters
+    ----------
+    u : 2d np.ndarray
+        a two dimensional array containing the u velocity component
+
+    v : 2d np.ndarray
+        a two dimensional array containing the v velocity component
+
+    ε : float
+        minimum normalization level (see the referenced article, eqn.2)
+
+    threshold : float
+        the threshold to determine whether the vector is valid or not
+
+    size: int
+        the representative size of the kernel of the median filter, the 
+        actual size of the kernel is (2*size+1, 2*size+1) - i.e., it's the
+        number of interrogation windows away from the interrogation 
+        window of interest
+
+    Returns
+    -------
+
+    ind : boolean 2d np.ndarray
+        a boolean array; true elements corresponds to outliers
+
+    """
+    if np.ma.is_masked(u):
+        masked_u = np.where(~u.mask, u.data, np.nan)
+        masked_v = np.where(~v.mask, v.data, np.nan)
+    else:
+        masked_u = u
+        masked_v = v
+
+    um = generic_filter(masked_u, 
+                        np.nanmedian, 
+                        mode='constant',
+                        cval=np.nan, 
+                        size=(2*size+1, 2*size+1)
+    )
+    vm = generic_filter(masked_v, 
+                        np.nanmedian, 
+                        mode='constant',
+                        cval=np.nan, 
+                        size=(2*size+1, 2*size+1)
+    )
+
+    def rfunc(x):
+        """
+        Implementation of r from the cited article (see the description of
+        the function above). x is the array within the filtering kernel. 
+        I.e., every element of x is a velocity vector ui or vi.
+        This function must return a scalar: https://stackoverflow.com/a/14060024/10073233
+        """
+        # copied from here: https://stackoverflow.com/a/60166608/10073233
+        y = x.copy() # need this step, because np.put() below changes the array in place,
+                     # and we can end up with a situation when the entire filtering kernel
+                     # is comprised of NaNs resulting in NumPy RuntimeWarning: All-NaN slice encountered
+        np.put(y, y.size//2, np.nan) # put NaN in the middle to avoid using
+                                     # the middle in the calculations
+        ym = np.nanmedian(y) # Um for the current filtering window
+        rm = np.nanmedian(np.abs(np.subtract(y,ym))) # median of |ui-um| or |vi-vm|
+        return rm
+
+    rm_u = generic_filter(masked_u, 
+                          rfunc, 
+                          mode='constant',
+                          cval=np.nan, 
+                          size=(2*size+1, 2*size+1)
+    )
+    rm_v = generic_filter(masked_v, 
+                          rfunc, 
+                          mode='constant',
+                          cval=np.nan, 
+                          size=(2*size+1, 2*size+1)
+    )
+
+    r0ast_u = np.divide(np.abs(np.subtract(masked_u,um)), np.add(rm_u,ε)) # r0ast stands for r_0^* - 
+                                                                          # see formula 2 in the 
+                                                                          # referenced article 
+                                                                          # (see description of the function)
+    r0ast_v = np.divide(np.abs(np.subtract(masked_v,vm)), np.add(rm_v,ε)) # r0ast stands for r_0^* - 
+                                                                          # see formula 2 in the 
+                                                                          # referenced article 
+                                                                          # (see description of the function)
+
+    ind = (np.sqrt(np.add(np.square(r0ast_u),np.square(r0ast_v)))) > threshold
+
+    return ind
+
+
 def typical_validation(
     u: np.ndarray,
     v: np.ndarray,