Python scipy.ndimage.filters.minimum_filter() Examples

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 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 get_distance(region, src):
    """
    Compute within-region distances from the src pixels.

    Parameters
    ----------
    region : np.ndarray(shape=(m, n), dtype=bool)
        mask of the region
    src : np.ndarray(shape=(m, n), dtype=bool)
        mask of the source pixels to compute distances from.

    Returns
    -------
    d : np.ndarray(shape=(m, n), dtype=float)
        approximate within-region distance from the nearest src pixel;
        (distances outside of the region are arbitrary).
    """

    dmax = float(region.size)
    d = np.full(region.shape, dmax)
    d[src] = 0
    for n in range(region.size):
        d_orth = minimum_filter(d, footprint=_ORTH2) + 1
        d_diag = minimum_filter(d, (3, 3)) + _SQRT2
        d_adj = np.minimum(d_orth[region], d_diag[region])
        d[region] = np.minimum(d_adj, d[region])
        if (d[region] < dmax).all():
            break
    return d 

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 morphOpen(V, footprint):
    ''' computes the morphological opening of V (correlation map) with circular footprint'''
    vrem   = filters.minimum_filter(V, footprint=footprint)
    vrem   = -filters.minimum_filter(-vrem, footprint=footprint)
    return vrem 

def find_centers(self,data):
        
        def fill_nan(A):
            #Interpolate to fill nan values
            A = np.array(A)
            inds = np.arange(len(A))
            good = np.where(np.isfinite(A))
            good_grad = np.gradient(good[0])
            if len(good[0])>=3:
                f = interp1d(inds[good], A[good],bounds_error=False,kind='quadratic')
                B = np.where(np.isfinite(A)[good[0][0]:good[0][-1]+1],
                             A[good[0][0]:good[0][-1]+1],
                             f(inds[good[0][0]:good[0][-1]+1]))
                return [np.nan]*good[0][0]+list(B)+[np.nan]*(inds[-1]-good[0][-1])
            else:
                return [np.nan]*len(A)
        
        #Check that sfc pressure spread is big enough to identify real minima
        if np.nanpercentile(data['p_sfc'],90)-np.nanpercentile(data['p_sfc'],10)>8:
            data['p_sfc'][:20]=[np.nan]*20 #NaN out the first 10 minutes of the flight
            p_sfc_interp = fill_nan(data['p_sfc']) #Interp p_sfc across missing data
            wspd_interp = fill_nan(data['wspd']) #Interp wspd across missing data
            #Smooth p_sfc and wspd
            p_sfc_smooth = [np.nan]*1+list(np.convolve(p_sfc_interp,[1/3]*3,mode='valid'))+[np.nan]*1
            wspd_smooth = [np.nan]*1+list(np.convolve(wspd_interp,[1/3]*3,mode='valid'))+[np.nan]*1
            #Add wspd to p_sfc to encourage finding p mins with wspd mins 
            #and prevent finding p mins in intense thunderstorms
            pw_test = np.array(p_sfc_smooth)+np.array(wspd_smooth)*.1
            #Find mins in 15-minute windows
            imin = np.nonzero(pw_test == minimum_filter(pw_test,30))[0]
            #Only use mins if below 15th %ile of mission p_sfc data and when plane p is 500-900mb
            imin = [i for i in imin if 800<p_sfc_interp[i]<np.nanpercentile(data['p_sfc'],15) and \
                    550<data['plane_p'][i]<900]
        else:
            imin=[]
        data['iscenter'] = np.zeros(len(data['p_sfc']))
        for i in imin:
            j = data.index.values[i]
            data['iscenter'][j] = 1
        return data 

def plot_fill_flat(roi, out, region, source, drain, dL, dH):
    from matplotlib import pyplot

    plot_detail = roi.size < 500
    cmap = 'Greens'

    pyplot.figure()

    ax = pyplot.subplot(221)
    pyplot.axis('off')
    pyplot.title('unfilled')
    im = pyplot.imshow(roi, interpolation='none')
    im.set_cmap(cmap)
    if plot_detail:
        y, x = np.where(region); pyplot.plot(x, y, 'k.')
        y, x = np.where(source); pyplot.plot(x, y, lw=0, color='k', marker='$H$', ms=12)
        y, x = np.where(drain);   pyplot.plot(x, y, lw=0, color='k', marker='$L$', ms=12)
    
    pyplot.subplot(222, sharex=ax, sharey=ax)
    pyplot.axis('off')
    pyplot.title('filled')
    im = pyplot.imshow(out, interpolation='none')
    im.set_cmap(cmap)
    if plot_detail:
        for elev in np.unique(out):
            y, x = np.where(out==elev)
            pyplot.plot(x, y, lw=0, color='k', marker='$%.3f$' % elev, ms=20)

    if plot_detail:
        flat = (minimum_filter(out, (3, 3)) >= out) & region
        y, x = np.where(flat); pyplot.plot(x, y, 'r_', ms=24)
        
        pyplot.subplot(223, sharex=ax, sharey=ax)
        pyplot.axis('off')
        pyplot.title('dL')
        im = pyplot.imshow(roi, interpolation='none')
        im.set_cmap(cmap)
        for d in np.unique(dL):
            if d == region.size: continue
            y, x = np.where(dL==d)
            pyplot.plot(x, y, lw=0, color='k', marker='$%.2f$' % d, ms=24)

        pyplot.subplot(224, sharex=ax, sharey=ax)
        pyplot.axis('off')
        pyplot.title('dH')
        im = pyplot.imshow(roi, interpolation='none')
        im.set_cmap(cmap)
        for d in np.unique(dH):
            if d == region.size: continue
            y, x = np.where(dH==d)
            pyplot.plot(x, y, lw=0, color='k', marker='$%.2f$' % d, ms=24)

    pyplot.tight_layout() 

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 sample_patches_from_multiple_stacks(datas, patch_size, n_samples, datas_mask=None, patch_filter=None, verbose=False):
    """ sample matching patches of size `patch_size` from all arrays in `datas` """

    # TODO: some of these checks are already required in 'create_patches'
    len(patch_size)==datas[0].ndim or _raise(ValueError())

    if not all(( a.shape == datas[0].shape for a in datas )):
        raise ValueError("all input shapes must be the same: %s" % (" / ".join(str(a.shape) for a in datas)))

    if not all(( 0 < s <= d for s,d in zip(patch_size,datas[0].shape) )):
        raise ValueError("patch_size %s negative or larger than data shape %s along some dimensions" % (str(patch_size), str(datas[0].shape)))

    if patch_filter is None:
        patch_mask = np.ones(datas[0].shape,dtype=np.bool)
    else:
        patch_mask = patch_filter(datas, patch_size)

    if datas_mask is not None:
        # TODO: Test this
        warnings.warn('Using pixel masks for raw/transformed images not tested.')
        datas_mask.shape == datas[0].shape or _raise(ValueError())
        datas_mask.dtype == np.bool or _raise(ValueError())
        from scipy.ndimage.filters import minimum_filter
        patch_mask &= minimum_filter(datas_mask, patch_size, mode='constant', cval=False)

    # get the valid indices

    border_slices = tuple([slice(s // 2, d - s + s // 2 + 1) for s, d in zip(patch_size, datas[0].shape)])
    valid_inds = np.where(patch_mask[border_slices])
    n_valid = len(valid_inds[0])

    if n_valid == 0:
        raise ValueError("'patch_filter' didn't return any region to sample from")

    sample_inds = choice(range(n_valid), n_samples, replace=(n_valid < n_samples))

    # valid_inds = [v + s.start for s, v in zip(border_slices, valid_inds)] # slow for large n_valid
    # rand_inds = [v[sample_inds] for v in valid_inds]
    rand_inds = [v[sample_inds] + s.start for s, v in zip(border_slices, valid_inds)]

    # res = [np.stack([data[r[0] - patch_size[0] // 2:r[0] + patch_size[0] - patch_size[0] // 2,
    #                  r[1] - patch_size[1] // 2:r[1] + patch_size[1] - patch_size[1] // 2,
    #                  r[2] - patch_size[2] // 2:r[2] + patch_size[2] - patch_size[2] // 2,
    #                  ] for r in zip(*rand_inds)]) for data in datas]

    res = [np.stack([data[tuple(slice(_r-(_p//2),_r+_p-(_p//2)) for _r,_p in zip(r,patch_size))] for r in zip(*rand_inds)]) for data in datas]

    return res



## Create training data