Esiste un equivalente di scipy.signal.deconvolve per gli array 2D?

Question

Vorrei deconvolgere un'immagine 2D con una funzione di diffusione puntiforme (PSF). Ho visto che esiste una funzione scipy.signal.deconvolve che funziona per gli array monodimensionali e scipy.signal.fftconvolve per la convalida di array multidimensionali. Esiste una funzione specifica in Scipy per deconvolgere gli array 2D?

ho definito una funzione fftdeconvolve sostituzione del prodotto in fftconvolve da un 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) 

Tuttavia, il codice appresso non recuperare il segnale originale dopo convoluzione e 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() 

L' gli array 2D risultanti sono mostrati nella riga in basso nella figura sottostante. Come posso recuperare il segnale originale usando la deconvoluzione FFT?

Answer 1

Queste funzioni utilizzando fftn, ifftn, fftshift e ifftshift dalla confezione fftpack del SciPy sembrano funzionare:

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

L'immagine in basso a sinistra mostra la convoluzione delle due funzioni gaussiane in la riga superiore e il retro degli effetti della convoluzione sono mostrati in basso a destra.

Answer 2

noti che deconvolving per divisione nel dominio di Fourier non è davvero utile per qualsiasi cosa, ma a scopo dimostrativo; qualsiasi tipo di rumore, anche numerico, può rendere il tuo risultato completamente inutilizzabile. Si può regolarizzare il rumore in vari modi; ma nella mia esperienza, un RL iteration è più facile da implementare e in molti modi più giustificabile dal punto di vista fisico.