docs: 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