1

I need to recover the sigma of a gaussian blur. I have both the original image and the smoothed image.

I tried a deconvolution type approach:

import numpy as np


nz = 141
nx = 681
h = 25

x = [i*h / 1000 for i in range(nx)]
z = [i*h / 1000 for i in range(nz)]

vp_true = np.fromfile("vp_true", np.float32).astype(np.float64).reshape(nz, nx, order='F')
vp_smooth = np.fromfile("vp_smooth", np.float32).astype(np.float64).reshape(nz, nx, order='F')

fft_vp_true = np.fft.fft2(vp_true)
fft_vp_smooth = np.fft.fft2(vp_smooth)

fft_g = np.fft.fftshift(fft_vp_smooth) / np.fft.fftshift(fft_vp_true)
g = np.real(np.fft.fftshift(np.fft.ifft2(np.fft.ifftshift(fft_g))))

fig = plt.figure()
ax = plt.gca()

im = ax.pcolorfast(x, z, g, cmap="jet")
ax.invert_yaxis()
ax.set_aspect('equal')

plt.tight_layout()
plt.savefig("test.pdf")
plt.show()

However, the kernel (g) obtained this way is not Gaussian at all enter image description here Am I doing something wrong here ?

Cris Luengo
  • 55,762
  • 10
  • 62
  • 120
Alexandre Hoffmann
  • 75
  • 6
  • 1
    The division you compute is very noise sensitive, your output is dominated by noise. Normally one would apply regularization in this problem. Look for Wiener deconvolution as a simple regularization. – Cris Luengo Sep 08 '22 at 13:15
  • 1
    Does this answer your question? [Finding for convolution kernel if many 0's for FFT?](https://stackoverflow.com/questions/48819161/finding-for-convolution-kernel-if-many-0s-for-fft) – Cris Luengo Sep 08 '22 at 13:16
  • it gives better results thanks – Alexandre Hoffmann Sep 08 '22 at 14:50

0 Answers0