Notch Reject Filtering in Python

user_13281578 0
Viewed 2534

I'm trying to implement notch-reject filtering in python for an assignment. I have tried using the notch reject filter formula from Rafael Gonzales book and all I got was a edge detected image. Then I tried ideal notch rejecting and here are the results:

Input image--Output of my program -- Expected output

Here is my code:

import cv2
import numpy as np
import matplotlib.pyplot as plt


def notch_reject_filter(shape, d0=9, u_k=0, v_k=0):
    P, Q = shape
    # Initialize filter with zeros
    H = np.zeros((P, Q))

    # Traverse through filter
    for u in range(0, P):
        for v in range(0, Q):
            # Get euclidean distance from point D(u,v) to the center
            D_uv = np.sqrt((u - P / 2 + u_k) ** 2 + (v - Q / 2 + v_k) ** 2)
            D_muv = np.sqrt((u - P / 2 - u_k) ** 2 + (v - Q / 2 - v_k) ** 2)

            if D_uv <= d0 or D_muv <= d0:
                H[u, v] = 0.0
            else:
                H[u, v] = 1.0

    return H


img = cv2.imread('input.png', 0)

img_shape = img.shape

original = np.fft.fft2(img) 
center = np.fft.fftshift(original)  

NotchRejectCenter = center * notch_reject_filter(img_shape, 32, 50, 50)  
NotchReject = np.fft.ifftshift(NotchRejectCenter)
inverse_NotchReject = np.fft.ifft2(NotchReject)  # Compute the inverse DFT of the result


plot_image = np.concatenate((img, np.abs(inverse_NotchReject)),axis=1)

plt.imshow(plot_image, "gray"), plt.title("Notch Reject Filter")
plt.show()
filter opencv image-processing numpy python
2 Answers
  • all I got was a edge detected image because your implementation was High pass filter which is a black circle in the middle, and that works as Edge detector.
  • Then I tried ideal notch rejecting This is correct if you applied that correctly.

The main concept is to filter the undesired Noise in the frequency domain, the noise can be seen as white spots, and your role is to suppress that white spots by multiplying them by black circles in frequency domain(known as filtering).

to improve this result add more notch filters (H5, H6, ...) to suppress the noise.

results

import cv2
import numpy as np
import matplotlib.pyplot as plt

#------------------------------------------------------
def notch_reject_filter(shape, d0=9, u_k=0, v_k=0):
    P, Q = shape
    # Initialize filter with zeros
    H = np.zeros((P, Q))

    # Traverse through filter
    for u in range(0, P):
        for v in range(0, Q):
            # Get euclidean distance from point D(u,v) to the center
            D_uv = np.sqrt((u - P / 2 + u_k) ** 2 + (v - Q / 2 + v_k) ** 2)
            D_muv = np.sqrt((u - P / 2 - u_k) ** 2 + (v - Q / 2 - v_k) ** 2)

            if D_uv <= d0 or D_muv <= d0:
                H[u, v] = 0.0
            else:
                H[u, v] = 1.0

    return H
#-----------------------------------------------------

img = cv2.imread('input.png', 0)

f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
phase_spectrumR = np.angle(fshift)
magnitude_spectrum = 20*np.log(np.abs(fshift))

img_shape = img.shape

H1 = notch_reject_filter(img_shape, 4, 38, 30)
H2 = notch_reject_filter(img_shape, 4, -42, 27)
H3 = notch_reject_filter(img_shape, 2, 80, 30)
H4 = notch_reject_filter(img_shape, 2, -82, 28)

NotchFilter = H1*H2*H3*H4
NotchRejectCenter = fshift * NotchFilter 
NotchReject = np.fft.ifftshift(NotchRejectCenter)
inverse_NotchReject = np.fft.ifft2(NotchReject)  # Compute the inverse DFT of the result


Result = np.abs(inverse_NotchReject)

plt.subplot(222)
plt.imshow(img, cmap='gray')
plt.title('Original')

plt.subplot(221)
plt.imshow(magnitude_spectrum, cmap='gray')
plt.title('magnitude spectrum')

plt.subplot(223)
plt.imshow(magnitude_spectrum*NotchFilter, "gray") 
plt.title("Notch Reject Filter")

plt.subplot(224)
plt.imshow(Result, "gray") 
plt.title("Result")


plt.show()
user_86182420

Drive by comment, using the for-loop for the notch filter generation is very slow. That operation can be optimized


def notch_reject_filter_vec(shape: tuple[int, int], d0: int, u_k: int, v_k: int):
    (M, N) = shape

    H_0_u = np.repeat(np.arange(M), N).reshape((M, N))
    H_0_v = np.repeat(np.arange(N), M).reshape((N, M)).transpose()

    D_uv = np.sqrt((H_0_u - M / 2 + u_k) ** 2 + (H_0_v - N / 2 + v_k) ** 2)
    D_muv = np.sqrt((H_0_u - M / 2 - u_k) ** 2 + (H_0_v - N / 2 - v_k) ** 2)

    selector_1 = D_uv <= d0
    selector_2 = D_muv <= d0

    selector = np.logical_or(selector_1, selector_2)

    H = np.ones((M, N))
    H[selector] = 0

    return H
user_51244540
Related