Deconvolution with Lucy-Richardson method

I am trying to deblur an image using MatLab. I have the point spread function (PSF) that the images was blurred with. Furthermore, I know there is noise that is Gaussian distributed, and the signal to noise ratio (SNR) is very high (>20).

Matlab has a couple of deconvolution functions that use direct filtering (regularized filter and Weiner filter), which do not yield satisfactory results.

MatLab has also the Lucy-Richardson (LR) iterative algorithm that, in my case, does a good job in deblurring the image (judged visually).

My question is: is it theoretically sound to use the LR method when the noise in the image has a Gaussian distribution ?

The LR assumes Poisson noise in the blurred image - but does this mean that it performs best with poisson noise, but may also be adequate to use for other types of noise - so I can expect similar results if I want to deblur similar images in the future?

Or does it mean that the LR may randomly yield non-sense results if used for images that have other types of noise than Poisson ?

Answer 1

[Image Analyst]

Well no and yes. Maybe not theoretically the same, however you know, or you should know, that Poisson noise for anything over an expected value of about 8-10 looks virtually identical to Gaussian noise. Just prove it to yourself by doing a Poisson curve for, say 100 and then check the differences between that and a Gaussian fitted to it. The differences will be very small.

And you pretty much never have just 10 photons per pixel unless you are doing very low exposure photon noise limited experiments like in astronomy or radiography. So for most intents and purposes, Poisson noise manifests itself as Gaussian noise, and so the Gaussian noise you have can be considered Poisson noise. The bottom line is how the image looks and if it improves it enough for you to get the required measurements out of it. And it sounds like it's working for you, so go with it.

[Rozh Al-Mashhdi]

Thank you for your answer.

My images are actually MRI (magnetic resonance) images acquired with a Philips MRI scanner. The average voxel/pixel values are app. 60-180. The standard deviation (noise) is much lower at app. 2-15. But these values have no absolute meaning - they are rather relative values that are actually normalised and rescaled (by the scanner software). The relative measure: the signal to noise ratio (i.e. mean / SD) is usually over 10.

I am confused about this: I know that a poisson with large mean value is app. gaussian. But in my case, the mean value of the signal is much higher than the variance (equivalently the SD) of the signal - as I said mean/SD >10. AFAIK, a poisson distribution has a variance that equals the average value. So does this mean that the noise in my data is not poisson ? Or am I missing something ?

But the images look very convincing after the deconvolution, so it seems that the "deconvlucy" function is doing something right.

PS: I am not using the "READOUT" parameter, so it is left at default 0. I am using the "DAMBPAR" parameter and the value is 2.

[John D'Errico]

I should point out that for a simple, basic Poisson random variable with parameter lambda, it is not the mean and standard deviation that are the same, but the mean and the VARIANCE are the same. For poisson parameter lambda, the mean and variance are both lambda.

https://en.wikipedia.org/wiki/Poisson_distribution

[Image Analyst]

Noise where the mean equals the SD or variance is not necessarily Poisson - it could be anything. In your case, with that high a noise and that low a coefficient of variation, it's almost certainly not Poisson. Ask your instrument maker for confirmation. However that doesn't mean that deconvlucy() can't give you a good result. If it improves your image, then why not go ahead and go with it?

[Rozh Al-Mashhdi]

Thanks for the answers, I think I am close now to understand the solution to my problem.

Just a brief followup, Image Analyst:

Why do you say the noise is high ? The values are not absolute and have therefore no isolated meaning in MRI, as they are arbitrarily rescaled/transformed. The signal to noise ratio is >10 so the noise is actually very low (compared with the signal that is)

At this point my question boils down to: knowing that the noise is in fact gaussian from a theoretical perspective. Not knowing the actual values for signal mean and SD in absolute terms but knowing that the mean of the signal is app. 10 times its SD. Is it still possible the noise may be considered to be app. possionian ?

I just want to make sure that, beside having good looking images (judged subjectively), I also have a theoretically sound reason to use "deconvlucy". This will go in a scientific paper, and it is highly probable the that referees may ask about this point.

[Image Analyst]

Rozh, noise and signal to noise ratio are different. In general a higher signal also means higher noise. Let's say I have a photon noise limited signal with an exposure of 1 ms. It will look very "noisy" even though the noise is very low, say around 5. Now let's say we increase the exposure to 100 ms. We're collecting 100 times as many photon so the signal will be 100 times more. Depending on what kind of noise we have, it could be 100 times greater also. Or it could be some lesser number, but it will most likely be more than the noise of 5 we had with the exposure of 1 ms. So more signal and more noise go hand-in-hand. But the signal to noise ratio may actually increase even though the noise is more. Usually it goes up because usually with more exposure you get a better signal and it "appears" less noisy.

If you had a true Poisson noise process the mean would equal the standard deviation even if the actual values are multiplied by some scaling factor such as in an electronic amplifier and a-to-d converter. However you said that your signal > 10 * SD, so we have like a narrow hump way out there from the origin, not like the humps you see for Poisson distributions where the hump is very broad and goes about way back to the origin. As an extreme to help you visualize it, what if your signal and noise were 1000 + or - 5. Well that's certainly not Poisson - a spread of +/- 5 way out at 1000 certainly does not look like a Poisson distribution. Let's say your signal is s and let's say it's Poisson so the variance is also s. So the SD is sqrt(s) and the SNR is s / sqrt(s) which equals sqrt(s). But your SNR is 10. So is your signal 100? So I say that no, your noise does not seem to be Poisson.

[Rozh Al-Mashhdi]

Thank you for the detailed answer, I appreciate it. That is also what along the lines of my understanding. Strange though that the "deconvlucy" seems to work so well - guess I maybe need to dive in to the algorithm to try figure out why.