Python scipy.ndimage.convolve() Examples

def test_correlate14(self):
        kernel = numpy.array([[1, 0],
                              [0, 1]])
        for type1 in self.types:
            array = numpy.array([[1, 2, 3],
                                 [4, 5, 6]], type1)
            for type2 in self.types:
                output = numpy.zeros(array.shape, type2)
                ndimage.correlate(array, kernel,
                                  output=output)
                assert_array_almost_equal([[2, 3, 5], [5, 6, 8]], output)
                assert_equal(output.dtype.type, type2)

                ndimage.convolve(array, kernel, output=output)
                assert_array_almost_equal([[6, 8, 9], [9, 11, 12]], output)
                assert_equal(output.dtype.type, type2) 

def compute_img_filter_response2d(img, filter_battery):
    """ compute image filter response in 2D

    :param [[float]] img: image
    :param [[[float]]] filter_battery: filters
    :return [[float]]:
    """
    if filter_battery.ndim != 3:
        raise ValueError('wrong battery dim %r' % filter_battery.shape)
    responses = np.array([ndimage.convolve(img, fl) for fl in filter_battery])
    if filter_battery.shape[0] > 1:
        # usually for rotational edge detectors and we tae the maximal response
        response = np.max(responses, axis=0)
    else:
        response = responses[0]
    return response 

def get_data(self):

        image_iter = self.ds_images.get_data()
        psf_iter = self.ds_psf.get_data()

        for dp_image in image_iter:

            # sample camera shake kernel
            dp_psf = next(psf_iter)

            # synthesize ego-motion
            for t, k in enumerate(dp_psf):
                blurry = dp_image[t]
                for c in range(3):
                    blurry[:, :, c] = ndimage.convolve(blurry[:, :, c], k, mode='constant', cval=0.0)
                dp_image[t] = blurry

            yield dp_image 

def test_correlate14(self):
        kernel = numpy.array([[1, 0],
                              [0, 1]])
        for type1 in self.types:
            array = numpy.array([[1, 2, 3],
                                 [4, 5, 6]], type1)
            for type2 in self.types:
                output = numpy.zeros(array.shape, type2)
                ndimage.correlate(array, kernel,
                                  output=output)
                assert_array_almost_equal([[2, 3, 5], [5, 6, 8]], output)
                assert_equal(output.dtype.type, type2)

                ndimage.convolve(array, kernel, output=output)
                assert_array_almost_equal([[6, 8, 9], [9, 11, 12]], output)
                assert_equal(output.dtype.type, type2) 

def convolve_mask(data, ksize=3, kernel=None, copy=True):
    """
    Convolve data over the missing regions of a mask

    Parameters
    ----------
    data : masked array_like
        Input field.
    ksize : int, optional
        Size of square kernel
    kernel : ndarray, optional
        Define a convolution kernel. Default is averaging
    copy : bool, optional
        If true, a copy of input array is made

    Returns
    -------
    fld : masked array
    """
    return convolve(data, ksize, kernel, copy, True) 

def test_correlate03(self):
        array = numpy.array([1])
        weights = numpy.array([1, 1])
        expected = [2]

        output = ndimage.correlate(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.convolve(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.correlate1d(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.convolve1d(array, weights)
        assert_array_almost_equal(output, expected) 

def compute_gabor(img, params):
    if not params["shared_dict"].get("gabor_kernels", False):
        gabor_theta = int(params.get("gabor_theta", 4))
        gabor_sigma = make_tuple(params.get("gabor_sigma", "(1,3)"))
        gabor_frequency = make_tuple(params.get("gabor_frequency", "(0.05, 0.25)"))

        kernels = []
        for theta in range(gabor_theta):
            theta = theta / 4. * np.pi
            for sigma in gabor_sigma:
                for frequency in gabor_frequency:
                    kernel = np.real(gabor_kernel(frequency, theta=theta,
                                                  sigma_x=sigma, sigma_y=sigma))
                    kernels.append(kernel)
        params["shared_dict"]["gabor_kernels"] = kernels

    kernels = params["shared_dict"]["gabor_kernels"]
    imgg = rgb2gray(img)
    feats = np.zeros((imgg.shape[0], imgg.shape[1], len(kernels)), dtype=np.double)
    for k, kernel in enumerate(kernels):
        filtered = ndi.convolve(imgg, kernel, mode='wrap')
        feats[:, :, k] = filtered
    return feats 

def test_correlate01(self):
        array = numpy.array([1, 2])
        weights = numpy.array([2])
        expected = [2, 4]

        output = ndimage.correlate(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.convolve(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.correlate1d(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.convolve1d(array, weights)
        assert_array_almost_equal(output, expected) 

def sobel(image):
    # type: (numpy.ndarray) -> numpy.ndarray
    """
    Find edges using the Sobel approximation to the derivatives.

    Inspired by the [Octave implementation](https://sourceforge.net/p/octave/image/ci/default/tree/inst/edge.m#l196).
    """
    from skimage.filters.edges import HSOBEL_WEIGHTS
    h1 = np.array(HSOBEL_WEIGHTS)
    h1 /= np.sum(abs(h1))  # normalize h1
    
    from scipy.ndimage import convolve
    strength2 = np.square(convolve(image, h1.T))

    # Note: https://sourceforge.net/p/octave/image/ci/default/tree/inst/edge.m#l59
    thresh2 = 2 * np.sqrt(np.mean(strength2))

    strength2[strength2 <= thresh2] = 0
    return _simple_thinning(strength2) 

def calcul_gradient(self):
        rac_mnt = self.racine + '_' + str(self.res) + 'm'
        print(rac_mnt)
        fic_mnt = rac_mnt + 'float.mnt'
        fic_hdr = rac_mnt + 'float.hdr'
        fic_dz_dl = rac_mnt + 'float.dz_dl'
        fic_dz_dc = rac_mnt + 'float.dz_dc'
        (nblig, nbcol, type_donnee, endian) = lire_entete_mnt(fic_hdr)

        srtm = (np.fromfile(fic_mnt, type_donnee)).reshape(nblig, nbcol)
        Noyau_horizontal = np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]])
        Noyau_vertical = np.array([[1, 2, 1], [0, 0, 0], [-1, -2, -1]])

        dz_dc = nd.convolve(srtm, Noyau_horizontal) / 8. / self.res
        dz_dl = nd.convolve(srtm, Noyau_vertical) / 8. / self.res

        dz_dl.tofile(fic_dz_dl)
        dz_dc.tofile(fic_dz_dc)
        return (fic_dz_dl, fic_dz_dc)

    #############################################################
    ########################Calcul_pentes########################
    ############################################################# 

def test_correlate13(self):
        kernel = numpy.array([[1, 0],
                              [0, 1]])
        for type1 in self.types:
            array = numpy.array([[1, 2, 3],
                                 [4, 5, 6]], type1)
            for type2 in self.types:
                output = ndimage.correlate(array, kernel,
                                                    output=type2)
                assert_array_almost_equal([[2, 3, 5], [5, 6, 8]], output)
                assert_equal(output.dtype.type, type2)

                output = ndimage.convolve(array, kernel,
                                          output=type2)
                assert_array_almost_equal([[6, 8, 9], [9, 11, 12]], output)
                assert_equal(output.dtype.type, type2) 

def ApplyAtoms(V,D,scale):
    out=[]
    for s in xrange(scale):
        if s!=0:
            print('scale='+str(s))
            V = pyr.pyramid_reduce_3d(V,downscale=2) # reduce the volume. e.g. from 512^3 to 256^3
        else: print('scale=0')
        for i in xrange(len(D)):
            print('s:'+str(s)+' i:'+str(i))
            conv = nd.convolve(V, D[i], mode='constant', cval=0.0)
            if s==0:
                out.append(conv)
            else:
                upscaled = pyr.pyramid_expand_3d(conv, upscale=2**s)
                out.append(upscaled)
    out=np.array(out)
    return out 

def test_correlate19(self):
        kernel = numpy.array([[1, 0],
                              [0, 1]])
        for type1 in self.types:
            array = numpy.array([[1, 2, 3],
                                 [4, 5, 6]], type1)
            output = ndimage.correlate(array, kernel,
                                       output=numpy.float32,
                                       mode='nearest', origin=[-1, 0])
            assert_array_almost_equal([[5, 6, 8], [8, 9, 11]], output)
            assert_equal(output.dtype.type, numpy.float32)

            output = ndimage.convolve(array, kernel,
                                      output=numpy.float32,
                                      mode='nearest', origin=[-1, 0])
            assert_array_almost_equal([[3, 5, 6], [6, 8, 9]], output)
            assert_equal(output.dtype.type, numpy.float32) 

def high_pass_filter(input_image, output_image): 

    I = imread(input_image)
    if len(I.shape)==3:
        kernel = np.array([[[-1, -1, -1],
                            [-1,  8, -1],
                            [-1, -1, -1]],
                           [[-1, -1, -1],
                            [-1,  8, -1],
                            [-1, -1, -1]],
                           [[-1, -1, -1],
                            [-1,  8, -1],
                            [-1, -1, -1]]])
    else:
        kernel = np.array([[-1, -1, -1],
                           [-1,  8, -1],
                           [-1, -1, -1]])


    If = ndimage.convolve(I, kernel)
    imsave(output_image, If)
# }}}

# {{{ low_pass_filter() 

def pyconv3d(signals, filters):
    Ns, Ts, C, Hs, Ws = signals.shape
    Nf, Tf, C, Hf, Wf = filters.shape

    Tf2 = Tf//2
    Hf2 = Hf//2
    Wf2 = Wf//2

    rval = numpy.zeros((Ns, Ts-Tf+1, Nf, Hs-Hf+1, Ws-Wf+1))
    for ns in xrange(Ns):
        for nf in xrange(Nf):
            for c in xrange(C):
                s_i = signals[ns, :, c, :, :]
                f_i = filters[nf, :, c, :, :]
                r_i = rval[ns, :, nf, :, :]
                o_i = ndimage.convolve(s_i, f_i, mode='constant', cval=1)
                o_i_sh0 = o_i.shape[0]
                # print s_i.shape, f_i.shape, r_i.shape, o_i.shape
                r_i += o_i[Tf2:o_i_sh0-Tf2, Hf2:-Hf2, Wf2:-Wf2]
    return rval 

def test_correlate19(self):
        kernel = numpy.array([[1, 0],
                              [0, 1]])
        for type1 in self.types:
            array = numpy.array([[1, 2, 3],
                                 [4, 5, 6]], type1)
            output = ndimage.correlate(array, kernel,
                                       output=numpy.float32,
                                       mode='nearest', origin=[-1, 0])
            assert_array_almost_equal([[5, 6, 8], [8, 9, 11]], output)
            assert_equal(output.dtype.type, numpy.float32)

            output = ndimage.convolve(array, kernel,
                                      output=numpy.float32,
                                      mode='nearest', origin=[-1, 0])
            assert_array_almost_equal([[3, 5, 6], [6, 8, 9]], output)
            assert_equal(output.dtype.type, numpy.float32) 

def pyconv3d(signals, filters):
    Ns, Ts, C, Hs, Ws = signals.shape
    Nf, Tf, C, Hf, Wf = filters.shape

    Tf2 = Tf//2
    Hf2 = Hf//2
    Wf2 = Wf//2

    rval = numpy.zeros((Ns, Ts-Tf+1, Nf, Hs-Hf+1, Ws-Wf+1))
    for ns in xrange(Ns):
        for nf in xrange(Nf):
            for c in xrange(C):
                s_i = signals[ns, :, c, :, :]
                f_i = filters[nf, :, c, :, :]
                r_i = rval[ns, :, nf, :, :]
                o_i = ndimage.convolve(s_i, f_i, mode='constant', cval=1)
                o_i_sh0 = o_i.shape[0]
                # print s_i.shape, f_i.shape, r_i.shape, o_i.shape
                r_i += o_i[Tf2:o_i_sh0-Tf2, Hf2:-Hf2, Wf2:-Wf2]
    return rval 

def get_ini_cor(cor_img, d1=21, d2=3):
    cor = convolve(cor_img, np.ones((d1, d1)), mode='constant', cval=0.0)
    cor_id = []
    X_loc = find_N_peaks(cor.sum(0), prominence=None,
                         distance=20, N=4)[0]
    for x in X_loc:
        x_ = int(np.round(x))

        V_signal = cor[:, max(0, x_-d2):x_+d2+1].sum(1)
        y1, y2 = find_N_peaks(V_signal, prominence=None,
                              distance=20, N=2)[0]
        cor_id.append((x, y1))
        cor_id.append((x, y2))

    cor_id = np.array(cor_id, np.float64)

    return cor_id 

def test_correlate03(self):
        array = numpy.array([1])
        weights = numpy.array([1, 1])
        expected = [2]

        output = ndimage.correlate(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.convolve(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.correlate1d(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.convolve1d(array, weights)
        assert_array_almost_equal(output, expected) 

def get_ini_cor(cor_img, d1=21, d2=3):
    cor = convolve(cor_img, np.ones((d1, d1)), mode='constant', cval=0.0)
    cor_id = []
    X_loc = find_N_peaks(cor.sum(0), prominence=None,
                         distance=20, N=4)[0]
    for x in X_loc:
        x_ = int(np.round(x))

        V_signal = cor[:, max(0, x_-d2):x_+d2+1].sum(1)
        y1, y2 = find_N_peaks(V_signal, prominence=None,
                              distance=20, N=2)[0]
        cor_id.append((x, y1))
        cor_id.append((x, y2))

    cor_id = np.array(cor_id, np.float64)

    return cor_id 

def test_correlate01(self):
        array = numpy.array([1, 2])
        weights = numpy.array([2])
        expected = [2, 4]

        output = ndimage.correlate(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.convolve(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.correlate1d(array, weights)
        assert_array_almost_equal(output, expected)

        output = ndimage.convolve1d(array, weights)
        assert_array_almost_equal(output, expected) 

def get_ini_cor(cor_img, d1=21, d2=3):
    cor = convolve(cor_img, np.ones((d1, d1)), mode='constant', cval=0.0)
    cor_id = []
    X_loc = find_N_peaks(cor.sum(0), prominence=None,
                         distance=20, N=4)[0]
    for x in X_loc:
        x_ = int(np.round(x))

        V_signal = cor[:, max(0, x_-d2):x_+d2+1].sum(1)
        y1, y2 = find_N_peaks(V_signal, prominence=None,
                              distance=20, N=2)[0]
        cor_id.append((x, y1))
        cor_id.append((x, y2))

    cor_id = np.array(cor_id, np.float64)

    return cor_id 

def low_pass_filter(input_image, output_image): 

    I = imread(input_image)
    if len(I.shape)==3:
        kernel = np.array([[[1, 1, 1],
                            [1, 1, 1],
                            [1, 1, 1]],
                           [[1, 1, 1],
                            [1, 1, 1],
                            [1, 1, 1]],
                           [[1, 1, 1],
                            [1, 1, 1],
                            [1, 1, 1]]])
    else:
        kernel = np.array([[1, 1, 1],
                           [1, 1, 1],
                           [1, 1, 1]])

    kernel = kernel.astype('float32')/9
    If = ndimage.convolve(I, kernel)
    imsave(output_image, If)
# }}}

# {{{ imgdiff() 

def circular_conv(signal1, signal2):
    """takes two 1D numpy array and do a circular convolution with them
    inputs :
        - signal1 : 1D numpy array
        - signal2 : 1D numpy array, same length as signal1
    output :
        - signal_conv : 1D numpy array, result of circular convolution of signal1 and signal2"""

    if signal1.shape != signal2.shape:
        raise Exception("The two signals for circular convolution do not have the same shape")

    signal2_extended = np.concatenate((signal2, signal2, signal2))  # replicate signal at both ends

    signal_conv_extended = np.convolve(signal1, signal2_extended, mode="same")  # median filtering

    signal_conv = signal_conv_extended[len(signal1):2*len(signal1)]  # truncate back the signal

    return signal_conv 

def test_correlate15(self):
        kernel = numpy.array([[1, 0],
                              [0, 1]])
        for type1 in self.types:
            array = numpy.array([[1, 2, 3],
                                 [4, 5, 6]], type1)
            output = ndimage.correlate(array, kernel,
                                       output=numpy.float32)
            assert_array_almost_equal([[2, 3, 5], [5, 6, 8]], output)
            assert_equal(output.dtype.type, numpy.float32)

            output = ndimage.convolve(array, kernel,
                                      output=numpy.float32)
            assert_array_almost_equal([[6, 8, 9], [9, 11, 12]], output)
            assert_equal(output.dtype.type, numpy.float32) 

def test_correlate11(self):
        array = numpy.array([[1, 2, 3],
                             [4, 5, 6]])
        kernel = numpy.array([[1, 1],
                              [1, 1]])
        output = ndimage.correlate(array, kernel)
        assert_array_almost_equal([[4, 6, 10], [10, 12, 16]], output)
        output = ndimage.convolve(array, kernel)
        assert_array_almost_equal([[12, 16, 18], [18, 22, 24]], output) 

def test_kernel_2D_as_kernel_3D(self):

        ## Shape in (z,y,x)-coordinates
        nda_shape = (40, 200, 200)

        # Define size of kernel
        N = 6

        # Create random data array
        nda = 255 * np.random.rand(nda_shape[0], nda_shape[1], nda_shape[2])

        # Create kernel with elements 1:N^2
        kernel = np.arange(1, N*N+1)

        # Define 2D- and equivalent 3D-kernel
        kernel_2D = kernel.reshape(N, N)
        kernel_3D = kernel.reshape(1, N, N)

        # Create data array copies
        nda_2D = np.array(nda)
        nda_3D = np.array(nda)

        # 2D
        for i in range(0, nda.shape[0]):
            nda_2D[i, :, :] = ndimage.convolve(nda_2D[i, :, :], kernel_2D)

        # 3D
        nda_3D = ndimage.convolve(nda_3D, kernel_3D)

        self.assertEqual(np.around(
            np.linalg.norm(nda_2D-nda_3D), decimals=self.accuracy), 0) 

def nudge_dataset(X, Y):
    """
    This produces a dataset 5 times bigger than the original one,
    by moving the 8x8 images in X around by 1px to left, right, down, up
    """
    direction_vectors = [
        [[0, 1, 0],
         [0, 0, 0],
         [0, 0, 0]],

        [[0, 0, 0],
         [1, 0, 0],
         [0, 0, 0]],

        [[0, 0, 0],
         [0, 0, 1],
         [0, 0, 0]],

        [[0, 0, 0],
         [0, 0, 0],
         [0, 1, 0]]]

    shift = lambda x, w: convolve(x.reshape((8, 8)), mode='constant',
                                  weights=w).ravel()
    X = np.concatenate([X] +
                       [np.apply_along_axis(shift, 1, X, vector)
                        for vector in direction_vectors])
    Y = np.concatenate([Y for _ in range(5)], axis=0)
    return X, Y 

def test_correlate17(self):
        array = numpy.array([1, 2, 3])
        tcor = [3, 5, 6]
        tcov = [2, 3, 5]
        kernel = numpy.array([1, 1])
        output = ndimage.correlate(array, kernel, origin=-1)
        assert_array_almost_equal(tcor, output)
        output = ndimage.convolve(array, kernel, origin=-1)
        assert_array_almost_equal(tcov, output)
        output = ndimage.correlate1d(array, kernel, origin=-1)
        assert_array_almost_equal(tcor, output)
        output = ndimage.convolve1d(array, kernel, origin=-1)
        assert_array_almost_equal(tcov, output) 

def test_correlate16(self):
        kernel = numpy.array([[0.5, 0],
                              [0, 0.5]])
        for type1 in self.types:
            array = numpy.array([[1, 2, 3], [4, 5, 6]], type1)
            output = ndimage.correlate(array, kernel, output=numpy.float32)
            assert_array_almost_equal([[1, 1.5, 2.5], [2.5, 3, 4]], output)
            assert_equal(output.dtype.type, numpy.float32)

            output = ndimage.convolve(array, kernel, output=numpy.float32)
            assert_array_almost_equal([[3, 4, 4.5], [4.5, 5.5, 6]], output)
            assert_equal(output.dtype.type, numpy.float32)