def immerkaer_local(input, size, output=None, mode="reflect", cval=0.0):
    r"""
    Estimate the local noise.
    
    The input image is assumed to have additive zero mean Gaussian noise. The Immerkaer
    noise estimation is applied to the image locally over a N-dimensional cube of
    side-length size. The size of the region should be sufficiently high for a stable
    noise estimation.
    
    Parameters
    ----------
    input : array_like
        Array of which to estimate the noise.
    size : integer
        The local region's side length.
    output : ndarray, 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        
        
    Returns
    -------
    sigmas : array_like
        Map of the estimated standard deviation of the images Gaussian noise per voxel.
    
    Notes
    -----
    Does not take the voxel spacing into account.
    Works good with medium to strong noise. Tends to underestimate for low noise levels.
    
    See also
    --------
    immerkaer
    """
    output = _ni_support._get_output(output, input)
    footprint = numpy.asarray([1] * size)
    
    # build nd-kernel to acquire square root of sum of squared elements
    kernel = [1, -2, 1]
    for _ in range(input.ndim - 1):
        kernel = numpy.tensordot(kernel, [1, -2, 1], 0)
    divider = numpy.square(numpy.abs(kernel)).sum() # 36 for 1d, 216 for 3D, etc.
    
    # compute laplace of input
    laplace = separable_convolution(input, [1, -2, 1], numpy.double, mode, cval)
    
    # compute factor
    factor = numpy.sqrt(numpy.pi / 2.) * 1. / ( numpy.sqrt(divider) * numpy.power(footprint.size, laplace.ndim) )
    
    # locally sum laplacian values
    separable_convolution(numpy.abs(laplace), footprint, output, mode, cval)
    
    output *= factor
    
    return output