2
$\begingroup$

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: enter image description here

$\endgroup$
9
  • $\begingroup$ 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... $\endgroup$
    – Jdip
    Commented Nov 23, 2022 at 8:55
  • 3
    $\begingroup$ please don't post pictures of code. Instead, copy&paste the code into your question, select it and click the "code format" button, labeled {}. $\endgroup$ Commented Nov 23, 2022 at 9:57
  • $\begingroup$ And you can also write equations using latex. $\endgroup$
    – Bob
    Commented Nov 23, 2022 at 11:25
  • $\begingroup$ 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. $\endgroup$ Commented Nov 23, 2022 at 16:33
  • $\begingroup$ 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. $\endgroup$ Commented Nov 23, 2022 at 16:35

1 Answer 1

1
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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! $\endgroup$ Commented Nov 26, 2022 at 2:39
  • $\begingroup$ @underdawg631 See updated answer. $\endgroup$ Commented Nov 27, 2022 at 6:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.