Python scipy.signal.gaussian() Examples

def preprocess(data, smoothing_sd=25, n_pcs=20):
    """
    Preprocess neural data for cca analysis with smoothing and pca

    :param data: array of shape (n_samples, n_features)
    :type data: array-like
    :param smoothing_sd: gaussian smoothing kernel standard deviation (ms)
    :type smoothing_sd: float
    :param n_pcs: number of pca dimensions to retain
    :type n_pcs: int
    :return: preprocessed neural data
    :rtype: array-like, shape (n_samples, pca_dims)
    """
    if smoothing_sd > 0:
        data = _smooth(data, sd=smoothing_sd)
    if n_pcs > 0:
        data = _pca(data, n_pcs=n_pcs)
    return data 

def smooth_dir_map(dir_map,sigma=2.0,mask = None):

    cos2Theta = np.cos(dir_map * 2)
    sin2Theta = np.sin(dir_map * 2)
    if mask is not None:
        assert (dir_map.shape[0] == mask.shape[0])
        assert (dir_map.shape[1] == mask.shape[1])
        cos2Theta[mask == 0] = 0
        sin2Theta[mask == 0] = 0

    cos2Theta = gaussian(cos2Theta, sigma, multichannel=False, mode='reflect')
    sin2Theta = gaussian(sin2Theta, sigma, multichannel=False, mode='reflect')

    dir_map = np.arctan2(sin2Theta,cos2Theta)*0.5


    return dir_map 

def test_smooth_1d():
    for edge in ['m', 'c']:
        for N in [20,21]:
            # values in [9.0,11.0]
            x = rand(N) + 10
            mn = 9.0
            mx = 11.0
            for M in range(18,27):
                print("1d", edge, "N=%i, M=%i" %(N,M))
                xsm = smooth(x, gaussian(M,2.0), edge=edge)
                assert len(xsm) == N
                # (N,1) case
                xsm2 = smooth(x[:,None], gaussian(M,2.0)[:,None], edge=edge)
                assert np.allclose(xsm, xsm2[:,0], atol=1e-14, rtol=1e-12)
                # Smoothed signal should not go to zero if edge effects are handled
                # properly. Also assert proper normalization (i.e. smoothed signal
                # is "in the middle" of the noisy original data).
                assert xsm.min() >= mn
                assert xsm.max() <= mx
                assert mn <= xsm[0] <= mx
                assert mn <= xsm[-1] <= mx
            # convolution with delta peak produces same data exactly
            assert np.allclose(smooth(x, np.array([0.0,1,0]), edge=edge),x, atol=1e-14,
                               rtol=1e-12) 

def test_smooth_nd():
    for edge in ['m', 'c']:
        a = rand(20, 2, 3) + 10
        for M in [5, 20, 123]:
            print("nd", edge, "M=%i" %M)
            kern = gaussian(M, 2.0)
            asm = smooth(a, kern[:,None,None], axis=0, edge=edge)
            assert asm.shape == a.shape
            for jj in range(asm.shape[1]):
                for kk in range(asm.shape[2]):
                    assert np.allclose(asm[:,jj,kk], smooth(a[:,jj,kk], kern,
                                                            edge=edge))
                    mn = a[:,jj,kk].min()
                    mx = a[:,jj,kk].max()
                    smn = asm[:,jj,kk].min()
                    smx = asm[:,jj,kk].max()
                    assert smn >= mn, "min: data=%f, smooth=%f" %(mn, smn)
                    assert smx <= mx, "max: data=%f, smooth=%f" %(mx, smx) 

def _gen_gaussian_kernel(self):
        sigma = self.gaussian_sigma
        kernel = th.zeros(self.kernel_size)
        delta_theta = self.kernel_rad / (self.kernel_size[0] - 1)
        sigma_idx = sigma / delta_theta
        gauss1d = signal.gaussian(2 * kernel.shape[0], sigma_idx)
        gauss2d = np.outer(gauss1d, np.ones(kernel.shape[1]))

        return gauss2d[-kernel.shape[0]:, :] 

def GaussTemporalFilter(x, order=11):
    """ x: energy-normalized spectrum """

    window = signal.gaussian(order, std=order/5)
    window = window / np.sum(window)
    
    _x = np.zeros_like(x)
    for i in range(513):
        _x[:, i] = np.convolve(np.log10(x[:, i]), window, 'same')

    _x = np.power(10., _x)
    _x = _x / np.sum(_x, 1).reshape([-1, 1])

    return _x 

def _gen_gaussian_kernel(self):
        sigma = self.gaussian_sigma
        kernel = th.zeros(self.kernel_size)
        delta_theta = self.kernel_rad / (self.kernel_size[0] - 1)
        sigma_idx = sigma / delta_theta
        gauss1d = signal.gaussian(2 * kernel.shape[0], sigma_idx)
        gauss2d = np.outer(gauss1d, np.ones(kernel.shape[1]))

        return gauss2d[-kernel.shape[0]:, :] 

def _gen_gaussian_kernel(self):
        sigma = self.gaussian_sigma
        kernel = th.zeros(self.kernel_size)
        delta_theta = self.kernel_rad / (self.kernel_size[0] - 1)
        sigma_idx = sigma / delta_theta
        gauss1d = signal.gaussian(2 * kernel.shape[0], sigma_idx)
        gauss2d = np.outer(gauss1d, np.ones(kernel.shape[1]))

        return gauss2d[-kernel.shape[0]:, :] 

def _gen_gaussian_kernel(self):
        sigma = self.gaussian_sigma
        kernel = th.zeros(self.kernel_size)
        delta_theta = self.kernel_rad / (self.kernel_size[0] - 1)
        sigma_idx = sigma / delta_theta
        gauss1d = signal.gaussian(2 * kernel.shape[0], sigma_idx)
        gauss2d = np.outer(gauss1d, np.ones(kernel.shape[1]))

        return gauss2d[-kernel.shape[0]:, :] 

def _gen_gaussian_kernel(self):
        sigma = self.gaussian_sigma
        kernel = th.zeros(self.kernel_size)
        delta_theta = self.kernel_rad / (self.kernel_size[0] - 1)
        sigma_idx = sigma / delta_theta
        gauss1d = signal.gaussian(2 * kernel.shape[0], sigma_idx)
        gauss2d = np.outer(gauss1d, np.ones(kernel.shape[1]))

        return gauss2d[-kernel.shape[0]:, :] 

def gen_gaussian_kernel(sigma_idx=8, kernel_size=(15, 30)):
        gauss1d = signal.gaussian(2 * kernel_size[0], sigma_idx)
        gauss2d = np.outer(gauss1d, np.ones(kernel_size[1]))

        return gauss2d[-kernel_size[0]:, :] 

def render(
    locs,
    info=None,
    oversampling=1,
    viewport=None,
    blur_method=None,
    min_blur_width=0,
):
    if viewport is None:
        try:
            viewport = [(0, 0), (info[0]["Height"], info[0]["Width"])]
        except TypeError:
            raise ValueError("Need info if no viewport is provided.")
    (y_min, x_min), (y_max, x_max) = viewport
    if blur_method is None:
        return render_hist(locs, oversampling, y_min, x_min, y_max, x_max)
    elif blur_method == "gaussian":
        return render_gaussian(
            locs, oversampling, y_min, x_min, y_max, x_max, min_blur_width
        )
    elif blur_method == "gaussian_iso":
        return render_gaussian_iso(
            locs, oversampling, y_min, x_min, y_max, x_max, min_blur_width
        )
    elif blur_method == "smooth":
        return render_smooth(locs, oversampling, y_min, x_min, y_max, x_max)
    elif blur_method == "convolve":
        return render_convolve(
            locs, oversampling, y_min, x_min, y_max, x_max, min_blur_width
        )
    else:
        raise Exception("blur_method not understood.") 

def _fftconvolve(image, blur_width, blur_height):
    kernel_width = 10 * int(_np.round(blur_width)) + 1
    kernel_height = 10 * int(_np.round(blur_height)) + 1
    kernel_y = _signal.gaussian(kernel_height, blur_height)
    kernel_x = _signal.gaussian(kernel_width, blur_width)
    kernel = _np.outer(kernel_y, kernel_x)
    kernel /= kernel.sum()
    return _signal.fftconvolve(image, kernel, mode="same") 

def _gkern2(kernlen=21, nsig=3):
    """Returns a 2D Gaussian kernel array."""

    # create nxn zeros
    inp = np.zeros((kernlen, kernlen))
    # set element at the middle to one, a dirac delta
    inp[kernlen // 2, kernlen // 2] = 1
    # gaussian-smooth the dirac, resulting in a gaussian filter mask
    return fi.gaussian_filter(inp, nsig) 

def timefilter(
    xr_in, steps, step_spec, timename="time", filtertype="gaussian", stdev=0.1
):
    timedim = xr_in.dims.index(timename)
    dt = np.diff(xr_in.time.data[0:2])[0]
    cut_dt = np.timedelta64(steps, step_spec)

    if filtertype == "gaussian":
        win_length = (cut_dt / dt).astype(int)
        a = [1.0]
        win = gaussian(win_length, std=(float(win_length) * stdev))
        b = win / win.sum()
        if np.nansum(win) == 0:
            raise RuntimeError("window to short for time interval")
            print("win_length", str(win_length))
            print("stddev", str(stdev))
            print("win", str(win))

    filtered = filtfilt(b, a, xr_in.data, axis=timedim, padtype=None, padlen=0)
    out = xr.DataArray(
        filtered, dims=xr_in.dims, coords=xr_in.coords, attrs=xr_in.attrs
    )
    out.attrs.update({"filterlength": (steps, step_spec), "filtertype": filtertype})
    if xr_in.name:
        out.name = xr_in.name + "_lowpassed"
    return out 

def gaussian_2dkernel(size=5, sig=1.):
    """
    Returns a 2D Gaussian kernel array with side length size and a sigma of sig
    """
    gkern1d = signal.gaussian(size, std=sig).reshape(size, 1)
    gkern2d = np.outer(gkern1d, gkern1d)
    return (gkern2d/gkern2d.sum()) 

def smooth_class_prior(sigma=5, do_plot=False):
    prior_prob = np.load(os.path.join(data_dir, "prior_prob.npy"))
    # add an epsilon to prior prob to avoid 0 vakues and possible NaN
    prior_prob += 1E-3 * np.min(prior_prob)
    # renormalize
    prior_prob = prior_prob / (1.0 * np.sum(prior_prob))

    # Smooth with gaussian
    f = interp1d(np.arange(prior_prob.shape[0]), prior_prob)
    xx = np.linspace(0, prior_prob.shape[0] - 1, 1000)
    yy = f(xx)
    window = gaussian(2000, sigma)  # 2000 pts in the window, sigma=5
    smoothed = convolve(yy, window / window.sum(), mode='same')
    fout = interp1d(xx, smoothed)
    prior_prob_smoothed = np.array([fout(i) for i in range(prior_prob.shape[0])])
    prior_prob_smoothed = prior_prob_smoothed / np.sum(prior_prob_smoothed)

    # Save
    file_name = os.path.join(data_dir, "prior_prob_smoothed.npy")
    np.save(file_name, prior_prob_smoothed)

    if do_plot:
        plt.plot(prior_prob)
        plt.plot(prior_prob_smoothed, "g--")
        plt.plot(xx, smoothed, "r-")
        plt.yscale("log")
        plt.show() 

def smooth_color_prior(size=64, sigma=5, do_plot=False):

    prior_prob = np.load(os.path.join(data_dir, "CelebA_%s_prior_prob.npy" % size))
    # add an epsilon to prior prob to avoid 0 vakues and possible NaN
    prior_prob += 1E-3 * np.min(prior_prob)
    # renormalize
    prior_prob = prior_prob / (1.0 * np.sum(prior_prob))

    # Smooth with gaussian
    f = interp1d(np.arange(prior_prob.shape[0]),prior_prob)
    xx = np.linspace(0,prior_prob.shape[0] - 1, 1000)
    yy = f(xx)
    window = gaussian(2000, sigma)  # 2000 pts in the window, sigma=5
    smoothed = convolve(yy, window / window.sum(), mode='same')
    fout = interp1d(xx,smoothed)
    prior_prob_smoothed = np.array([fout(i) for i in range(prior_prob.shape[0])])
    prior_prob_smoothed = prior_prob_smoothed / np.sum(prior_prob_smoothed)

    # Save
    file_name = os.path.join(data_dir, "CelebA_%s_prior_prob_smoothed.npy" % size)
    np.save(file_name, prior_prob_smoothed)

    if do_plot:
        plt.plot(prior_prob)
        plt.plot(prior_prob_smoothed, "g--")
        plt.plot(xx, smoothed, "r-")
        plt.yscale("log")
        plt.show() 

def orientation(image, stride=8, window=17):
    with K.tf.name_scope('orientation'):
        assert image.get_shape().as_list()[3] == 1, 'Images must be grayscale'
        strides = [1, stride, stride, 1]
        E = np.ones([window, window, 1, 1])
        sobelx = np.reshape(np.array([[-1, 0, 1], [-2, 0, 2], [-1, 0, 1]], dtype=float), [3, 3, 1, 1])
        sobely = np.reshape(np.array([[-1, -2, -1], [0, 0, 0], [1, 2, 1]], dtype=float), [3, 3, 1, 1])
        gaussian = np.reshape(gaussian2d((5, 5), 1), [5, 5, 1, 1])
        with K.tf.name_scope('sobel_gradient'):
            Ix = K.tf.nn.conv2d(image, sobelx, strides=[1,1,1,1], padding='SAME', name='sobel_x')
            Iy = K.tf.nn.conv2d(image, sobely, strides=[1,1,1,1], padding='SAME', name='sobel_y')
        with K.tf.name_scope('eltwise_1'):
            Ix2 = K.tf.multiply(Ix, Ix, name='IxIx')
            Iy2 = K.tf.multiply(Iy, Iy, name='IyIy')
            Ixy = K.tf.multiply(Ix, Iy, name='IxIy')
        with K.tf.name_scope('range_sum'):
            Gxx = K.tf.nn.conv2d(Ix2, E, strides=strides, padding='SAME', name='Gxx_sum')
            Gyy = K.tf.nn.conv2d(Iy2, E, strides=strides, padding='SAME', name='Gyy_sum')
            Gxy = K.tf.nn.conv2d(Ixy, E, strides=strides, padding='SAME', name='Gxy_sum')
        with K.tf.name_scope('eltwise_2'):
            Gxx_Gyy = K.tf.subtract(Gxx, Gyy, name='Gxx_Gyy')
            theta = atan2([2*Gxy, Gxx_Gyy]) + np.pi
        # two-dimensional low-pass filter: Gaussian filter here
        with K.tf.name_scope('gaussian_filter'):
            phi_x = K.tf.nn.conv2d(K.tf.cos(theta), gaussian, strides=[1,1,1,1], padding='SAME', name='gaussian_x')
            phi_y = K.tf.nn.conv2d(K.tf.sin(theta), gaussian, strides=[1,1,1,1], padding='SAME', name='gaussian_y')
            theta = atan2([phi_y, phi_x])/2
    return theta 

def gaussian2d(shape=(5,5),sigma=0.5):
    """
    2D gaussian mask - should give the same result as MATLAB's
    fspecial('gaussian',[shape],[sigma])
    """
    m,n = [(ss-1.)/2. for ss in shape]
    y,x = np.ogrid[-m:m+1,-n:n+1]
    h = np.exp( -(x*x + y*y) / (2.*sigma*sigma) )
    h[ h < np.finfo(h.dtype).eps*h.max() ] = 0
    sumh = h.sum()
    if sumh != 0:
        h /= sumh
    return h 

def gausslabel(length=180, stride=2):
    gaussian_pdf = signal.gaussian(length+1, 3)
    label = np.reshape(np.arange(stride/2, length, stride), [1,1,-1,1])
    y = np.reshape(np.arange(stride/2, length, stride), [1,1,1,-1])
    delta = np.array(np.abs(label - y), dtype=int)
    delta = np.minimum(delta, length-delta)+length/2
    return gaussian_pdf[delta] 

def get_quality_map_ori_dict(img, dict, spacing, dir_map = None, block_size = 16):
    if img.dtype=='uint8':
        img = img.astype(np.float)
    img = FastEnhanceTexture(img)
    h, w = img.shape
    blkH, blkW = dir_map.shape

    quality_map = np.zeros((blkH,blkW),dtype=np.float)
    fre_map = np.zeros((blkH,blkW),dtype=np.float)
    ori_num = len(dict)
    #dir_map = math.pi/2 - dir_map
    dir_ind = dir_map*ori_num/math.pi
    dir_ind = dir_ind.astype(np.int)
    dir_ind = dir_ind%ori_num

    patch_size = np.sqrt(dict[0].shape[1])
    patch_size = patch_size.astype(np.int)
    pad_size = (patch_size-block_size)//2
    img = np.lib.pad(img, (pad_size, pad_size), 'symmetric')
    for i in range(0,blkH):
        for j in range(0,blkW):
            ind = dir_ind[i,j]
            patch = img[i*block_size:i*block_size+patch_size,j*block_size:j*block_size+patch_size]

            patch = patch.reshape(patch_size*patch_size,)
            patch = patch - np.mean(patch)
            patch = patch / (np.linalg.norm(patch)+0.0001)
            patch[patch>0.05] = 0.05
            patch[patch<-0.05] = -0.05

            simi = np.dot(dict[ind], patch)
            similar_ind = np.argmax(abs(simi))
            quality_map[i,j] = np.max(abs(simi))
            fre_map[i,j] = 1./spacing[ind][similar_ind]

    quality_map = gaussian(quality_map,sigma=2)
    return quality_map, fre_map 

def _smooth(data, sd):
    from scipy.signal import gaussian
    from scipy.signal import convolve
    n_bins = data.shape[0]
    w = n_bins - 1 if n_bins % 2 == 0 else n_bins
    window = gaussian(w, std=sd)
    for j in range(data.shape[1]):
        data[:, j] = convolve(data[:, j], window, mode='same', method='auto')
    return data 

def velocity_smoothed(pos, freq, smooth_size=0.03):
    """
    Compute wheel velocity from uniformly sampled wheel data

    Parameters
    ----------
    pos : array_like
        Array of wheel positions
    smooth_size : float
        Size of Gaussian smoothing window in seconds
    freq : float
        Sampling frequency of the data

    Returns
    -------
    vel : np.ndarray
        Array of velocity values
    acc : np.ndarray
        Array of acceleration values
    """
    # Define our smoothing window with an area of 1 so the units won't be changed
    stdSamps = np.round(smooth_size * freq)  # Standard deviation relative to sampling frequency
    N = stdSamps * 6  # Number of points in the Gaussian
    gauss_std = (N - 1) / 6  # @fixme magic number everywhere!
    win = gaussian(N, gauss_std)
    win = win / win.sum()  # Normalize amplitude

    # Convolve and multiply by sampling frequency to restore original units
    vel = np.insert(convolve(np.diff(pos), win, mode='same'), 0, 0) * freq
    acc = np.insert(convolve(np.diff(vel), win, mode='same'), 0, 0) * freq

    return vel, acc 

def __init__(
            self,
            action_space=(0, 1, 2, 3),
            time_threshold=5,
            pips_threshold=10,
            pips_scale=1e-4,
            kernel_size=5,
            kernel_stddev=1
    ):
        """

        Args:
            action_space:       actions to advice: 0 - hold, 1- buy, 2 - sell, 3 - close
            time_threshold:     how many points (in number of ENVIRONMENT timesteps) on each side to use
                                for the comparison to consider comparator(n, n+x) to be True
            pips_threshold:     int, minimal peaks difference in pips
                                to consider comparator(n, n+x) to be True
            pips_scale:         actual single pip value wrt signal value
            kernel_size:        gaussian convolution kernel size (used to compute distribution over actions)
            kernel_stddev:      gaussian kernel standard deviation
        """
        self.action_space = action_space
        self.time_threshold = time_threshold
        self.value_threshold = pips_threshold * pips_scale
        self.kernel_size = kernel_size
        self.kernel = signal.gaussian(kernel_size, std=kernel_stddev)
        self.data = None 

def fit(self, episode_data, resampling_factor=1):
        """
        Estimates `advised` actions probabilities distribution based on data received.

        Args:
            episode_data:           1D np.array of unscaled price values in OHL[CV] format
            resampling_factor:      factor by which to resample given data
                                    by taking min/max values inside every resampled bar

        Returns:
             Np.array of size [resampled_data_size, actions_space_size] of probabilities of advised actions, where
             resampled_data_size = int(len(episode_data) / resampling_factor) + 1/0

        """
        # Vector of advised actions:
        data = self.resample_data(episode_data, resampling_factor)
        signals = self.estimate_actions(data)
        signals = self.adjust_signals(signals)

        # One-hot actions encoding:
        actions_one_hot = np.zeros([signals.shape[0], len(self.action_space)])
        actions_one_hot[np.arange(signals.shape[0]), signals] = 1

        # Want a bit relaxed discrete distribution over actions instead of one hot (heuristic):
        actions_distr = np.zeros(actions_one_hot.shape)

        # For all actions except 'hold' (due to heuristic skewness):
        actions_distr[:, 0] = actions_one_hot[:, 0]

        # ...spread out actions probabilities by convolving with gaussian kernel :
        for channel in range(1, actions_one_hot.shape[-1]):
            actions_distr[:, channel] = np.convolve(actions_one_hot[:, channel], self.kernel, mode='same') + 0.1

        # Normalize:
        actions_distr /= actions_distr.sum(axis=-1)[..., None]

        return actions_distr 

def lorentz(M, std=1.0, sym=True):
    r"""Lorentz window (same as Cauchy function). Function skeleton stolen from
    scipy.signal.gaussian().

    The Lorentz function is

    .. math::

        L(x) = \frac{\Gamma}{(x-x_0)^2 + \Gamma^2}

    Here :math:`x_0 = 0` and `std` = :math:`\Gamma`.
    Some definitions use :math:`1/2\,\Gamma` instead of :math:`\Gamma`, but
    without 1/2 we get comparable peak width to Gaussians when using this
    window in convolutions, thus ``scipy.signal.gaussian(M, std=5)`` is similar
    to ``lorentz(M, std=5)``.

    Parameters
    ----------
    M : int
        number of points
    std : float
        spread parameter :math:`\Gamma`
    sym : bool

    Returns
    -------
    w : (M,)
    """
    if M < 1:
        return np.array([])
    if M == 1:
        return np.ones(1,dtype=float)
    odd = M % 2
    if not sym and not odd:
        M = M+1
    n = np.arange(0, M) - (M - 1.0) / 2.0
    w = std / (n**2.0 + std**2.0)
    w /= w.max()
    if not sym and not odd:
        w = w[:-1]
    return w 

def smooth1d(array, window_size=None, kernel='gaussian'):
    """Apply a centered window smoothing to a 1D array.

    Parameters
    ----------
    array : ndarray
        the array to apply the smoothing to
    window_size : int
        the size of the smoothing window
    kernel : str
        the type of smoothing (`gaussian`, `mean`)

    Returns
    -------
    the smoothed array (same dim as input)
    """

    # some defaults
    if window_size is None:
        if len(array) >= 9:
            window_size = 9
        elif len(array) >= 7:
            window_size = 7
        elif len(array) >= 5:
            window_size = 5
        elif len(array) >= 3:
            window_size = 3

    if window_size % 2 == 0:
        raise ValueError('Window should be an odd number.')

    if isinstance(kernel, str):
        if kernel == 'gaussian':
            kernel = gaussian(window_size, 1)
        elif kernel == 'mean':
            kernel = np.ones(window_size)
        else:
            raise NotImplementedError('Kernel: ' + kernel)
    kernel = kernel / np.asarray(kernel).sum()
    return filters.convolve1d(array, kernel, mode='mirror') 

def _generate_noise_temporal_task(stimfunction_tr,
                                  motion_noise='gaussian',
                                  ):
    """Generate the signal dependent noise

    Create noise specific to the signal, for instance there is variability
    in how the signal manifests on each event

    Parameters
    ----------

    stimfunction_tr : 1 Dimensional array
        This is the timecourse of the stimuli in this experiment,
        each element represents a TR

    motion_noise : str
        What type of noise will you generate? Can be gaussian or rician

    Returns
    ----------

    noise_task : one dimensional array, float
        Generates the temporal task noise timecourse

    """

    # Make the noise to be added
    stimfunction_tr = stimfunction_tr != 0
    if motion_noise == 'gaussian':
        noise = stimfunction_tr * np.random.normal(0, 1,
                                                   size=stimfunction_tr.shape)
    elif motion_noise == 'rician':
        noise = stimfunction_tr * stats.rice.rvs(0, 1,
                                                 size=stimfunction_tr.shape)

    noise_task = stimfunction_tr + noise

    # Normalize
    noise_task = stats.zscore(noise_task).flatten()

    return noise_task 

def isogkern(kernlen, std):
    gkern1d = signal.gaussian(kernlen, std=std).reshape(kernlen, 1)
    gkern2d = np.outer(gkern1d, gkern1d)
    gkern2d = gkern2d/np.sum(gkern2d)
    return gkern2d