def separable_convolution(input, weights, output=None, mode="reflect", cval=0.0, origin=0):
    r"""
    Calculate a n-dimensional convolution of a separable kernel to a n-dimensional input.
    
    Achieved by calling convolution1d along the first axis, obtaining an intermediate
    image, on which the next convolution1d along the second axis is called and so on.
    
    Parameters
    ----------
    input : array_like
        Array of which to estimate the noise.
    weights : ndarray
        One-dimensional sequence of numbers.          
    output : array, optional
        The `output` parameter passes an array in which to store the
        filter output.
    mode : {'reflect', 'constant', 'nearest', 'mirror', 'wrap'}, optional
        The `mode` parameter determines how the array borders are
        handled, where `cval` is the value when mode is equal to
        'constant'. Default is 'reflect'
    cval : scalar, optional
        Value to fill past edges of input if `mode` is 'constant'. Default
        is 0.0
    origin : scalar, optional
        The `origin` parameter controls the placement of the filter.
        Default 0.0.
        
    Returns
    -------
    output : ndarray
        Input image convolved with the supplied kernel.
    """
    input = numpy.asarray(input)
    output = _ni_support._get_output(output, input)
    axes = list(range(input.ndim))
    if len(axes) > 0:
        convolve1d(input, weights, axes[0], output, mode, cval, origin)
        for ii in range(1, len(axes)):
            convolve1d(output, weights, axes[ii], output, mode, cval, origin)
    else:
        output[...] = input[...]
    return output