Use Richardson-Lucy deconvolution from `scikit-image`

[wtbarnes]

Description

In the xrt_deconvolve function, a custom implementation of Richardonson-Lucy deconvolution is used. However, there is an existing implementation of this algorithm in scikit-image: https://scikit-image.org/docs/stable/api/skimage.restoration.html#skimage.restoration.richardson_lucy. We should investigate whether this implementation is suitable for our needs and if so, remove the custom implementation and rely on this instead. In this case, there is little cost in adding this package as a dependency as it is a widely used package in the scientific Python community and is also an optional dependency of sunpy already.

I've opened a similar issue on aiapy as that package also has a custom implementation of RL: https://gitlab.com/LMSAL_HUB/aia_hub/aiapy/-/issues/136.

[jslavin]

I looked at the scikit-image routine when I first started to port the IDL version of xrt_deconvolve to python. The results did not match the IDL results closely and appeared worse (noisier). It might be worth looking into it again though.

[wtbarnes]

Ah ok that is good to know. I haven't ever actually used this implementation as I don't think it existed when I implemented the version that is in aiapy.

Either way, it doesn't seem like there is a need for each package to implement version of RL deconvolution. I wonder if a more general deconvolution function that operates on Map objects could live in the sunkit-image package. This function could then have multiple different ways to perform the deconvolution and could also include methods for speeding up the deconvolution (e.g. with a GPU as is done in aiapy).

[jslavin]

I also recall that Kathy was not happy with some results she got using the aiapy deconvolution routine (on AIA data, of course). I think she said that she created an issue on it.
I seem to recall that part of the problem with the using the scikit-image routine was the definition of the point spread function which is off-centered by a half pixel in the XRT point spread images. I'd like to look into this more, but it's down the list a bit for me.

[wtbarnes]

I actually just responded with a follow-up to that issue on aiapy regarding the results of the deconvolution algorithm. I'm unable to reproduce those results currently, but I think this is even more of an argument for having one canonical implementation of the deconvolution: we can have one implementation of an algorithm (or even a set of implementations) in one package such that bugs need only be fixed in one place. Even if we don't end up using what is in scikit-image, we could at least have an implementation that works for solar images given some PSF.

I seem to recall that part of the problem with the using the scikit-image routine was the definition of the point spread function which is off-centered by a half pixel in the XRT point spread images.

I don't think I fully understand this. Could this not be solved by shifting the XRT point spread function by half a pixel?

[jslavin]

As I recall, I did find improved results by shifting the point spread function - actually just cutting off a column and row such that the convolution was done using an image with an odd number of pixels and the central one being the peak of the point spread function. However this is all a bit hazy in my memory so I'd have to look back at what I did. I'll try to do that as soon as I get the chance.