def METHOD_NAME(data, filter_type, filter_par=None):
    """Filter 2D matrix with selected filter.
    Parameters: data        - 2D np.array, matrix to be filtered
                filter_type - string, filter type
                filter_par  - (list of) int/float, optional, parameter for low/high pass filter
                              for low/highpass_avg, it's kernel size in int
                              for low/highpass_gaussain, it's sigma in float
                              for double_difference, it's local and regional kernel sizes in int
    Returns:    data_filt   - 2D np.array, matrix after filtering.
    """
    if filter_type == "sobel":
        data_filt = filters.sobel(data)
    elif filter_type == "roberts":
        data_filt = filters.roberts(data)
    elif filter_type == "canny":
        data_filt = feature.canny(data)
    elif filter_type == "lowpass_avg":
        p = int(filter_par)
        kernel = np.ones((p, p), np.float32)/(p*p)
        data_filt = ndimage.convolve(data, kernel)
    elif filter_type == "highpass_avg":
        p = int(filter_par)
        kernel = np.ones((p, p), np.float32)/(p*p)
        lp_data = ndimage.convolve(data, kernel)
        data_filt = data - lp_data
    elif filter_type == "lowpass_gaussian":
        # ORIGINAL: data_filt = filters.gaussian(data, sigma=filter_par)
        #   nan pixels can enlarge to big holes depending on the size of your gaussian kernel
        #   we can do normalized convolution (https://stackoverflow.com/a/36307291/7128154) as below:
        V=np.array(data)
        V[np.isnan(data)]=0
        VV=filters.gaussian(V, sigma=filter_par)
        W=np.ones_like(data)
        W[np.isnan(data)]=0
        WW=filters.gaussian(W, sigma=filter_par)
        WW[WW==0]=np.nan
        data_filt = VV/WW
        data_filt[np.isnan(data)] = np.nan
    elif filter_type == "highpass_gaussian":
        # Use the existing logic for lowpass_gaussian, then subtract from original data
        return data - METHOD_NAME(data, "lowpass_gaussian", filter_par=filter_par)
    elif filter_type == "median":
        data_filt = filters.median(data, morphology.disk(filter_par))
    elif filter_type == "double_difference":
        """Amplifies the local deformation signal by reducing the influence
        of regional deformation trends from atmospheric artifacts, tectonic
        deformation, and other sources. Intended use is to identify landslide-related
        deformation. Filter has the form:
            result =  regional mean of data - local mean of data
        where both the regional and local kernel size can be set using the
        filter_par argument. The regional kernel has a doughnut shape
        because the pixels used to calculate the local mean are not
        included in the regional mean. This ensures that the local mean
        and regional mean results are separate.
        """
        local_kernel = morphology.disk(filter_par[0], dtype=np.float32)
        local_kernel = np.pad(local_kernel, pad_width=int(filter_par[1]-filter_par[0]), mode='constant')
        regional_kernel = morphology.disk(filter_par[1], dtype=np.float32)
        regional_kernel[local_kernel == 1] = 0
        # normalize
        local_kernel /= local_kernel.sum(axis=(0,1))
        regional_kernel /= regional_kernel.sum(axis=(0,1))
        # double-difference kernel
        combined_kernel = regional_kernel - local_kernel
        data_filt = ndimage.convolve(data, combined_kernel)
    else:
        raise Exception('Un-recognized filter type: '+filter_type)
    return data_filt