CODEWITHSUNDEEP
pythonscipyfftdeconvolution
Tony Power
Tony Power

Reputation: 1178

Deconvolution doesn't seem to work (fft2; scipy)

I'm convolving a signal f, with a kernel h, in spatial domain, then, I'm deconvolving in frequency domain. Since I know the kernel h, in theory I should be able to get f back without any problem, however, it's not working. This is what I have:

    from scipy.fft import fft, ifft, fftfreq, fft2, ifft2, fftshift, ifftshift
    from scipy.ndimage import convolve

    d = 11  # diameter of image
    f = np.zeros((d, d))
    r = 1  # radius of line
    f[d//2-1: d//2+r+1] = 1  # image with a line
    f[d//2, d//2] = 0  # and 0 at the center
    f[0, 0] = 2  # and a 2 at top left corner
    
    print("\nf =")
    print(f)
    
    h = np.zeros_like(f)
    h[d//2-1, d//2-1] = h[d//2, d//2] = h[d//2+1, d//2+1] = 1  # identity filter
    print("\nh =")
    print(h)

    g = convolve(f, h)  # conv spatial domain
    G = fft2(g)
    
    H = fft2(h)  # fft of filter

    f_approx = fftshift(ifft2(G/H))  # deconvolution
    
    print("\nf^approx =")
    print(np.round(f_approx.real, 2))

This is the output:

f =
[[2. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 0. 1. 1. 1. 1. 1.]
 [1. 1. 1. 1. 1. 1. 1. 1. 1. 1. 1.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

h =
[[0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 1. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]
 [0. 0. 0. 0. 0. 0. 0. 0. 0. 0. 0.]]

f^approx =
[[ 0.   -0.67 -0.   -0.   -0.   -0.   -0.   -0.   -0.   -0.    1.33]
 [-0.67  2.    1.33  0.   -0.   -0.   -0.    0.   -0.    0.   -0.  ]
 [-0.    1.33  0.   -0.67 -0.   -0.    0.    0.   -0.    0.   -0.  ]
 [-0.    0.   -0.67 -2.   -0.67  0.   -0.   -0.   -0.   -0.   -0.  ]
 [ 1.    1.    1.    0.33  3.    2.33  1.    1.    1.    1.    1.  ]
 [ 1.    1.    1.    1.    2.33  0.    0.33  1.    1.    1.    1.  ]
 [ 1.    1.    1.    1.    1.    0.33 -1.    0.33  1.    1.    1.  ]
 [-0.    0.    0.    0.   -0.    0.   -0.67  2.    1.33 -0.    0.  ]
 [-0.   -0.    0.   -0.   -0.   -0.   -0.    1.33 -0.   -0.67 -0.  ]
 [ 0.    0.    0.   -0.    0.   -0.    0.    0.   -0.67 -2.   -0.67]
 [ 1.33 -0.   -0.   -0.    0.    0.    0.   -0.   -0.   -0.67  2.  ]]

Why am I not getting the original f back? Am I using the functions incorrectly, or is there some problem with the approach? Thanks

PS: by the way, if the h is the identity, then the deconvolution works fine, the problem is when the filter is something else.

Upvotes: 2

Views: 393

Answers (1)

Tony Power
Tony Power

Reputation: 1178

Thanks @chris-luengo. Also, thanks to Ahmed.

I guess this would be it, then:

    from scipy.fft import fft, ifft, fftfreq, fft2, ifft2, fftshift, ifftshift
    from scipy.ndimage import convolve

    d = 11  # diameter of image
    f = np.zeros((d, d))
    r = 1  # radius of line
    f[d//2-1: d//2+r+1] = 1  # image with a line
    f[d//2, d//2] = 0  # and 0 at the center
    f[0, 0] = 2  # and a 2 at top left corner
    
    print("\nf =")
    print(f)
    
    h = np.zeros_like(f)
    h[d//2-1, d//2-1] = h[d//2, d//2] = h[d//2+1, d//2+1] = 1
    print("\nh =")
    print(h)

    g = convolve(f, h, mode="wrap")  # conv spatial domain
    G = fft2(g)
    
    H = fft2(ifftshift(h))  # fft of filter

    f_approx = ifft2(G/H)  # deconvolution
    
    print("\nf^approx =")
    print(np.round(f_approx.real, 2))

I hope it's useful to someone ;)

Upvotes: 2

Related Questions