我想用点扩散函数(PSF)去卷积2D图像.我已经看到有一个scipy.signal.deconvolve函数适用于一维数组,而scipy.signal.fftconvolve适用于多维数组. scipy中是否有特定的函数来解卷积2D数组？
我已经定义了一个fftdeconvolve函数替换fftconvolve中的乘积除以：
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)
但是,以下代码在解卷和解卷后不会恢复原始信号：
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()
解决方法
使用来自SciPy的fftpack包的fftn,ifftn,fftshift和ifftshift的这些函数似乎有效：
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()
