Python scipy.ndimage.filters.maximum_filter() Examples

def get_2D_peaks(arr2D, amp_min=DEFAULT_AMP_MIN):
    struct = generate_binary_structure(2, 1)
    neighborhood = iterate_structure(struct, PEAK_NEIGHBORHOOD_SIZE)

    local_max = maximum_filter(arr2D, footprint=neighborhood) == arr2D
    background = (arr2D == 0)
    eroded_background = binary_erosion(background, structure=neighborhood,
                                       border_value=1)

    detected_peaks = local_max ^ eroded_background
    amps = arr2D[detected_peaks]
    j, i = np.where(detected_peaks)
    amps = amps.flatten()
    
    peaks = zip(i, j, amps)
    peaks_filtered = [x for x in peaks if x[2] > amp_min] 

    frequency_idx = [x[1] for x in peaks_filtered]
    time_idx = [x[0] for x in peaks_filtered]
    
    return zip(frequency_idx, time_idx) 

def findpeaks(image, thresh):
    """
    Return positions of all peaks in image above threshold thresh
    Based on `"detect_peaks" Stack Overflow discussion <https://stackoverflow.com/questions/3684484/peak-detection-in-a-2d-array/3689710#3689710>`_
    
    :param image: array of values to search
    :param thresh: threshold for peaks
    :type image: numpy.ndarray
    :type thresh: float
    :returns: index array (equivalent of where output)
    """
    # define an 8-connected neighborhood
    neighborhood = generate_binary_structure(2,2)
    # find local maximum for each pixel
    amax = maximum_filter(image, footprint=neighborhood)
    w = numpy.where((image == amax) & (image >= thresh))
    return w 

def filter_local_maxima_with_std(self, ft2d):
        data = np.abs(np.fft.fftshift(ft2d))
        data /= (np.max(data) + 1e-7)

        data_max = maximum_filter(data, self.neighborhood_size)
        data_mean = uniform_filter(data, self.neighborhood_size)
        data_mean_sq = uniform_filter(data ** 2, self.neighborhood_size)
        data_std = np.sqrt(data_mean_sq - data_mean ** 2) + 1e-7

        maxima = ((data_max - data_mean) / data_std)
        fraction_of_local_max = (data == data_max)
        maxima *= fraction_of_local_max
        maxima = maxima.astype(float)
        maxima /= (np.max(maxima) + 1e-7)

        maxima = np.maximum(maxima, np.fliplr(maxima), np.flipud(maxima))
        maxima = np.fft.ifftshift(maxima)

        background_ft2d = np.multiply(maxima, ft2d)
        foreground_ft2d = np.multiply(1 - maxima, ft2d)

        return background_ft2d, foreground_ft2d 

def localMax(img, size = 5):
    """Calculates local maxima of an image
        
    Arguments:
        img (array): image
        size (float or None): size of volume to search for maxima
        
    Returns:
        array: mask that is True at local maxima
    """
    
    if size is None:
        return img;
    
    if not isinstance(size, tuple):
       size = (size, size, size);

    #return regmin(-img, regionalMaxStructureElement);
    #return (maximum_filter(img, footprint = regionalMaxStructureElement) == img);
    return (maximum_filter(img, size = size) == img)
    
      
#def regionalMax(img, regionalMaxStructureElement = numpy.ones((3,3,3), dtype = bool)):
#    """Calculates regional maxima of an image."""
#   
#    return regmin(-img, regionalMaxStructureElement); 

def compute_line_seeds(self, binary, bottom, top, colseps, scale):
        """Base on gradient maps, computes candidates for baselines
        and xheights.  Then, it marks the regions between the two
        as a line seed."""
        t = self.parameter['threshold'] 
        vrange = int(self.parameter['vscale']*scale)
        bmarked = maximum_filter(bottom == maximum_filter(bottom, (vrange, 0)), (2, 2))
        bmarked *= np.array((bottom > t*np.amax(bottom)*t)*(1-colseps), dtype=bool)
        tmarked = maximum_filter(top == maximum_filter(top, (vrange, 0)), (2, 2))
        tmarked *= np.array((top > t*np.amax(top)*t/2)*(1-colseps), dtype=bool)
        tmarked = maximum_filter(tmarked, (1, 20))
        testseeds = np.zeros(binary.shape, 'i')
        seeds = np.zeros(binary.shape, 'i')
        delta = max(3, int(scale/2))
        for x in range(bmarked.shape[1]):
            transitions = sorted([(y, 1) for y in psegutils.find(bmarked[:, x])]+[(y, 0) for y in psegutils.find(tmarked[:, x])])[::-1]
            transitions += [(0, 0)]
            for l in range(len(transitions)-1):
                y0, s0 = transitions[l]
                if s0 == 0:
                    continue
                seeds[y0-delta:y0, x] = 1
                y1, s1 = transitions[l+1]
                if s1 == 0 and (y0-y1) < 5*scale:
                    seeds[y1:y0, x] = 1
        seeds = maximum_filter(seeds, (1, int(1+scale)))
        seeds *= (1-colseps)
        seeds, _ = morph.label(seeds)
        return seeds

    
    ################################################################
    # The complete line segmentation process.
    ################################################################ 

def detect_objects_heatmap(heatmap):
    data = 256 * heatmap
    data_max = filters.maximum_filter(data, 3)
    maxima = (data == data_max)
    data_min = filters.minimum_filter(data, 3)
    diff = ((data_max - data_min) > 0.3)
    maxima[diff == 0] = 0
    labeled, num_objects = ndimage.label(maxima)
    slices = ndimage.find_objects(labeled)
    objects = np.zeros((num_objects, 2), dtype=np.int32)
    pidx = 0
    for (dy, dx) in slices:
        pos = [(dy.start + dy.stop - 1) // 2, (dx.start + dx.stop - 1) // 2]
        if heatmap[pos[0], pos[1]] > config.CENTER_TR:
            objects[pidx, :] = pos
            pidx += 1
    return objects[:pidx] 

def find_peaks(param, img):
    """
    Given a (grayscale) image, find local maxima whose value is above a given
    threshold (param['thre1'])
    :param img: Input image (2d array) where we want to find peaks
    :return: 2d np.array containing the [x,y] coordinates of each peak found
    in the image
    """

    peaks_binary = (maximum_filter(img, footprint=generate_binary_structure(
        2, 1)) == img) * (img > param['thre1'])
    # Note reverse ([::-1]): we return [[x y], [x y]...] instead of [[y x], [y
    # x]...]
    return np.array(np.nonzero(peaks_binary)[::-1]).T 

def localMax(V, footprint, thres):
    ''' find local maxima of V (correlation map) using a filter with (usually circular) footprint
        inputs:
            V, footprint, thres
        outputs:
            i,j: indices of local max greater than thres
    '''
    maxV = filters.maximum_filter(V, footprint=footprint, mode = 'reflect')
    imax = V > np.maximum(thres, maxV - 1e-10)
    i,j  = imax.nonzero()
    i    = i.astype(np.int32)
    j    = j.astype(np.int32)
    return i,j 

def local_maxima(arr):
    # http://stackoverflow.com/questions/3684484/peak-detection-in-a-2d-array/3689710#3689710
    """
    Takes an array and detects the troughs using the local maximum filter.
    Returns a boolean mask of the troughs (i.e. 1 when
    the pixel's value is the neighborhood maximum, 0 otherwise)
    """
    # define an connected neighborhood
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.morphology.html#generate_binary_structure
    neighborhood = morphology.generate_binary_structure(len(arr.shape),2)
    # apply the local maximum filter; all locations of maximal value 
    # in their neighborhood are set to 1
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.filters.html#maximum_filter
    local_max = (filters.maximum_filter(arr, footprint=neighborhood)==arr)
    # local_max is a mask that contains the peaks we are 
    # looking for, but also the background.
    # In order to isolate the peaks we must remove the background from the mask.
    # 
    # we create the mask of the background
    background = (arr==arr.min())           # mxu: in the original version, was         background = (arr==0)
    # 
    # a little technicality: we must erode the background in order to 
    # successfully subtract it from local_max, otherwise a line will 
    # appear along the background border (artifact of the local maximum filter)
    # http://www.scipy.org/doc/api_docs/SciPy.ndimage.morphology.html#binary_erosion
    eroded_background = morphology.binary_erosion(
        background, structure=neighborhood, border_value=1)
    # 
    # we obtain the final mask, containing only peaks, 
    # by removing the background from the local_max mask
    #detected_maxima = local_max - eroded_backround             # mxu: this is the old version, but the boolean minus operator is deprecated
    detected_maxima = np.bitwise_and(local_max, np.bitwise_not(eroded_background))          # Material nonimplication, see http://en.wikipedia.org/wiki/Material_nonimplication
    return np.where(detected_maxima) 

def r_dilation(image, size, origin=0):
    """Dilation with rectangular structuring element using maximum_filter"""
    return filters.maximum_filter(image, size, origin=origin) 

def find_peaks(param, img):
    """
    Given a (grayscale) image, find local maxima whose value is above a given
    threshold (param['thre1'])
    :param img: Input image (2d array) where we want to find peaks
    :return: 2d np.array containing the [x,y] coordinates of each peak found
    in the image
    """

    peaks_binary = (maximum_filter(img, footprint=generate_binary_structure(
        2, 1)) == img) * (img > param['thre1'])
    # Note reverse ([::-1]): we return [[x y], [x y]...] instead of [[y x], [y
    # x]...]
    return np.array(np.nonzero(peaks_binary)[::-1]).T 

def find_N_peaks(signal, r=29, min_v=0.05, N=None):
    max_v = maximum_filter(signal, size=r, mode='wrap')
    pk_loc = np.where(max_v == signal)[0]
    pk_loc = pk_loc[signal[pk_loc] > min_v]
    # check for odd case, remove one 
    if (pk_loc.shape[0]%2)!=0:
        pk_id = np.argsort(-signal[pk_loc])
        pk_loc = pk_loc[pk_id[:-1]]
        pk_loc = np.sort(pk_loc)
    if N is not None:
        order = np.argsort(-signal[pk_loc])
        pk_loc = pk_loc[order[:N]]
        pk_loc = pk_loc[np.argsort(pk_loc)]
    return pk_loc, signal[pk_loc] 

def _grow(self, arr):
        return maximum_filter(arr, footprint=self._footprint) 

def filter_local_maxima(self, ft2d):
        data = np.abs(np.fft.fftshift(ft2d))
        data /= (np.max(data) + 1e-7)
        threshold = np.std(data)

        data_max = maximum_filter(data, self.neighborhood_size)
        maxima = (data == data_max)
        data_min = minimum_filter(data, self.neighborhood_size)
        diff = ((data_max - data_min) > threshold)
        maxima[diff == 0] = 0
        maxima = np.maximum(maxima, np.fliplr(maxima), np.flipud(maxima))
        maxima = np.fft.ifftshift(maxima)

        background_ft2d = np.multiply(maxima, ft2d)
        foreground_ft2d = np.multiply(1 - maxima, ft2d)

        return background_ft2d, foreground_ft2d 

def nms(score, ks):
    assert ks % 2 == 1
    ret_score = score.copy()
    maxpool = maximum_filter(score, footprint=np.ones((ks, ks)))
    ret_score[score < maxpool] = 0.
    return ret_score 

def find_peaks(param, img):
    '''
    Given a (grayscale) image, find local maxima whose value is above a given threshold (param['thre1'])
    Args:
        param (dict): 
        img (ndarray): Input image (2d array) where we want to find peaks
    Returns: 
        2d np.array containing the [x,y] coordinates of each peak foun in the image
    '''

    peaks_binary = (maximum_filter(img, footprint=generate_binary_structure(2, 1)) == img) * (img > param['thre1'])
    # Note reverse ([::-1]): we return [[x y], [x y]...] instead of [[y x], [y x]...]
    return np.array(np.nonzero(peaks_binary)[::-1]).T 

def find_N_peaks(signal, r=29, min_v=0.05, N=None):
    max_v = maximum_filter(signal, size=r, mode='wrap')
    pk_loc = np.where(max_v == signal)[0]
    pk_loc = pk_loc[signal[pk_loc] > min_v]
    if N is not None:
        order = np.argsort(-signal[pk_loc])
        pk_loc = pk_loc[order[:N]]
        pk_loc = pk_loc[np.argsort(pk_loc)]
    return pk_loc, signal[pk_loc] 

def test_all():
    print("~"*40, " maximum filter")
    _test_some(max_filter, spf.maximum_filter, cval = -np.inf)
    print("~" * 40, " minimum filter")
    _test_some(min_filter, spf.minimum_filter, cval = np.inf)
    print("~" * 40, " uniform filter")
    _test_some(uniform_filter, spf.uniform_filter, cval = 0.) 

def r_erosion(image, size, origin=0):
    """Erosion with rectangular structuring element using maximum_filter"""
    return filters.minimum_filter(image, size, origin=origin) 

def detect(self, image):
        # define an 8-connected neighborhood
        neighborhood = generate_binary_structure(2, 2)

        # apply the local maximum filter; all pixel of maximal value
        # in their neighborhood are set to 1
        local_max = maximum_filter(image, footprint=neighborhood) == image
        # local_max is a mask that contains the peaks we are
        # looking for, but also the background.
        # In order to isolate the peaks we must remove the background from the mask.

        # we create the mask of the background
        background = (image < self.min_th)

        # a little technicality: we must erode the background in order to
        # successfully subtract it form local_max, otherwise a line will
        # appear along the background border (artifact of the local maximum filter)
        eroded_background = binary_erosion(background, structure=neighborhood, border_value=1)

        # we obtain the final mask, containing only peaks,
        # by removing the background from the local_max mask (xor operation)
        detected_peaks = local_max ^ eroded_background

        detected_peaks[image < self.min_th] = False
        peaks = np.array(np.nonzero(detected_peaks)).T

        if len(peaks) == 0:
            return peaks, np.array([])

        # nms
        if len(peaks) == 1:
            clusters = [0]
        else:
            clusters = fclusterdata(peaks, self.min_dist, criterion="distance")
        peak_groups = {}
        for ind_junc, ind_group in enumerate(clusters):
            if ind_group not in peak_groups.keys():
                peak_groups[ind_group] = []
                peak_groups[ind_group].append(peaks[ind_junc])
        peaks_nms = []
        peaks_score = []
        for peak_group in peak_groups.values():
            values = [image[y, x] for y, x in peak_group]
            ind_max = np.argmax(values)
            peaks_nms.append(peak_group[int(ind_max)])
            peaks_score.append(values[int(ind_max)])

        return np.float32(np.array(peaks_nms)), np.float32(np.array(peaks_score)) 

def get_2D_peaks(arr2D, plot=False, amp_min=DEFAULT_AMP_MIN):
    # http://docs.scipy.org/doc/scipy/reference/generated/scipy.ndimage.morphology.iterate_structure.html#scipy.ndimage.morphology.iterate_structure
    struct = generate_binary_structure(2, 1)
    neighborhood = iterate_structure(struct, PEAK_NEIGHBORHOOD_SIZE)

    # find local maxima using our fliter shape
    local_max = maximum_filter(arr2D, footprint=neighborhood) == arr2D
    background = (arr2D == 0)
    eroded_background = binary_erosion(background, structure=neighborhood,
                                       border_value=1)

    # Boolean mask of arr2D with True at peaks
    detected_peaks = local_max ^ eroded_background

    # extract peaks
    amps = arr2D[detected_peaks]
    j, i = np.where(detected_peaks)

    # filter peaks
    amps = amps.flatten()
    peaks = zip(i, j, amps)
    peaks_filtered = [x for x in peaks if x[2] > amp_min]  # freq, time, amp

    # get indices for frequency and time
    frequency_idx = [x[1] for x in peaks_filtered]
    time_idx = [x[0] for x in peaks_filtered]

    # scatter of the peaks
    if plot:
      fig, ax = plt.subplots()
      ax.imshow(arr2D)
      ax.scatter(time_idx, frequency_idx)
      ax.set_xlabel('Time')
      ax.set_ylabel('Frequency')
      ax.set_title("Spectrogram")
      plt.gca().invert_yaxis()
      plt.show()

    return zip(frequency_idx, time_idx)

# Hash list structure: sha1_hash[0:20] time_offset
# example: [(e05b341a9b77a51fd26, 32), ... ] 

def slp_contour(fig, m, slp, lons, lats, window=100):
    """
    Add sea-level pressure labels to a contour map. I don't remember where I found the code for this function
    some time in the past, but I wish I could attribute it.

    :param fig:
    :param m:
    :param slp:
    :param lons:
    :param lats:
    :param window:
    :return:
    """
    def extrema(mat, mode='wrap', w=10):
        """
        Find the indices of local extrema (min and max)
        in the input array.
        """

        from scipy.ndimage.filters import minimum_filter, maximum_filter

        mn = minimum_filter(mat, size=w, mode=mode)
        mx = maximum_filter(mat, size=w, mode=mode)
        return np.nonzero(mat == mn), np.nonzero(mat == mx)

    caxisP = np.arange(900, 1050, 4)
    c2 = m.contour(lons, lats, slp, caxisP, latlon=True, linewidth=1.0, colors='black')
    plt.clabel(c2, c2.levels, inline=True, fmt='%0.0f')
    # Plot highs and lows for slp
    local_min, local_max = extrema(slp, mode='wrap', w=window)
    x, y = m(lons, lats)
    xlows = x[local_min]
    xhighs = x[local_max]
    ylows = y[local_min]
    yhighs = y[local_max]
    lowvals = slp[local_min]
    highvals = slp[local_max]
    # Plot lows
    xyplotted = []
    yoffset = 0.022 * (m.ymax - m.ymin)
    dmin = 20.0 * yoffset
    for x, y, p in zip(xlows, ylows, lowvals):
        if (m.xmax - dmin > x > m.xmin + dmin and m.ymax - dmin > y > m.ymin + dmin):
            dist = [np.sqrt((x - x0) ** 2 + (y - y0) ** 2) for x0, y0 in xyplotted]
            if not dist or min(dist) > dmin:
                plt.text(x, y, 'L', fontsize=14, fontweight='bold', ha='center', va='center', color='r')
                plt.text(x, y - yoffset, repr(int(p)), fontsize=9, ha='center', va='top', color='r',
                         bbox=dict(boxstyle="square", ec='None', fc=(1, 1, 1, 0.5)))
                xyplotted.append((x, y))
    # Plot highs
    xyplotted = []
    for x, y, p in zip(xhighs, yhighs, highvals):
        if (m.xmax - dmin > x > m.xmin + dmin and m.ymax - dmin > y > m.ymin + dmin):
            dist = [np.sqrt((x - x0) ** 2 + (y - y0) ** 2) for x0, y0 in xyplotted]
            if not dist or min(dist) > dmin:
                plt.text(x, y, 'H', fontsize=14, fontweight='bold', ha='center', va='center', color='b')
                plt.text(x, y - yoffset, repr(int(p)), fontsize=9, ha='center', va='top', color='b',
                         bbox=dict(boxstyle="square", ec='None', fc=(1, 1, 1, 0.5)))
                xyplotted.append((x, y))
    return fig 

def no_background_patches(threshold=0.4, percentile=99.9):

    """Returns a patch filter to be used by :func:`create_patches` to determine for each image pair which patches
    are eligible for sampling. The purpose is to only sample patches from "interesting" regions of the raw image that
    actually contain a substantial amount of non-background signal. To that end, a maximum filter is applied to the target image
    to find the largest values in a region.

    Parameters
    ----------
    threshold : float, optional
        Scalar threshold between 0 and 1 that will be multiplied with the (outlier-robust)
        maximum of the image (see `percentile` below) to denote a lower bound.
        Only patches with a maximum value above this lower bound are eligible to be sampled.
    percentile : float, optional
        Percentile value to denote the (outlier-robust) maximum of an image, i.e. should be close 100.

    Returns
    -------
    function
        Function that takes an image pair `(y,x)` and the patch size as arguments and
        returns a binary mask of the same size as the image (to denote the locations
        eligible for sampling for :func:`create_patches`). At least one pixel of the
        binary mask must be ``True``, otherwise there are no patches to sample.

    Raises
    ------
    ValueError
        Illegal arguments.
    """

    (np.isscalar(percentile) and 0 <= percentile <= 100) or _raise(ValueError())
    (np.isscalar(threshold)  and 0 <= threshold  <=   1) or _raise(ValueError())

    from scipy.ndimage.filters import maximum_filter
    def _filter(datas, patch_size, dtype=np.float32):
        image = datas[0]
        if dtype is not None:
            image = image.astype(dtype)
        # make max filter patch_size smaller to avoid only few non-bg pixel close to image border
        patch_size = [(p//2 if p>1 else p) for p in patch_size]
        filtered = maximum_filter(image, patch_size, mode='constant')
        return filtered > threshold * np.percentile(image,percentile)
    return _filter



## Sample patches