numpy - vectors - ¿Hay un equivalente de scipy.signal.deconvolucionar para matrices 2D?

numpy zeros

Question

Me gustaría desconvocar una imagen 2D con una función de dispersión de puntos (PSF). He visto que hay una función scipy.signal.deconvolve que funciona para matrices unidimensionales, y scipy.signal.fftconvolve para scipy.signal.fftconvolve matrices multidimensionales. ¿Hay alguna función específica en scipy para desconvertir las matrices 2D?

He definido una función fftdeconvolve que reemplaza el producto en fftconvolve por una división:

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)

Sin embargo, el código siguiente no recupera la señal original después de la convolución y la desconvolución:

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

Las matrices 2D resultantes se muestran en la fila inferior en la figura siguiente. ¿Cómo podría recuperar la señal original usando desconvolución de FFT?

Answer 1

Estas funciones que usan fftn, ifftn, fftshift e ifftshift del paquete fftpack de SciPy parecen funcionar:

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

La imagen en la parte inferior izquierda muestra la convolución de las dos funciones gaussianas en la fila superior, y el reverso de los efectos de la convolución se muestra en la parte inferior derecha.

Answer 2

Tenga en cuenta que desconvocar por división en el dominio de Fourier no es realmente útil para nada más que fines de demostración; cualquier tipo de ruido, incluso numérico, puede hacer que su resultado sea completamente inutilizable. Uno puede regularizar el ruido de varias maneras; pero en mi experiencia, una iteración de RL es más fácil de implementar y, en muchos sentidos, más justificable físicamente.