def noise_power_spectrum(data, ground_truth, radial=False,
                         radial_binning_factor=2.0):
    """Return the Noise Power Spectrum (NPS).

    The NPS is given by the squared magnitude of the Fourier transform of the
    noise.

    Parameters
    ----------
    data : `DiscreteLpElement` or `array-like`
        Input data to compare to the ground truth. If not a
        `DiscreteLpElement`, a default space with cell size 1 will be assumed.
    ground_truth : `array-like`
        Reference to which ``data`` should be compared.
    radial : bool
        If ``True``, compute the radial NPS.
    radial_binning_factor : positive float, optional
        Reduce the number of radial bins by this factor. Increasing this
        number can help reducing fluctuations due to the variance of points
        that fall in a particular annulus.
        A binning factor of ``1`` corresponds to a bin size equal to
        image pixel size for images with square pixels, otherwise ::

            max(norm2(c)) / norm2(shape)

        where the maximum is taken over all corners of the image domain.

    Returns
    -------
    noise_power_spectrum : `DiscreteLp`-element
        The space is the Fourier space corresponding to ``space``, and
        hence the axes indicate frequency.
        If ``radial`` is ``True``, an average over concentric annuli is
        taken. The result is an element of a one-dimensional space with
        domain ``[0, rmax]``, where ``rmax`` is the radius of the smallest
        ball containing ``space.domain``. Its shape is ``(N,)`` with ::

            N = int(sqrt(sum(n ** 2 for n in image.shape)) / binning_factor)
    """
    try:
        space = data.space
        assert isinstance(space, odl.DiscreteLp)
    except (AttributeError, AssertionError):
        data = np.asarray(data)
        space = odl.uniform_discr(
            [0] * data.ndim, data.shape, data.shape, data.dtype
        )
        data = space.element(data)

    ft = odl.trafos.FourierTransform(space, halfcomplex=False)
    nps = np.abs(ft(data - ground_truth)).real ** 2

    if radial:
        return spherical_sum(nps, binning_factor=radial_binning_factor)
    else:
        return nps
Return the Noise Power Spectrum (NPS).

    The NPS is given by the squared magnitude of the Fourier transform of the
    noise.

    Parameters
    ----------
    data : `DiscreteLpElement` or `array-like`
        Input data to compare to the ground truth. If not a
        `DiscreteLpElement`, a default space with cell size 1 will be assumed.
    ground_truth : `array-like`
        Reference to which ``data`` should be compared.
    radial : bool
        If ``True``, compute the radial NPS.
    radial_binning_factor : positive float, optional
        Reduce the number of radial bins by this factor. Increasing this
        number can help reducing fluctuations due to the variance of points
        that fall in a particular annulus.
        A binning factor of ``1`` corresponds to a bin size equal to
        image pixel size for images with square pixels, otherwise ::

            max(norm2(c)) / norm2(shape)

        where the maximum is taken over all corners of the image domain.

    Returns
    -------
    noise_power_spectrum : `DiscreteLp`-element
        The space is the Fourier space corresponding to ``space``, and
        hence the axes indicate frequency.
        If ``radial`` is ``True``, an average over concentric annuli is
        taken. The result is an element of a one-dimensional space with
        domain ``[0, rmax]``, where ``rmax`` is the radius of the smallest
        ball containing ``space.domain``. Its shape is ``(N,)`` with ::

            N = int(sqrt(sum(n ** 2 for n in image.shape)) / binning_factor)