Python ndimage.convolve1d函数代码示例

def test_sg_coeffs_exact():
    polyorder = 4
    window_length = 9
    halflen = window_length // 2

    x = np.linspace(0, 21, 43)
    delta = x[1] - x[0]

    # The data is a cubic polynomial.  We'll use an order 4
    # SG filter, so the filtered values should equal the input data
    # (except within half window_length of the edges).
    y = 0.5 * x ** 3 - x
    h = savgol_coeffs(window_length, polyorder)
    y0 = convolve1d(y, h)
    assert_allclose(y0[halflen:-halflen], y[halflen:-halflen])

    # Check the same input, but use deriv=1.  dy is the exact result.
    dy = 1.5 * x ** 2 - 1
    h = savgol_coeffs(window_length, polyorder, deriv=1, delta=delta)
    y1 = convolve1d(y, h)
    assert_allclose(y1[halflen:-halflen], dy[halflen:-halflen])

    # Check the same input, but use deriv=2. d2y is the exact result.
    d2y = 3.0 * x
    h = savgol_coeffs(window_length, polyorder, deriv=2, delta=delta)
    y2 = convolve1d(y, h)
    assert_allclose(y2[halflen:-halflen], d2y[halflen:-halflen])

def move_nanvar_filter(arr, window, axis=-1):
    "Moving window variance ignoring NaNs, implemented with a filter."
    arr = np.array(arr, copy=False)
    global convolve1d
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis is None:
        raise ValueError("An `axis` value of None is not supported.")
    if window < 1:
        raise ValueError("`window` must be at least 1.")
    if window > arr.shape[axis]:
        raise ValueError("`window` is too long.")
    arr = arr.astype(float)
    nrr = np.isnan(arr)
    arr[nrr] = 0
    nrr = nrr.astype(int)
    w = np.ones(window, dtype=int)
    x0 = (1 - window) // 2
    convolve1d(nrr, w, axis=axis, mode="constant", cval=0, origin=x0, output=nrr)
    y = convolve1d(arr, w, axis=axis, mode="constant", cval=np.nan, origin=x0)
    y /= window - nrr
    y *= y
    arr *= arr
    convolve1d(arr, w, axis=axis, mode="constant", cval=np.nan, origin=x0, output=arr)
    arr /= window - nrr
    arr -= y
    arr[nrr == window] = np.nan
    return arr

 def test_03_02_adaptive_threshold_different(self):
     r = np.random.RandomState()
     r.seed(31)
     block = r.uniform(size=(10,10))
     i,j = np.mgrid[0:10:2,0:10:2]
     block[i,j] *= .5
     i,j = np.mgrid[0:50,0:50]
     img = block[i%10, j%10] * .5
     #
     # Make the middle higher in intensity
     #
     img[20:30, 20:30] *= 2
     global_threshold = T.get_global_threshold(T.TM_OTSU, block)
     adaptive_threshold = T.get_adaptive_threshold(
         T.TM_OTSU, img, global_threshold,
         adaptive_window_size = 10)
     #
     # Check that the gradients are positive for i,j<15 and negative
     # for i,j>=15
     #
     gradient = convolve1d(adaptive_threshold, [-1, 0, 1], 0)
     self.assertTrue(np.all(gradient[20:25, 20:30] < 0))
     self.assertTrue(np.all(gradient[25:30, 20:30] > 0))
     gradient = convolve1d(adaptive_threshold, [-1, 0, 1], 1)
     self.assertTrue(np.all(gradient[20:30, 20:25] < 0))
     self.assertTrue(np.all(gradient[20:30, 25:30] > 0))

	def conv_2_sep_Y_circ(self,image,xkernel,thetakernel):#,upsample=10):
		# up_image = np.repeat(image,upsample)
		# x_convolved = ndimage.convolve1d(image,xkernel,mode='mirror',axis=1)
		# both_convolved = ndimage.convolve1d(x_convolved,thetakernel,mode='mirror',axis=0)
		x_convolved = ndimage.convolve1d(image,xkernel,mode='constant',axis=1)
		both_convolved = ndimage.convolve1d(x_convolved,thetakernel,mode='wrap',axis=0)
		return both_convolved

def unsharp_masking(X):
    lp = np.array(X)
    for i, ws in zip([0, 1, 2], [50, 50, 25]):
        h = hamming(ws)
        h /= h.sum()
        convolve1d(lp, h, axis=i, output=lp)
    return X - lp

def second_derivatives(array, smooth=2):
    """
    Compute the second derivatives of all dimensions pairs of the input array
    
    :Inputs:
        array:  any ndarray
        smooth: the second derivative are computed by convolving the input array
                by [1,-2,1]. The smooth parameter set how many times this basic
                filter is convoluted with [1,2,1]/2, which smooth it.
    
    :Output:
        Return a tuple of the second derivative arrays in the order (where dd_ij 
        is the the second derivative d2(array)/didj for a N-dimensional array):
        (dd_00, dd_01, ..., dd_0N, dd_11, dd_12, ..., dd_1N, ..., dd_N-1N, dd_NN)

    :Example:
       for 3d array 'volume_array'
       dv_00, dv_01, dv_02, dv_11, dv_12, dv_22 = second_derivatives(volume_array)
       
    See also:
      numpy.gradient    
    """
    # compute the derivative filter
    dd  = [1,-1]
    for i in xrange(smooth):
        dd = _np.convolve(dd,[1,2,1])/2. 
    
    # compute the second derivatives
    res = ()
    for i in xrange(array.ndim):
        tmp = _nd.convolve1d(array,dd,axis=i)
        for j in xrange(i,array.ndim):
            res += _nd.convolve1d(tmp,dd,axis=j),
    
    return res

def move_var_filter(arr, window, axis=-1):
    "Moving window variance implemented with a filter."
    global convolve1d
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis == None:
        raise ValueError, "An `axis` value of None is not supported."
    if window < 1:  
        raise ValueError, "`window` must be at least 1."
    if window > arr.shape[axis]:
        raise ValueError, "`window` is too long."  
    arr = arr.astype(float)
    w = np.empty(window)
    w.fill(1.0 / window)
    x0 = (1 - window) // 2
    y = convolve1d(arr, w, axis=axis, mode='constant', cval=np.nan, origin=x0)
    y *= y
    arr *= arr
    convolve1d(arr, w, axis=axis, mode='constant', cval=np.nan, origin=x0,
               output=arr)
    arr -= y 
    return arr

def extract_oriented_patches2D( img, r, coordinates, nb_proj=100 ):

    img = img.astype('float64')

    projection_matrix = np.zeros( (nb_proj,2),
                                  dtype='float64' )
    for i in xrange(nb_proj):
        theta = float(i) * 2.0 * math.pi / nb_proj
        projection_matrix[i,0] = math.sin(theta)
        projection_matrix[i,1] = math.cos(theta)

    print "computing gradients..."
    # grad = np.dstack( ( nd.sobel( img, mode='constant', axis=0 ),
    #                     nd.sobel( img, mode='constant', axis=1 ) ) )

    grad = np.dstack( ( nd.convolve1d( img, [-1,1], mode='constant', axis=0 ),
                        nd.convolve1d( img, [-1,1], mode='constant', axis=1 ) ) )
    
    print "projecting gradients..."
    hist = _patches.project_gradients2D( grad, projection_matrix )
    hist = integral_image2D( hist )

    print hist

    print "extracting patches..."
    Y = coordinates[:,0].copy().astype('int32')
    X = coordinates[:,1].copy().astype('int32')
    return _patches.extract_oriented_patches2D( img,
                                                hist,
                                                projection_matrix,
                                                X,
                                                Y,
                                                r )

def move_nanmax_filter(arr, window, axis=-1):
    "Moving window maximium ignoring NaNs, implemented with a filter."
    global maximum_filter1d, convolve1d
    if maximum_filter1d is None:
        try:
            from scipy.ndimage import maximum_filter1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis == None:
        raise ValueError, "An `axis` value of None is not supported."
    if window < 1:  
        raise ValueError, "`window` must be at least 1."
    if window > arr.shape[axis]:
        raise ValueError, "`window` is too long."
    arr = arr.astype(float)
    nrr = np.isnan(arr)
    arr[nrr] = -np.inf
    x0 = (window - 1) // 2
    maximum_filter1d(arr, window, axis=axis, mode='constant', cval=np.nan,
                     origin=x0, output=arr)
    w = np.ones(window, dtype=int)
    nrr = nrr.astype(int)
    x0 = (1 - window) // 2
    convolve1d(nrr, w, axis=axis, mode='constant', cval=0, origin=x0,
               output=nrr)
    arr[nrr == window] = np.nan
    return arr

def seperable_gaussian_convolution():
    F = plt.imread('../images/cameraman.png')

    # Example of Gaussian convolution.
    s = 3
    G = convolve1d(F, Gauss1(s), axis=-1, mode='nearest')
    G = convolve1d(G, Gauss1(s), axis=0, mode='nearest')
    plt.subplot(1, 2, 1)
    plt.imshow(F, cmap=plt.cm.gray)
    plt.subplot(1, 2, 2)
    plt.imshow(G, cmap=plt.cm.gray)
    plt.show()

    # Gaussian convolution function run 100 times for average execution time.
    n = 10
    s_values = np.array([1, 2, 3, 5, 7, 9, 11, 15, 19])
    times = avg_runnning(F, s_values, n)
    times2d = gauss2d.avg_runnning(F, s_values, n)
    plt.clf()
    plt.title("Average execution time versus the scale\n(run {0} times)"
              .format(n))
    plt.xlabel("s")
    plt.ylabel("time (ms)")
    plt.plot(s_values, times, label="1d convolution")
    plt.plot(s_values, times2d, label="2d convolution")
    plt.legend()
    plt.show()

def move_nansum_filter(arr, window, axis=-1):
    """
    Moving sum (ignoring NaNs) along specified axis using the filter method.
    
    Parameters
    ----------
    arr : array_like
        Input array.
    window : int
        The number of elements in the moving window.
    axis : int, optional
        The axis over which to perform the moving sum. By default the moving
        sum is taken over the last axis (-1).
    
    Returns
    -------
    y : ndarray
        The moving sum (ignoring NaNs) of the input array along the specified
        axis.(A window with all NaNs returns NaN for the window sum.) The
        output has the same shape as the input.

    Notes
    -----
    The calculation of the sums uses scipy.ndimage.convolve1d. 

    Examples
    --------
    >>> from bottlechest.slow.move import move_sum_filter
    >>> arr = np.array([1, 2, np.nan, 4, 5, 6, 7])
    >>> move_nansum_filter(arr, window=2, axis=0)
    array([ NaN,   3.,   2.,   4.,   9.,  11.,  13.])

    """
    arr = np.array(arr, copy=False)
    global convolve1d
    if convolve1d is None:
        try:
            from scipy.ndimage import convolve1d
        except ImportError:
            raise ValueError("'filter' method requires SciPy.")
    if axis == None:
        raise ValueError("An `axis` value of None is not supported.")
    if window < 1:  
        raise ValueError("`window` must be at least 1.")
    if window > arr.shape[axis]:
        raise ValueError("`window` is too long.")
    arr = arr.astype(float)
    nrr = np.isnan(arr)
    arr[nrr] = 0
    nrr = nrr.astype(int)
    w = np.ones(window, dtype=int)
    x0 = (1 - window) // 2
    convolve1d(arr, w, axis=axis, mode='constant', cval=np.nan, origin=x0,
               output=arr)
    convolve1d(nrr, w, axis=axis, mode='constant', cval=0, origin=x0,
               output=nrr)
    arr[nrr == window] = np.nan
    return arr

def main():
    s = 2.0
    m = gauss_1(s)

    img = imread('cameraman.png')
    img2 = convolve1d(img, m, axis=0, mode='nearest')
    img2 = convolve1d(img2, m, axis=1, mode='nearest')
    imshow(img2, cmap=cm.gray)
    show()

def simple_gradient(Im):
    g1=np.array([-0.5, 0, 0.5])
    
    Imx=nd.convolve1d(Im, g1, axis=0)
    Imy=nd.convolve1d(Im, g1, axis=1)
    
    ImMag=(Imx**2 +Imy**2)**0.5
    
    return ImMag, Imx, Imy

def local_standardize(X, kernelsize=(17, 17, 15)):
    local_sq = X ** 2
    local_mean = np.asarray(X)
    for axis, ks in enumerate(kernelsize):
        # w = np.ones(ks) / ks
        w = np.hamming(ks)
        w /= w.sum()
        local_sq = convolve1d(local_sq, w, axis=axis, mode="reflect")
        local_mean = convolve1d(local_mean, w, axis=axis, mode="reflect")
    return (X - local_mean) / np.sqrt(local_sq - local_mean ** 2)

def GodinTypeFilter(data, n, axis=0):
    ''' perform 3 times moving average over the specified array axis.
    suitable for time averaging
    '''
    weights = sp.ones((n),sp.float32)/n
    weights2 = sp.ones((n+1),sp.float32)/(n+1)
    data=nd.convolve1d(data, weights, axis=axis, mode='constant')
    data=nd.convolve1d(data, weights, axis=axis, mode='constant')
    data=nd.convolve1d(data, weights2, axis=axis, mode='constant')
    return data