2D阵列是否有scipy.signal.deconvolve的等效项? [英] Is there a equivalent of scipy.signal.deconvolve for 2D arrays?

问题描述

我想用点扩展函数(PSF)对2D图像进行反卷积.我已经看到有一个scipy.signal.deconvolve函数可用于一维数组，而scipy.signal.fftconvolve可用于卷积多维数组. scipy中是否有特定功能可以对2D阵列进行反卷积?

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?

我定义了fftdeconvolve函数，以除以fftconvolve中的乘积:

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()

生成的2D数组显示在下图中的下一行.如何使用FFT反卷积恢复原始信号?

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

推荐答案

使用SciPy的fftpack软件包中的fftn，ifftn，fftshift和ifftshift的这些函数似乎起作用:

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.

这篇关于2D阵列是否有scipy.signal.deconvolve的等效项?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！