2013-07-04 8 views
7

2次元画像を点広がり関数(PSF)でデコンボリューションしたいと思います。私は、1次元配列のために働くscipy.signal.deconvolve関数があり、多次元配列を畳み込むのにscipy.signal.fftconvolveがあることを見ました。 scipyに2D配列のデコンボリューションを行うための特定の関数はありますか?2次元アレイに対応するscipy.signal.deconvolveはありますか?

Iはdivistionによってfftconvolveに製品を交換fftdeconvolve関数を定義した:

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), 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デコンボリューションを使用して元の信号を復元するにはどうすればよいですか? scipyのダウンロードのFFTPACKパッケージからfftn、関数ifftn場合、fftshiftとifftshiftを使用してこれらの機能が動作するように見える

enter image description here

答えて

7

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

画像左下にある2つのガウス関数の畳み込みを示し上の行、畳み込みの影響の逆が右下に表示されます。

enter image description here

3

フーリエ領域に分割して逆畳み込みするデモの目的が、何のために本当に便利ではないことに注意してください。どんな種類のノイズでも数値であっても、結果は完全に使用できなくなる可能性があります。さまざまな方法でノイズを正規化することができます。私の経験では、RL iterationは実装が簡単で、多くの点で物理的に正当なものです。

+0

水域を持つfドメインの分割について、 'spec_out = spec_signal/spec_to_deconv'ではなく、' spec_out = spec_signal /(spec_to_deconv + waterlevel) 'とはどういう意味ですか? 'spec_to_deconv'の小さな値の摂動を避け、' spec_out'に大きな変化をもたらします。 – Lee

+1

これは基本的にウィーナーフィルタの背後にあるアイデアです。それは最高の激しい物理的正当性を有しており、閾値を設定することは、常に、雑音抑圧と導入される歪みとの間の妥協を伴うことになる。 –

+0

RL繰り返しは1d時系列にも適用できますか? – Lee

関連する問題