PIXINSIGHT : Deconvolve with a "moving" radial PSF function ?

Hi John,
Your response is very interesting.
To be honest I haven’t practiced this sort of maths for a long time (I studied engineering… long time ago), so I can’t really challenge the basis hypothesis of convolution (for it is just an integral between two functions, with a change of variable in one of them - to make that moving along the integration range :-) ).
At first I think I could do sort of what you said : calculate several deconvolutions of the same image with various kernels, and combine these results together by a proper mask function. It is probably mathematically incorrect as Deconv(image,kernel 1 + kernel 2) is probably different from Deconv(image, kernel 1) + Deconv(image, kernel 2) - I havent checked :-). On top of that it may works for a one dimensional function of the kernel : you would interpolate between two « useful » extremums, but for a 2 dimensional functions, you need 4 extremums, which complicates a bit more the combination.
I will make some tests, just to see…
By curiosity I have had a look at the Lucy Richardson algorithm : at each iteration you use the kernel to compute a new value for each pixel, using the surrounding pixels of the result of the last iteration weighted by the kernel. I guess it would be possible to call a function that gives you a specific kernel depending on the pixel you are updating ? This function would be calibrated prior to the application of the deconvolution, just like when you calculate your own PSF based on your image before using the deconvolution function of your software. It would the optical center of the image, and a few extremums PSF in a few direction, and then the kernels would be interpolated in between…
Not sure of course it wouldn’t mess the convergence of the algorithm  toward the correct / maximum likelihood solution, but I think it is worth a test : just for me at least to understand why it is wrong if it is :-)
I was wondering what would be my next little Python project, now I know !

Thanks @barrabclaw , very interesting paper, will have a deeper look at it later. I don’t know if they talk about modeling a PSF function that is practical for computing, as they say, or about the algorithm that would use such a function.
By the way, I have seen they mention what @John Hayes suggested : splitting the image in areas where the PSF is considered constant, apply deconvolution with a proper PSF per region, and composite the results by interpolating between regions to avoid artifact at boundaries. This is what I clumsily tried to explain when talking about combining them with a proper mask function, and it seems this method has been explored and is sometime used. I think it is worth a try on my images, but the polar / radial configuration of my PSF does not make it practical, I have to think about it.
Thanks again, definitely worth diving in it !

@barrabclaw I plan to make some tests in both methods but in any case first : model the variation of PSF across the field of view.
I will use PixInsight and save those PSF as images and see how I can distort the center one using regular openCV rotation matrices and skew to obtain the other ones quickly (and therefore have a function that calculate the PSF on the fly given the center PSF, the optical center coordinates and the position at which we evaluate the function). And find an approximation law to that purpose : I guess that for coma something to the square of the distance from the optical center will be required.

I was also wondering if the PSF was the correct approach for coma : when you deconvolve an image that is produced by a corrected / flat field optical device, you try to compensate for dispersion of light in the lenses and diffraction mainly I guess ? Maybe a tiny bit of defocus ? But the physical image at least matches the sensor plane, while in coma, the « metric » is changed, as the image is formed on a curved plane, but recorded / projected on a flat plane. I am wondering if prior to deconvolution, one should try, first, to compensate for that « metric » change ? Like when you try to make planar projection of a spherical map ? I think that can be done very quickly by geometrically and non linearly stretching the image toward the optical center, and resampling the result. And then stick to classical PSF ?
It may be completely stupid, as I am by no way an optical expert nor a mathematician, but would love to read your comments.

Just to join the conversation, your suggestion of using a different PDF for different parts of the image is a very good one, and is completely doable.      Just takes the time and effort to determine how to determine and interpolate the PDF over various areas of the image and actually code it.   Of course we would then have to create PDFs over different parts of the picture, which would be more time-consuming.

For my fellow geeks, deconvolution is simply trying to determine what the signal really is, given the response that we record on our image, assuming that image is a linear (sum of response) to all of the individual signals over the field.   In other words, if we see a response (brightness in a pixel) the only assumption is  that all of the various potential contributors (ie, stars, nebula) that are in the signal contributed to the overall response in a linear way.  For example, the brightness of a given pixel may have 4% contribution from the sky,  40% contribution from the star to the east, and 56% contribution from a star to the west.   Now I just made these numbers up, but the key to deconvolution is determining what these numbers are.   Again the only assumption is that the contributions of every signal (every star) adds up (linearly) to the total response that is seen in the image.

To perform deconvolution, we  need to know what the response is be to a single, known signal.  The details on how the transformation (convolution) occurs is not important.   Once we know that, we can "deconvolve" what the signal needs to be to give the response on our image.   In otherwords what combination of "stars" or signals do we need that when turned into a response by our atmosphere, telescope and camera and added all together, yields the result we see on the image.

In engineering this is often done via a calibration.   Put a known, simple single signal through whatever thing it is you are investigating and see what the response is like.  For most engineering examples, this "calibration" involves putting a simple signal (point source of light, pressure pulse, Dirac-delta, etc.) through a process and seeing what the response is like.

In astronomy, we are fortunate to already have such calibrations in our image.   The signal consists of a star, which is a single point source of light (close enough to have 0 diameter, yet emitting photons) which should appear as a single pixel, of a given brightness in our image if our optical system, camera, atmosphere, whatever else, were perfect.    The response we actually get is a fuzzy circle that is brightest in the centre and tails off with distance away from this.    So we know the signal - a single pixel of given brightness and the response, a fuzzy circle we call the PDF.  That, is the machine (our telescope and camera) show us the point source of light as a fuzzy circle (the PDF).    Note our PDFs also have noise added, but we can get rid of most noise by averaging the calibration - using multiple responses, by averagine the PDFs from multiple stars.  Deconvolution doesn't care why the respose is fuzzy or oval or wide,  only how fuzzy, oval, or wide it is.

Finally, knowing that any brightness (other than noise) in an image must be the sum of all the contributions from all of the stars in the image, we can determine what the original signal must look like to generate the image (response) from the camera.   In a perfect world, the deconvolved image should only show stars that occupy a single pixel, and the image would be perfect recording of the signal.   Instead, for a host of reasons, what we get in practice is an image with smaller stars, sharper nebulae, etc that is more representative of the actual signal, - still not perfect but closer to the signal than the original image.  The "host of reasons" includes residual noise, camera resolution, numerical dispersion, non-linearity of camera response, light diffraction (non-linear), etc. etc.

Using different PDFs for different parts of the picture is no different that using a different PDF for a different picture, a different filter, or a different camera - as long as the signal-response relationship, or calibration distance, or PDF diameter,  is small compared to the overall scale of the picture.   Does a point source of light give a circular image or an oval image and how much spread (standard deviation) results.   Any nebulosity in the image is just the linear sum of all of the "point sources" that make up the nebula.

The beauty of deconvolution is that is doen't need to know anything about the nature of what happens to the signal to get the response and doesn't need to know anything about what created the signal,  just what the transformation from signal to reponse is and that is the PDF.

Hope this entertained you a little.

Dave