Wiener deconvolution and psf

[acycliq]

Component
PyDIP

Describe the problem that this feature would solve
Hi Cris

Well, this is not really a feature request but it doesnt fit under bugs either. I am trying to understand how to use the WienerDeconvolution function with a toy example but apparently I am missing something quite basic.

More specifically I saw your answer to this post on SO and I want to replicate it with a single call to dip.WienerDeconvolution passing-in the image and the psf only as follows:

import diplib as dip

image = dip.ImageRead('trui.ics');
kernel = dip.CreateGauss([3,3])

smooth = dip.ConvolveFT(image, kernel)
smooth = dip.GaussianNoise(smooth, 5.0)  # variance = 5.0
out = dip.WienerDeconvolution(smooth, kernel)

My out however doesnt look anything close to yours (the out from stackOverflow)

Is it the kernel (psf) I have misunderstood?

[crisluengo]

@acycliq There's a regularization parameter that you should set. The default value apparently is too small for this case. If it's too small, the output will be too noisy, as you see in your output. If it's too large, the output will be too smooth (i.e. not sharpened enough). Usually you find the right value with a bit of trial-and-error, change the value in factors of ten first.

[acycliq]

Yes, thanks! Apparently setting it to 0.1 yields something really close to your result , dip.WienerDeconvolution(smooth, kernel, 0.1)

Is it only by trial and error that you can set this please? Any pointers to help getting a descent estimate?

Thanks so much!

[crisluengo]

You can estimate the parameter if you know the SNR of the image.

You could uses dip.EstimateNoiseVariance() to attempt to estimate the noise variance (which is the noise power). The signal power can be estimated from the input image assuming it has little noise:

K = dip.EstimateNoiseVariance(smooth) / dip.MeanSquareModulus(smooth)[0][0]

This doesn't seem to provide the right values though, I'd have to look into it in more detail to see why.

[acycliq]

Yes thanks Cris, either way you (ie the user) should somehow come up with some reasonable values for the signal power (or maybe the gaussian blur) or the regularization parameter for the data at hand, Makes sense..

Last question please, these functions like dip.WienerDeconvolution will work for a 3d image, right?

I would find these few lines from your SO answer quite useful in the deconvolution notebook you have written under /examples/python. There was also another post on SO about segmentation where there was a rectangular on a disk and you had come up with a very nice solution

[crisluengo]

Yes, all the deconvolution functions work for any number of dimensions. Note that there are other, better, deconvolution algorithms implemented. The Wiener filter is rather simplistic...

[acycliq]

Thanks again Cris for all your help!