Python scipy.ndimage.gaussian_filter() Examples

def printlike_fibrous(image, blur=1.0, blotches=5e-5, inverted=None):
    if inverted:
        selector = image
    elif inverted is None:
        selector = autoinvert(image)
    else:
        selector = 1 - image

    selector = random_blotches(selector, 3*blotches, blotches)
    paper = make_multiscale_noise(image.shape, [1.0, 5.0, 10.0, 50.0], weights=[1.0, 0.3, 0.5, 0.3], span=(0.7, 1.0))
    paper -= make_fibrous_image(image.shape, 300, 500, 0.01, span=(0.0, 0.25), blur=0.5)
    ink = make_multiscale_noise(image.shape, [1.0, 5.0, 10.0, 50.0], span=(0.0, 0.5))
    blurred = ndi.gaussian_filter(selector, blur)
    printed = blurred * ink + (1-blurred) * paper
    if inverted:
        return 1 - printed
    else:
        return printed 

def GaussianFilter(in_dem, sigma=1, out_file=None):

    print("Gaussian filtering ...")
    start_time = time.time()
    dem = rd.LoadGDAL(in_dem)
    no_data = dem.no_data
    projection = dem.projection
    geotransform = dem.geotransform

    gau = ndimage.gaussian_filter(dem, sigma=sigma)
    gau = np2rdarray(gau, no_data, projection, geotransform)
    print("Run time: {:.4f} seconds".format(time.time() - start_time))

    if out_file is not None:
        print("Saving dem ...")
        rd.SaveGDAL(out_file, gau)
        return out_file

    return gau


# #####################################  main script 

def ssim_exact(img1, img2, sd=1.5, C1=0.01**2, C2=0.03**2):

    mu1 = gaussian_filter(img1, sd)
    mu2 = gaussian_filter(img2, sd)
    mu1_sq = mu1 * mu1
    mu2_sq = mu2 * mu2
    mu1_mu2 = mu1 * mu2
    sigma1_sq = gaussian_filter(img1 * img1, sd) - mu1_sq
    sigma2_sq = gaussian_filter(img2 * img2, sd) - mu2_sq
    sigma12 = gaussian_filter(img1 * img2, sd) - mu1_mu2

    ssim_num = ((2 * mu1_mu2 + C1) * (2 * sigma12 + C2))

    ssim_den = ((mu1_sq + mu2_sq + C1) * (sigma1_sq + sigma2_sq + C2))

    ssim_map = ssim_num / ssim_den
    return numpy.mean(ssim_map) 

def subtract_background_dog(z, sigma_min, sigma_max):
    """Difference of gaussians method for background removal.

    Parameters
    ----------
    sigma_max : float
        Large gaussian blur sigma.
    sigma_min : float
        Small gaussian blur sigma.

    Returns
    -------
        Denoised diffraction pattern as np.array
    """
    blur_max = ndi.gaussian_filter(z, sigma_max)
    blur_min = ndi.gaussian_filter(z, sigma_min)

    return np.maximum(np.where(blur_min > blur_max, z, 0) - blur_max, 0) 

def circular_filter_1d(signal, window_size, kernel='gaussian'):

    """ This function filters circularly the signal inputted with a median filter of inputted size, in this context
    circularly means that the signal is wrapped around and then filtered
    inputs :
        - signal : 1D numpy array
        - window_size : size of the kernel, an int
    outputs :
        - signal_smoothed : 1D numpy array, same size as signal"""

    signal_extended = np.concatenate((signal, signal, signal))  # replicate signal at both ends
    if kernel == 'gaussian':
        signal_extended_smooth = ndimage.gaussian_filter(signal_extended, window_size)  # gaussian
    elif kernel == 'median':
        signal_extended_smooth = medfilt(signal_extended, window_size)  # median filtering
    else:
        raise Exception("Unknow type of kernel")

    signal_smoothed = signal_extended_smooth[len(signal):2*len(signal)]  # truncate back the signal

    return signal_smoothed 

def random_blobs(shape, blobdensity, size, roughness=2.0):
    from random import randint
    from builtins import range  # python2 compatible
    h, w = shape
    numblobs = int(blobdensity * w * h)
    mask = np.zeros((h, w), 'i')
    for i in range(numblobs):
        mask[randint(0, h-1), randint(0, w-1)] = 1
    dt = ndi.distance_transform_edt(1-mask)
    mask =  np.array(dt < size, 'f')
    mask = ndi.gaussian_filter(mask, size/(2*roughness))
    mask -= np.amin(mask)
    mask /= np.amax(mask)
    noise = np.random.rand(h, w)
    noise = ndi.gaussian_filter(noise, size/(2*roughness))
    noise -= np.amin(noise)
    noise /= np.amax(noise)
    return np.array(mask * noise > 0.5, 'f') 

def find_beam_center_blur(z, sigma):
    """Estimate direct beam position by blurring the image with a large
    Gaussian kernel and finding the maximum.

    Parameters
    ----------
    sigma : float
        Sigma value for Gaussian blurring kernel.

    Returns
    -------
    center : np.array
        np.array containing indices of estimated direct beam positon.
    """
    blurred = ndi.gaussian_filter(z, sigma, mode="wrap")
    center = np.unravel_index(blurred.argmax(), blurred.shape)
    return np.array(center) 

def threshold_segmentation(img):
    # calculate the overview level size and retrieve the image
    img_hsv = img.convert('HSV')
    img_hsv_np = np.array(img_hsv)

    # dilate image and then threshold the image
    schannel = img_hsv_np[:, :, 1]
    mask = np.zeros(schannel.shape)

    schannel = dilation(schannel, star(3))
    schannel = ndimage.gaussian_filter(schannel, sigma=(5, 5), order=0)
    threshold_global = threshold_otsu(schannel)

    mask[schannel > threshold_global] = FOREGROUND
    mask[schannel <= threshold_global] = BACKGROUND

    return mask 

def tf(self, img, k=0):
        if self.num > 0 and k >= self.num:
            return img

        # image is nhwtc
        for n in range(img.shape[0]):
            sig = self.sigma.sample()
            # sample each channel saperately to avoid correlations
            if sig > self.eps:
                if len(img.shape) == self.dim+2:
                    C = img.shape[-1]
                    for c in range(C):
                        img[n,..., c] = ndimage.gaussian_filter(img[n, ..., c], sig)
                elif len(img.shape) == self.dim+1:
                    img[n] = ndimage.gaussian_filter(img[n], sig)
                else:
                    raise ValueError('image shape is not supported')

        return img 

def get_saliency_for_shallownet(image_url,sal_url):
    arr_files = glob.glob(image_url+"*.jpg")
    for i in range(len(arr_files)):  
        url_image = arr_files[i]
        image = io.imread(url_image)       
        img = misc.imresize(image,(96,96))
        img = np.asarray(img, dtype = 'float32') / 255.
        img = img.transpose(2,0,1).reshape(3, 96, 96)
        xt = np.zeros((1, 3, 96, 96), dtype='float32')
        xt[0]=img
        y = juntingnet.predict(xt)
        tmp = y.reshape(48,48)
        blured= ndimage.gaussian_filter(tmp, sigma=3)
        sal_map = cv2.resize(tmp,(image.shape[1],image.shape[0]))
        sal_map -= np.min(sal_map)
        sal_map /= np.max(sal_map)
        #saliency = misc.imresize(y,(img.shape[0],img.shape[1]))
        aux = url_image.split("/")[-1].split(".")[0]
        misc.imsave(sal_url+'/'+aux+'.png', sal_map) 

def test_generic_laplace01(self):
        def derivative2(input, axis, output, mode, cval, a, b):
            sigma = [a, b / 2.0]
            input = numpy.asarray(input)
            order = [0] * input.ndim
            order[axis] = 2
            return ndimage.gaussian_filter(input, sigma, order,
                                           output, mode, cval)
        for type in self.types:
            array = numpy.array([[3, 2, 5, 1, 4],
                                    [5, 8, 3, 7, 1],
                                    [5, 6, 9, 3, 5]], type)
            output = numpy.zeros(array.shape, type)
            tmp = ndimage.generic_laplace(array, derivative2,
                    extra_arguments=(1.0,), extra_keywords={'b': 2.0})
            ndimage.gaussian_laplace(array, 1.0, output)
            assert_array_almost_equal(tmp, output) 

def test_gaussian_gradient_magnitude02(self):
        for type in [numpy.int32, numpy.float32, numpy.float64]:
            array = numpy.array([[3, 2, 5, 1, 4],
                                    [5, 8, 3, 7, 1],
                                    [5, 6, 9, 3, 5]], type) * 100
            tmp1 = ndimage.gaussian_filter(array, 1.0, [1, 0])
            tmp2 = ndimage.gaussian_filter(array, 1.0, [0, 1])
            output = numpy.zeros(array.shape, type)
            ndimage.gaussian_gradient_magnitude(array, 1.0,
                                                           output)
            expected = tmp1 * tmp1 + tmp2 * tmp2
            expected = numpy.sqrt(expected).astype(type)
            assert_array_almost_equal(expected, output) 

def filter(self, signal):
        """Filter a given signal with a choice of filter type (self.lres_filter).
        """
        signal = signal.copy()
        filter_size = [1, self.downsamp_t*2-1, self.downsamp_xz*2-1, self.downsamp_xz*2-1]

        if self.lres_filter == 'none' or (not self.lres_filter):
            output = signal
        elif self.lres_filter == 'gaussian':
            sigma = [0, int(self.downsamp_t/2), int(self.downsamp_xz/2), int(self.downsamp_xz/2)]
            output = ndimage.gaussian_filter(signal, sigma=sigma)
        elif self.lres_filter == 'uniform':
            output = ndimage.uniform_filter(signal, size=filter_size)
        elif self.lres_filter == 'median':
            output = ndimage.median_filter(signal, size=filter_size)
        elif self.lres_filter == 'maximum':
            output = ndimage.maximum_filter(signal, size=filter_size)
        else:
            raise NotImplementedError(
                "lres_filter must be one of none/gaussian/uniform/median/maximum")
        return output 

def gaussian_filter_with_nan(U,sigma):
    """
    Apply Gaussian filter when the data contain NaN
    INPUT:
        U           : a 2-D array (matrix)
        sigma       : std for the Gaussian kernel
    OUTPUT:
        Z           : filtered matrix
    """
    V=U.copy()
    V[np.isnan(U)]=0
    VV= ndimage.gaussian_filter(V,sigma=sigma)

    W=0*U.copy()+1
    W[np.isnan(U)]=0
    WW= ndimage.gaussian_filter(W,sigma=sigma)
    Z=VV/WW
    return(Z)

#===============================================================================
#=============================================================================== 

def chunk_lcn(chunk, sigma_mean, sigma_std, std_bias=0.0, rescale=1.0):
    """
    based on matlab code by Guanglei Xiong, see http://www.mathworks.com/matlabcentral/fileexchange/8303-local-normalization
    assuming chunk.shape == (num_examples, x, y, channels)

    'rescale' is an additional rescaling constant to get the variance of the result in the 'right' range.
    """
    means = np.zeros(chunk.shape, dtype=chunk.dtype)
    for k in xrange(len(chunk)):
        means[k] = skimage.filter.gaussian_filter(chunk[k], sigma_mean, multichannel=True)

    chunk = chunk - means # centering
    del means # keep memory usage in check

    variances = np.zeros(chunk.shape, dtype=chunk.dtype)
    chunk_squared = chunk**2
    for k in xrange(len(chunk)):
        variances[k] = skimage.filter.gaussian_filter(chunk_squared[k], sigma_std, multichannel=True)

    chunk = chunk / np.sqrt(variances + std_bias)

    return chunk / rescale

    # TODO: make this 100x faster lol. otherwise it's not usable. 

def test_generic_laplace01(self):
        def derivative2(input, axis, output, mode, cval, a, b):
            sigma = [a, b / 2.0]
            input = numpy.asarray(input)
            order = [0] * input.ndim
            order[axis] = 2
            return ndimage.gaussian_filter(input, sigma, order,
                                           output, mode, cval)
        for type_ in self.types:
            array = numpy.array([[3, 2, 5, 1, 4],
                                 [5, 8, 3, 7, 1],
                                 [5, 6, 9, 3, 5]], type_)
            output = numpy.zeros(array.shape, type_)
            tmp = ndimage.generic_laplace(array, derivative2,
                                          extra_arguments=(1.0,),
                                          extra_keywords={'b': 2.0})
            ndimage.gaussian_laplace(array, 1.0, output)
            assert_array_almost_equal(tmp, output) 

def predict():
    NumSample_test=2000

    names_test = loadNameList(MAT_TEST,NumSample_test,'testing')
    X_test=newload2('data_iSun_Test.cPickle')
    y_pred_test = net2.predict(X_test)

    print("X.shape == {}; X.min == {:.3f}; X.max == {:.3f}".format(
    X_test.shape, X_test.min(), X_test.max()))
    print("y.shape == {}; y.min == {:.3f}; y.max == {:.3f}".format(
    y_pred_test.shape, y_pred_test.min(), y_pred_test.max()))
      
    NumSample_val=926

    names_val = loadNameList(MAT_VAL,NumSample_val,'validation')
    X_val=newload('data_iSun_VAL.cPickle')
    y_pred_val = net2.predict(X_val)

    print("X.shape == {}; X.min == {:.3f}; X.max == {:.3f}".format(
    X_val.shape, X_val.min(), X_val.max()))
    print("y.shape == {}; y.min == {:.3f}; y.max == {:.3f}".format(
    y_pred_val.shape, y_pred_val.min(), y_pred_val.max()))

    for i in range(NumSample_val):
        img = misc.imread('/imatge/jpan/work/iSUN/images/'+names_val[i]+'.jpg')
        tmp = y_pred_val[i].reshape(48,48)
        blured= ndimage.gaussian_filter(tmp, sigma=3)
        y = misc.imresize(blured,(img.shape[0],img.shape[1]))/255.
        d = {'I':y} 
        scipy.io.savemat('/imatge/jpan/work/backup/val_isun/'+names_val[i],d)
  
    for i in range(NumSample_test):
        img = misc.imread('/imatge/jpan/work/iSUN/images/'+names_test[i]+'.jpg')
        #imsave('validation_test/'+names[i]+'.png', y)
        tmp = y_pred_test[i].reshape(48,48)
        blured= ndimage.gaussian_filter(tmp, sigma=3)
        y = misc.imresize(blured,(img.shape[0],img.shape[1]))/255.
        d = {'I':y} 
        scipy.io.savemat('/imatge/jpan/work/backup/result_isun/'+names_test[i],d) 

def saliencyPredictor(xt):
    y_pred_test = net2.predict(xt)
    i=0
    for file in glob.glob(url+"*.jpg"):   
      img = misc.imread(file)  
      tmp = y_pred_test[i].reshape(48,48)
      blured= ndimage.gaussian_filter(tmp, sigma=3)
      y = misc.imresize(blured,(img.shape[0],img.shape[1]))/255.
      misc.imsave('./saliency/'+os.path.basename(file), y)  
      i= i+1 

def _augment_images(self, images, random_state, parents, hooks):
        result = images
        nb_images = len(images)
        samples = self.sigma.draw_samples((nb_images,), random_state=random_state)
        for i in range(nb_images):
            nb_channels = images[i].shape[2]
            sig = samples[i]
            if sig > 0 + self.eps:
                # note that while gaussian_filter can be applied to all channels
                # at the same time, that should not be done here, because then
                # the blurring would also happen across channels (e.g. red
                # values might be mixed with blue values in RGB)
                for channel in range(nb_channels):
                    result[i][:, :, channel] = ndimage.gaussian_filter(result[i][:, :, channel], sig)
        return result 

def image_blurred(image, sig):
    return ndimage.gaussian_filter(image, sigma=sig)
    


# In[22]:

#Apply the transformations on all the images in a particular folder. New images will be created in aug_images sub-folder 

def im_lcn_bias(img, sigma_mean, sigma_std, std_bias):
    """
    LCN with an std bias to avoid noise amplification
    """
    means = ndimage.gaussian_filter(img, sigma_mean)
    img_centered = img - means
    stds = np.sqrt(ndimage.gaussian_filter(img_centered**2, sigma_std) + std_bias)
    return img_centered / stds 

def im_lcn(img, sigma_mean, sigma_std):
    """
    based on matlab code by Guanglei Xiong, see http://www.mathworks.com/matlabcentral/fileexchange/8303-local-normalization
    """
    means = ndimage.gaussian_filter(img, sigma_mean)
    img_centered = img - means
    stds = np.sqrt(ndimage.gaussian_filter(img_centered**2, sigma_std))
    return img_centered / stds 

def smooth_image(data, sigma):
    from scipy import ndimage
    return ndimage.gaussian_filter(data, sigma=sigma) 

def lowpass_gaussian(image, sigma=2.0):
    return nd.gaussian_filter(image, sigma=sigma) 

def highpass_gaussian(image, sigma=2.0):
    return image - nd.gaussian_filter(image, sigma=sigma) 

def _compute_gaussian_blur_ksize(sigma):
    if sigma < 3.0:
        ksize = 3.3 * sigma  # 99% of weight
    elif sigma < 5.0:
        ksize = 2.9 * sigma  # 97% of weight
    else:
        ksize = 2.6 * sigma  # 95% of weight

    # we use 5x5 here as the minimum size as that simplifies
    # comparisons with gaussian_filter() in the tests
    # TODO reduce this to 3x3
    ksize = int(max(ksize, 5))
    if ksize % 2 == 0:
        ksize += 1
    return ksize 

def test_gh_5430():
    # At least one of these raises an error unless gh-5430 is
    # fixed. In py2k an int is implemented using a C long, so
    # which one fails depends on your system. In py3k there is only
    # one arbitrary precision integer type, so both should fail.
    sigma = np.int32(1)
    out = sndi._ni_support._normalize_sequence(sigma, 1)
    assert_equal(out, [sigma])
    sigma = np.int64(1)
    out = sndi._ni_support._normalize_sequence(sigma, 1)
    assert_equal(out, [sigma])
    # This worked before; make sure it still works
    sigma = 1
    out = sndi._ni_support._normalize_sequence(sigma, 1)
    assert_equal(out, [sigma])
    # This worked before; make sure it still works
    sigma = [1, 1]
    out = sndi._ni_support._normalize_sequence(sigma, 2)
    assert_equal(out, sigma)
    # Also include the OPs original example to make sure we fixed the issue
    x = np.random.normal(size=(256, 256))
    perlin = np.zeros_like(x)
    for i in 2**np.arange(6):
        perlin += sndi.filters.gaussian_filter(x, i, mode="wrap") * i**2
    # This also fixes gh-4106, show that the OPs example now runs.
    x = np.int64(21)
    sndi._ni_support._normalize_sequence(x, 0) 

def test_gaussian_truncate():
    # Test that Gaussian filters can be truncated at different widths.
    # These tests only check that the result has the expected number
    # of nonzero elements.
    arr = np.zeros((100, 100), float)
    arr[50, 50] = 1
    num_nonzeros_2 = (sndi.gaussian_filter(arr, 5, truncate=2) > 0).sum()
    assert_equal(num_nonzeros_2, 21**2)
    num_nonzeros_5 = (sndi.gaussian_filter(arr, 5, truncate=5) > 0).sum()
    assert_equal(num_nonzeros_5, 51**2)

    # Test truncate when sigma is a sequence.
    f = sndi.gaussian_filter(arr, [0.5, 2.5], truncate=3.5)
    fpos = f > 0
    n0 = fpos.any(axis=0).sum()
    # n0 should be 2*int(2.5*3.5 + 0.5) + 1
    assert_equal(n0, 19)
    n1 = fpos.any(axis=1).sum()
    # n1 should be 2*int(0.5*3.5 + 0.5) + 1
    assert_equal(n1, 5)

    # Test gaussian_filter1d.
    x = np.zeros(51)
    x[25] = 1
    f = sndi.gaussian_filter1d(x, sigma=2, truncate=3.5)
    n = (f > 0).sum()
    assert_equal(n, 15)

    # Test gaussian_laplace
    y = sndi.gaussian_laplace(x, sigma=2, truncate=3.5)
    nonzero_indices = np.where(y != 0)[0]
    n = nonzero_indices.ptp() + 1
    assert_equal(n, 15)

    # Test gaussian_gradient_magnitude
    y = sndi.gaussian_gradient_magnitude(x, sigma=2, truncate=3.5)
    nonzero_indices = np.where(y != 0)[0]
    n = nonzero_indices.ptp() + 1
    assert_equal(n, 15) 

def test_multiple_modes_sequentially():
    # Test that the filters with multiple mode cababilities for different
    # dimensions give the same result as applying the filters with
    # different modes sequentially
    arr = np.array([[1., 0., 0.],
                    [1., 1., 0.],
                    [0., 0., 0.]])

    modes = ['reflect', 'wrap']

    expected = sndi.gaussian_filter1d(arr, 1, axis=0, mode=modes[0])
    expected = sndi.gaussian_filter1d(expected, 1, axis=1, mode=modes[1])
    assert_equal(expected,
                 sndi.gaussian_filter(arr, 1, mode=modes))

    expected = sndi.uniform_filter1d(arr, 5, axis=0, mode=modes[0])
    expected = sndi.uniform_filter1d(expected, 5, axis=1, mode=modes[1])
    assert_equal(expected,
                 sndi.uniform_filter(arr, 5, mode=modes))

    expected = sndi.maximum_filter1d(arr, size=5, axis=0, mode=modes[0])
    expected = sndi.maximum_filter1d(expected, size=5, axis=1, mode=modes[1])
    assert_equal(expected,
                 sndi.maximum_filter(arr, size=5, mode=modes))

    expected = sndi.minimum_filter1d(arr, size=5, axis=0, mode=modes[0])
    expected = sndi.minimum_filter1d(expected, size=5, axis=1, mode=modes[1])
    assert_equal(expected,
                 sndi.minimum_filter(arr, size=5, mode=modes)) 

def test_multiple_modes():
    # Test that the filters with multiple mode cababilities for different
    # dimensions give the same result as applying a single mode.
    arr = np.array([[1., 0., 0.],
                    [1., 1., 0.],
                    [0., 0., 0.]])

    mode1 = 'reflect'
    mode2 = ['reflect', 'reflect']

    assert_equal(sndi.gaussian_filter(arr, 1, mode=mode1),
                 sndi.gaussian_filter(arr, 1, mode=mode2))
    assert_equal(sndi.prewitt(arr, mode=mode1),
                 sndi.prewitt(arr, mode=mode2))
    assert_equal(sndi.sobel(arr, mode=mode1),
                 sndi.sobel(arr, mode=mode2))
    assert_equal(sndi.laplace(arr, mode=mode1),
                 sndi.laplace(arr, mode=mode2))
    assert_equal(sndi.gaussian_laplace(arr, 1, mode=mode1),
                 sndi.gaussian_laplace(arr, 1, mode=mode2))
    assert_equal(sndi.maximum_filter(arr, size=5, mode=mode1),
                 sndi.maximum_filter(arr, size=5, mode=mode2))
    assert_equal(sndi.minimum_filter(arr, size=5, mode=mode1),
                 sndi.minimum_filter(arr, size=5, mode=mode2))
    assert_equal(sndi.gaussian_gradient_magnitude(arr, 1, mode=mode1),
                 sndi.gaussian_gradient_magnitude(arr, 1, mode=mode2))
    assert_equal(sndi.uniform_filter(arr, 5, mode=mode1),
                 sndi.uniform_filter(arr, 5, mode=mode2))