Potential Enhancements and Clarifications for Wiener Deconvolution (restoration.wiener)

[nombac]

Description:

I would like to suggest a couple of enhancements or clarifications for the restoration.wiener function within scikit-image, especially in regard to the treatment of the transfer function (PSF) in Fourier space and the handling of complex Point Spread Functions (PSFs). However, I'd like to preface these suggestions with the acknowledgment that there may be aspects I've misunderstood, and I am open to clarification.

Firstly, the documentation currently specifies that the transfer function's shape should be (N1, N2, ..., ND) when the is_real parameter is set to True, and (N1, N2, ..., ND // 2 + 1) otherwise. This guidance seems somewhat counterintuitive because, based on the Hermitian symmetry principle for real-valued signals, only half of the frequency components (up to ND // 2 + 1 for the last dimension) should be necessary due to redundancy. I wonder if this might be a misinterpretation on my part or if there's a potential oversight in the documentation. I propose revising the documentation or the function's behavior to reflect what seems to be a more efficient approach: (N1, N2, ..., ND // 2 + 1) when is_real=True.

Secondly, it has come to my attention that the implementation may discard the imaginary part of a complex PSF, focusing only on the real component (psf = psf.real.astype(float_type, copy=False)). This approach could potentially omit critical phase information inherent in the imaginary part, which might be essential for accurate deconvolution. I am curious if this is an intentional aspect of the design for specific use cases, or if it might be an area for improvement. Enhancing the function to fully utilize the complex nature of PSFs could ensure that all vital information is preserved for more precise deconvolution outcomes.

These suggestions aim to enhance the utility and accuracy of scikit-image for image restoration tasks, which are crucial in various scientific research domains. It is essential to address these points to ensure that users can achieve reliable and precise results. I am open to the possibility that my understanding may not be complete, and I welcome any clarifications or discussions on these matters.

[imagesc-bot]

This issue has been mentioned on Image.sc Forum. There might be relevant details there:

https://forum.image.sc/t/wiener-deconvolution-implementation-handling-complex-psf-in-scikit-image/94170/4

[lagru]

[...] the implementation may discard the imaginary part of a complex PSF, focusing only on the real component (psf = psf.real.astype(float_type, copy=False)). This approach could potentially omit critical phase information inherent in the imaginary part, which might be essential for accurate deconvolution. [...]

As stated in my comment, I’m not really familiar with that function but after some Git sleuthing it looks like it might have been an oversight that was introduced in single precision support in skimage.restoration module #5219. I haven’t checked the PR and context in detail though.

It seems like instead of always casting the real part to float_type it should check if reg and psf are complex and then use the corresponding complex type instead. E.g. float32 → complex64…

Code snippet for investigation

Not sure if this actually makes sense in practice actually!

import skimage as ski
import scipy as sp
import numpy as np

img = ski.color.rgb2gray(ski.data.astronaut())
psf = np.ones((5, 5)) / 25
img = sp.signal.convolve2d(img, psf, 'same')
rng = np.random.default_rng()
img += 0.1 * img.std() * rng.standard_normal(img.shape)

psf = psf * (1+1j)

deconvolved_img = ski.restoration.wiener(img, psf, 0.1)

[lagru]

Optimistically pinging @jni, @grlee77 and @FirefoxMetzger in the hope they have a deeper understanding of this and some quick insights into the potential bug. 🙏

[imagesc-bot]

This issue has been mentioned on Image.sc Forum. There might be relevant details there:

https://forum.image.sc/t/clarification-needed-shape-of-transfer-function-in-restoration-wiener-with-is-real-parameter/94173/2

[nombac]

It appears that the issue can be addressed by making the following two modifications:

1. Comment out the following two lines:

   # psf = psf.real.astype(float_type, copy=False)
   # reg = reg.real.astype(float_type, copy=False)

2. Modify this line:

   if is_real:
       deconv = uft.uirfftn(wiener_filter * uft.urfftn(image), shape=image.shape)  
   else:
       deconv = uft.uifftn(wiener_filter * uft.ufftn(image))

   to:

   if not is_real:
       deconv = np.real(uft.uifftn(wiener_filter * uft.ufftn(image)))
   else:
       deconv = uft.uirfftn(wiener_filter * uft.urfftn(image), shape=image.shape)

The second modification includes logical clarification.

[nombac]

I've identified another potential bug. When the "psf" (Point Spread Function) is provided with the same shape as the "image", it's directly passed to "trans_func" without being transformed by "ir2tf":

if psf.shape != reg.shape:
    trans_func = uft.ir2tf(psf, image.shape, is_real=is_real)
else:
    trans_func = psf

To ensure that the "psf" is always correctly transformed, regardless of its shape, I suggest modifying the condition to check if the "psf" is a complex object. If it's not, it should be transformed using "ir2tf". Here's the proposed change:

if not np.iscomplexobj(psf):
    trans_func = uft.ir2tf(psf, image.shape, is_real=is_real)
else:
    trans_func = psf

This change ensures that the "psf" is always in the correct format for processing, whether it initially matches the "image" shape or not.

[nombac]

Below is the "git diff" highlighting the minimal changes needed. From my testing, these modifications work well for my specific use cases.

--- a/skimage/restoration/deconvolution.py
+++ b/skimage/restoration/deconvolution.py
@@ -116,10 +116,8 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
         reg = uft.ir2tf(reg, image.shape, is_real=is_real)
     float_type = _supported_float_type(image.dtype)
     image = image.astype(float_type, copy=False)
-    psf = psf.real.astype(float_type, copy=False)
-    reg = reg.real.astype(float_type, copy=False)
 
-    if psf.shape != reg.shape:
+    if not np.iscomplexobj(psf):
         trans_func = uft.ir2tf(psf, image.shape, is_real=is_real)
     else:
         trans_func = psf
@@ -130,7 +128,7 @@ def wiener(image, psf, balance, reg=None, is_real=True, clip=True):
         deconv = uft.uirfftn(wiener_filter * uft.urfftn(image),
                              shape=image.shape)
     else:
-        deconv = uft.uifftn(wiener_filter * uft.ufftn(image))
+        deconv = np.real(uft.uifftn(wiener_filter * uft.ufftn(image)))
 
     if clip:
         deconv[deconv > 1] = 1