Numpy – is there an equivalent of SciPy for 2D arrays signal. Deconvolve stuff?

I want to deconvolute 2D images with point spread function (PSF) I've seen a SciPy signal. The deconvolve function applies to one-dimensional arrays, while SciPy signal. Fftconvolve applies to multidimensional arrays Is there a specific function in SciPy to deconvolute 2D arrays?

I have defined a fftdeconvolve function to replace the product in fftconvolve divided by:

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 following code will not recover the original signal after unwinding:

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),4)
psf = stats.norm.pdf(np.sqrt((X - sx/2)**2 + (Y - sy/2)**2),10)

star_conv = fftconvolve(star,psf,mode="same")
star_deconv = fftdeconvolve(star_conv,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 array is shown in the next row in the figure below How to use FFT deconvolution to recover the original signal?

Solution

These functions of fftn, ifftn, fftshift and ifftshift using the fftpack package from SciPy seem to be valid:

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,0].imshow(np.real(star_conv))
axes[1,1].imshow(np.real(star_deconv))
plt.show()

The image in the lower left corner shows the convolution of two Gaussian functions in the uplink, and the lower right corner shows the inversion of the convolution effect