# Implementation of Wiener filter to deblur an image using Python and OpenCV

EDIT: I have debugged the runtime warning, and now I am able to get an output image. However, the output image is still blurry. Increasing the constant value heavily distorts the image. Setting constant within a range of .001-.0001 (as recommended in the text) produces the results shown below. My code and output image reflect the changes/the fix for the runtime warning. I believe I have some issue with my implementation of the theory, but am still lost

I am attempting to implement a Wiener Filter to deblur an image based on a specific blurring transfer function (defined in the frequency domain). I know I am doing something wrong either based on the theory or based on my implementation, but I am just not sure what I am missing. My procedure so far is to take the 2D FFT of a blurred image, and to take the fftshift of this (which is $$G(u,v)$$ in the equation). Then to plug this in, as well as the filter transfer function (which is defined in the main method) into the equation provided. From here I take the inverse shift and take the 2D inverse Fourier transform to obtain the original image in the spatial domain.

The equations and definitions I am using are from the 4th edition of Digital Image Processing by Gonzalez and Woods.

The blurring transfer function is defined as (with $$T = 1$$, and $$a=b=0.1$$)

$$H(u,v) = \frac{T}{\pi(ua+vb)} \sin[\pi(ua+vb)]e^{-j\pi(ua+vb)}$$

The equation I am following to obtain the original image is :

$$\hat{F}(u,v) = \left[ \frac{1}{H(u,v)} \frac{|H(u,v)|^2}{|H(u,v)|^2+ K}\right]G(u,v)$$

From the equation for $$F(u,v)$$, $$|H(u,v)|^2$$ is defined as the conjugate of the frequency domain transfer function times the transfer function.

This is my code:

Wiener Function:

def pWienerTF4e(image,H,K):

transfer_func = H
constant = K
H_abs_sq= np.multiply(np.conjugate(transfer_func),transfer_func)
blurred_image = image
blurred_image_fft = np.fft.fft2(blurred_image)
blurred_image_fft_shift = np.fft.fftshift(blurred_image_fft)
fft_rows = blurred_image_fft_shift.shape
fft_cols = blurred_image_fft_shift.shape
first_func = np.zeros((fft_rows,fft_cols), dtype = complex)

for i in range(fft_rows):
for k in range(fft_cols):
if denom[i][k] != 0+0*1j:
first_func[i][k] = H_abs_sq[i][k]/denom[i][k]

unblurred_image_f= np.multiply(first_func,blurred_image_fft_shift)
unblurred_image_f_reshift = np.fft.ifftshift(unblurred_image_f)
unblurred_image = np.fft.ifft2(unblurred_image_f_reshift)
unblurred_image_abs = np.abs(unblurred_image)

return unblurred_image_abs


Main Method:

def main():
img_rows = img.shape
img_cols = img.shape
transfer = np.zeros((img_rows,img_cols), dtype = complex)
for i in range(1,img_rows):
for k in range(1,img_cols):
amplitude = (1/(math.pi*(i*.1+k*.1)))*(math.sin(math.pi*(i*.1+k*.1)))
angle = math.pi*(i*.1+k*.1)
x = amplitude*math.cos(angle)
y = amplitude*math.sin(angle)

transfer[i][k] = x+y*1j

unblurred = pWienerTF4e(img,transfer,.00035)

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(unblurred, cmap = 'gray')
plt.title('Output Image'), plt.xticks([]), plt.yticks([])
plt.show()

plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(unblurred, cmap = 'gray')
plt.title('Output Image'), plt.xticks([]), plt.yticks([])
plt.show()

• Have you tried debugging and looking at denom? You have a problem there as per the runtime warning. From there you can backtrack to the problem in your code...
– Jdip
Nov 23, 2022 at 8:55
• please don't post pictures of code. Instead, copy&paste the code into your question, select it and click the "code format" button, labeled {}. Nov 23, 2022 at 9:57
• And you can also write equations using latex.
– Bob
Nov 23, 2022 at 11:25
• Why is constant an imaginary value? If you need to take the abs of the IFFT, you’re doing something wrong. The imaginary part should be approximately zero, then you take the real part. Nov 23, 2022 at 16:33
• Sorry, it’s not imaginary. It’s hard to read code from images. Still, no reason to make it complex. Also, if you use + and * instead of np.add and np.multiply you’ll get much better readable code. Nov 23, 2022 at 16:35

First of all, note that

$$\frac{1}{H(u,v)} \frac{|H(u,v)|^2}{|H(u,v)|^2+ K} = \frac{H^*(u,v)}{|H(u,v)|^2+ K}$$

If implemented this way, you are not going to have a division by zero, as long as $$K>0$$.

Next, you need to tweak $$K$$ to find the optimal filter. If it’s too small, noise will dominate the result. If it’s too big, not much filtering will be applied. So, if your output is too smooth, decrease the value of $$K$$ further.

Finally, double-check your transfer function. You’re using u, v equal to indices, meaning you only define the transfer function for positive frequencies. It is not symmetric, and therefore will not yield a real-valued inverse transform.

I would use np.fft.fftfreq(img_cols) as the $$u$$ and $$v$$ values. This produces values in the range -0.5 to 0.5, but keeps the origin in the first element. We notice that $$\frac{\sin(\pi x)}{\pi x}$$ is the sinc function, where we divide by 0 we must fill in 1. The transfer function blurs in the $$ua+vb$$ direction only, which matches the blur we observe in the image.

I would use the following code to build the transfer function:

u = np.fft.fftfreq(img_cols)[None, :]
v = np.fft.fftfreq(img_rows)[:, None]
T = 1
a = 0.1
b = 0.1
c = np.pi * (u * a + v * b)  # common part
transfer = T / c * np.sin(c) * np.exp(1j * c)
transfer[c == 0] = 1  # this is where the division by 0 happened


Note that the origin is in the top-left of the image, matching the origin for the FFT. So, do not use np.fft.fftshift and np.fft.ifftshift when computing the Wiener filter.

• Sorry for responding so late. I tried to use u and v as frequencies by taking the the 1D fft along axis 0 of the image for x and along axis 1 for y, and used these values for u and v. However I began to run into an error regarding double scalars when I tried this. Do you have another idea or any thoughts on my method of looking at u and v as frequencies? I really appreciate it! Nov 26, 2022 at 2:39
• @underdawg631 See updated answer. Nov 27, 2022 at 6:57