Improving convolution

[AlecThomson]

Hi all,

It's great to see that spectral cube now has a convolve_to method - this is going to be super useful. I see, though, that by default astropy convolution is used. I've encountered issues with this implementation, especially when convolving to a common resolution along a spectral axis. Both astropy and scipy (which I believe is much faster than astropy, but not robust to NaNs) require that the convolution kernel is defined in the image plane - even if using Fourier convolution. When the convolution kernel is small (e.g. when only marginally increasing the beam size), an image-plane kernel can be under-sampled. This will then wreak havoc on the convolved image.

We've been dealing with this issue in our own project. We do both a check for an undersampled beam, and provide an analytic convolution method that avoids the above problem.

I think it'd be important to add a check for an undersampled convolving beam, and at least issue a warning to the user (or even raise an error). If the maintainers are also interested, I'd be happy to help in adding our convolution method to spectral-cube as an option as well.

Let me know what you think!

[keflavich]

Yes, we'd welcome the contribution!

The astropy convolution is about as fast as the scipy methods (with a little extra overhead) if you disable the normalized convolution feature; see astropy/astropy#11242 (comment). (but scipy's real convolution can be quite a bit faster and more memory efficient)

However, the small-kernel discretization issue is an interesting one that I haven't encountered in practice - do you have any examples of when this is a problem? Is there a good heuristic for when we should be raising that warning?

[akleroy]

(sorry somewhat out of the blue)

This would be amazing!

If you are looking for a use case: we regularly convolve our ALMA or VLA data products to have a round synthesized beam before analysis. This is an ideal case for the analytic method described above because using image-based cases you have to convolve the major axis as well and pad to avoid issues as @AlecThomson describes, so you end up losing resolution for no reason other than software. The analytic case solves this nicely (and there was a prototype version in the PHANGS pipeline from Erik Rosolowsky at one point). Having a tested, production version of this would be very handy!

[low-sky]

gist FWIW
Regardless, the barrier we ran up into was the aliasing of the kernel for small (or zero) width kernels in the image plane. This is nominally fine, but the image needs to be well sampled. Otherwise, there's non-zero power out near where the kernel hits the edge of the FT domain and nominally wraps.

This would be great to get into a fully fledged and tested method. There are places where it can go wrong, but most radio applications won't encounter this.

[keflavich]

Thanks all, this is clearly a valuable addition.

@low-sky If I read your gist correctly, figure 1 shows that the FFT-based kernel is identical to an analytic kernel that is larger by (1/N), where N is the size of the kernel (the text in the gist says image, but I think it should be the kernel)? For a well-sampled beam, this is then a generally small effect, but one that is totally avoidable.

[AlecThomson]

Hi all,

Here's another little gist to demonstrate the issue. Even for a well-sampled beam, the convolving beam (specifically the image kernel) can be under-sampled. I think the most likely case for this occur is during mosicacking, where PSF rotations between fields can mean very small convolving beam. We could also imagine a case where we have narrowly-spaced channels - with a very similar PSF - that need to be convolved to a common resolution.

Here is a check that we use in RACS-tools for beam sampling - we check if the minor axis is sampled by at least 2 pixels on the grid (following the Nyquist-Shannon theorem). This is possibly a bit too harsh - since its probably BMIN*sin(BPA) that needs to be critically sampled.

w = astropy.wcs.WCS(target_header)
pixelscales = astropy.wcs.utils.proj_plane_pixel_scales(w)
dx = pixelscales[0] * u.deg
dy = pixelscales[1] * u.deg
if not dx == dy:
    raise Exception("GRID MUST BE SAME IN X AND Y")
grid = dy
if conv_mode != "robust": # i.e. not analytic
    # Get the minor axis of the convolving beams
    minorcons = []
    for beam in beams:
        minorcons += [cmn_beam.deconvolve(beam).minor.to(u.arcsec).value]
    minorcons = np.array(minorcons) * u.arcsec
    samps = minorcons / grid.to(u.arcsec)
    # Check that convolving beam will be Nyquist sampled
    if any(samps.value < 2):
        # Set the convolving beam to be Nyquist sampled
        nyq_con_beam = Beam(major=grid * 2, minor=grid * 2, pa=0 * u.deg)
        # Find new target based on common beam * Nyquist beam
        # Not sure if this is best - but it works
        nyq_beam = cmn_beam.convolve(nyq_con_beam)
        nyq_beam = Beam(
            major=my_ceil(nyq_beam.major.to(u.arcsec).value, precision=1)
            * u.arcsec,
            minor=my_ceil(nyq_beam.minor.to(u.arcsec).value, precision=1)
            * u.arcsec,
            pa=round_up(nyq_beam.pa.to(u.deg), decimals=2),
        )
        logger.info(f"Smallest common Nyquist sampled beam is: {nyq_beam!r}")
        if target_beam is not None:
            if target_beam < nyq_beam:
                logger.warning("TARGET BEAM WILL BE UNDERSAMPLED!")
                raise Exception("CAN'T UNDERSAMPLE BEAM - EXITING")
        else:
            logger.warning("COMMON BEAM WILL BE UNDERSAMPLED!")
            logger.warning("SETTING COMMON BEAM TO NYQUIST BEAM")
            cmn_beam = nyq_beam