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

Question

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[0]
    fft_cols = blurred_image_fft_shift.shape[1]
    first_func = np.zeros((fft_rows,fft_cols), dtype = complex)
            
    denom = np.multiply(np.add(H_abs_sq,constant),transfer_func)
    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 = cv.imread("blurred.tif",0)
    img_rows = img.shape[0]
    img_cols = img.shape[1]
    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()

Output:

Answer 1

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.