Is there a equivalent of scipy.signal.deconvolve for 2D arrays?

Question

I would like to deconvolve a 2D image with a point spread function (PSF). I've seen there is a scipy.signal.deconvolve function that works for one-dimensional arrays, and scipy.signal.fftconvolve to convolve multi-dimensional arrays. Is there a specific function in scipy to deconvolve 2D arrays?

I have defined a fftdeconvolve function replacing the product in fftconvolve by a divistion:

def fftdeconvolve(in1, in2, mode="full"):
    """Deconvolve two N-dimensional arrays using FFT. See convolve.

    """
    s1 = np.array(in1.shape)
    s2 = np.array(in2.shape)
    complex_result = (np.issubdtype(in1.dtype, np.complex) or
                      np.issubdtype(in2.dtype, np.complex))
    size = s1+s2-1

    # Always use 2**n-sized FFT
    fsize = 2**np.ceil(np.log2(size))
    IN1 = fftpack.fftn(in1,fsize)
    IN1 /= fftpack.fftn(in2,fsize)
    fslice = tuple([slice(0, int(sz)) for sz in size])
    ret = fftpack.ifftn(IN1)[fslice].copy()
    del IN1
    if not complex_result:
        ret = ret.real
    if mode == "full":
        return ret
    elif mode == "same":
        if np.product(s1,axis=0) > np.product(s2,axis=0):
            osize = s1
        else:
            osize = s2
        return _centered(ret,osize)
    elif mode == "valid":
        return _centered(ret,abs(s2-s1)+1)

However, the code below does not recover the original signal after convolving and deconvolving:

sx, sy = 100, 100
X, Y = np.ogrid[0:sx, 0:sy]
star = stats.norm.pdf(np.sqrt((X - sx/2)**2 + (Y - sy/2)**2), 0, 4)
psf = stats.norm.pdf(np.sqrt((X - sx/2)**2 + (Y - sy/2)**2), 0, 10)

star_conv = fftconvolve(star, psf, mode="same")
star_deconv = fftdeconvolve(star_conv, psf, mode="same")

f, axes = plt.subplots(2,2)
axes[0,0].imshow(star)
axes[0,1].imshow(psf)
axes[1,0].imshow(star_conv)
axes[1,1].imshow(star_deconv)
plt.show()

The resulting 2D arrays are shown in the lower row in the figure below. How could I recover the original signal using FFT deconvolution?

enter image description here

score 11 · Answer 1 · edited Jan 14 '14 at 20:41

These functions using fftn, ifftn, fftshift and ifftshift from the SciPy's fftpack package seem to work:

from scipy import fftpack

def convolve(star, psf):
    star_fft = fftpack.fftshift(fftpack.fftn(star))
    psf_fft = fftpack.fftshift(fftpack.fftn(psf))
    return fftpack.fftshift(fftpack.ifftn(fftpack.ifftshift(star_fft*psf_fft)))

def deconvolve(star, psf):
    star_fft = fftpack.fftshift(fftpack.fftn(star))
    psf_fft = fftpack.fftshift(fftpack.fftn(psf))
    return fftpack.fftshift(fftpack.ifftn(fftpack.ifftshift(star_fft/psf_fft)))

star_conv = convolve(star, psf)
star_deconv = deconvolve(star_conv, psf)

f, axes = plt.subplots(2,2)
axes[0,0].imshow(star)
axes[0,1].imshow(psf)
axes[1,0].imshow(np.real(star_conv))
axes[1,1].imshow(np.real(star_deconv))
plt.show()

The image in the left bottom shows the convolution of the two Gaussian functions in the upper row, and the reverse of the effects of convolution is shown in the bottom right.

enter image description here

score 5 · Answer 2 · edited Aug 13 '15 at 12:30

5

Note that deconvolving by division in the fourier domain isn't really useful for anything but demonstration purposes; any kind of noise, even numerical, may render your outcome completely unusable. One may regularize the noise in various ways; but in my experience, an RL iteration is easier to implement, and in many ways more physically justifiable.

edited Aug 13 '15 at 12:30

firegurafiku

3,017
1
28
37

answered Jan 14 '14 at 20:56

Eelco Hoogendoorn

10,459
1
44
42

What about division in f domain with a waterlevel, `spec_out = spec_signal / (spec_to_deconv+waterlevel)` rather than `spec_out = spec_signal / spec_to_deconv`? Avoid small perturbations of small values in `spec_to_deconv` resulting in large variations in `spec_out`. – wsdzbm Jul 18 '16 at 10:31
1

That's essentially the idea behind a Wiener filter; but it has at best strenuous physical justification, and setting the threshold will always involve a compromise between the noise suppression and the distortion which is introduced. – Eelco Hoogendoorn Jul 18 '16 at 10:40
Is the RL iteration applicable to 1d time series too? – wsdzbm Jul 18 '16 at 15:34
Its applicable to any linear inversion problem where both the right-hand side, the solution, and the operator have strictly non-negative coefficients. – Eelco Hoogendoorn Jul 18 '16 at 18:02

Is there a equivalent of scipy.signal.deconvolve for 2D arrays?

2 Answers2

Linked