[OpenCV complete routine] 90 Frequency domain notch filter

Posted by dirkdetken on Sat, 05 Feb 2022 01:05:49 +0100

[OpenCV complete routine] 90 Frequency domain notch filter

Welcome to pay attention "100 complete OpenCV routines" Series, continuously updating
Welcome to pay attention "Python Xiaobai's OpenCV learning course" Series, continuously updating

5.2 Notch Filter

Notch filter blocks or passes through the frequency in the rectangular neighborhood of the predetermined frequency, which is an important selective filter.

Notch filter is a filter that can rapidly attenuate the input signal at a certain frequency point to achieve the filtering effect of hindering the passage of this frequency signal. Notch filter is a kind of band stop filter. Its stop band is very narrow, and the starting order must be more than second order (including second order).

The zero phase shift filter must give symmetry to the origin (the center of the frequency rectangle), so the center is ( u 0 , v 0 ) (u_0,v_0) The notch filter transfer function of (u0, v0) is ( − u 0 , − v 0 ) (-u_0, -v_0) (− u0, − v0) position must have a corresponding notch.

The transfer function of the notch band stop filter can be constructed by the product of the high pass filter whose center has been shifted to the center of the notch filter:
H N R ( u , v ) = ∏ k = 1 Q H k ( u , v ) H − k ( u , v ) H_{NR}(u,v) = \prod_{k=1}^Q H_k(u,v) H_{-k}(u,v) HNR(u,v)=k=1∏QHk(u,v)H−k(u,v)

The distance calculation formula of the filter is:
D k ( u , v ) = ( u − M / 2 − u k ) 2 + ( v − N / 2 − v k ) 2 D − k ( u , v ) = ( u − M / 2 + u k ) 2 + ( v − N / 2 + v k ) 2 D_k(u,v) = \sqrt{(u-M/2-u_k)^2 + (v-N/2-v_k)^2} \\ D_{-k}(u,v) = \sqrt{(u-M/2+u_k)^2 + (v-N/2+v_k)^2} Dk(u,v)=(u−M/2−uk)2+(v−N/2−vk)2 D−k(u,v)=(u−M/2+uk)2+(v−N/2+vk)2
Therefore, the third-order Butterworth notch band stop filter is:
H N R ( u , v ) = ∏ k = 1 3 [ 1 1 + [ D 0 k / D k ( u , v ) ] n ] [ 1 1 + [ D − k / D k ( u , v ) ] n ] H_{NR}(u,v) = \prod_{k=1}^3 [\frac{1}{1+[D_{0k}/D_k(u,v)]^n}] [\frac{1}{1+[D_{-k}/D_k(u,v)]^n}] HNR(u,v)=k=1∏3[1+[D0k/Dk(u,v)]n1][1+[D−k/Dk(u,v)]n1]
The transfer function of notch bandpass filter can be obtained by Notch bandstop filter:
H N P ( u , v ) = 1 − H N R ( u , v ) H_{NP}(u,v) = 1 - H_{NR}(u,v) HNP(u,v)=1−HNR(u,v)

Routine 8.29 uses notch filtering to delete moire patterns in digital printed images

This example uses notch filtering to reduce the moire ripple in the digital printed image.

Moire pattern is the expression of beat difference principle. After the superposition of two equal amplitude sine waves with similar frequencies, the signal amplitude will change according to the difference between the two frequencies. If the spatial frequency of the pixel of the photosensitive element is close to the spatial frequency of the stripe in the image, the moire will be generated.

Notch filtering is to selectively modify the local region of DFT. The typical processing method is interactive operation. Directly use the mouse to select the rectangular region in the Fourier spectrum and find the maximum point as (uk,vk). In order to simplify the program, this routine deletes the mouse interaction part and only retains the filtering process after (uk,vk).

# OpenCVdemo08.py
# Demo08 of OpenCV
# 8. Image frequency domain filtering
# Copyright 2021 Youcans, XUPT
# Crated: 2021-12-30  
   
    # Routine 8.29 uses notch filtering to delete moire patterns in digital printed images

    def gaussLowPassFilter(shape, radius=10):  # Gaussian low pass filter
        u, v = np.mgrid[-1:1:2.0 / shape[0], -1:1:2.0 / shape[1]]
        D = np.sqrt(u ** 2 + v ** 2)
        D0 = radius / shape[0]
        kernel = np.exp(- (D ** 2) / (2 * D0 ** 2))
        return kernel

    def butterworthNRFilter(shape, radius=9, uk=60, vk=80, n=2):  # Butterworth notch band stop filter
        M, N = shape[1], shape[0]
        u, v = np.meshgrid(np.arange(M), np.arange(N))
        Dm = np.sqrt((u - M // 2 - uk) ** 2 + (v - N // 2 - vk) ** 2)
        Dp = np.sqrt((u - M // 2 + uk) ** 2 + (v - N // 2 + vk) ** 2)
        D0 = radius
        n2 = n * 2
        kernel = (1 / (1 + (D0 / (Dm + 1e-6)) ** n2)) * (1 / (1 + (D0 / (Dp + 1e-6)) ** n2))
        return kernel

    def imgFrequencyFilter(img, lpTyper="GaussLP", radius=10):
        normalize = lambda x: (x - x.min()) / (x.max() - x.min() + 1e-8)

        # (1) Edge fill
        imgPad = np.pad(img, ((0, img.shape[0]), (0, img.shape[1])), mode="reflect")
        # (2) Centralization: F (x, y) * - 1 ^ (x + y)
        mask = np.ones(imgPad.shape)
        mask[1::2, ::2] = -1
        mask[::2, 1::2] = -1
        imgPadCen = imgPad * mask
        # (3) Fourier transform
        fft = np.fft.fft2(imgPadCen)

        # (4) Construct the frequency domain filter transfer function:
        if lpTyper == "GaussLP":
            print(lpTyper)
            freFilter = gaussLowPassFilter(imgPad.shape, radius=60)
        elif lpTyper == "GaussHP":
            freFilter = gaussLowPassFilter(imgPad.shape, radius=60)
        elif lpTyper == "ButterworthNR":
            print(lpTyper)
            freFilter = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=80, n=2)  # Butterworth notch band stop filter
        elif lpTyper == "MButterworthNR":
            print(lpTyper)
            BNRF1 = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=80, n=2)  # Butterworth notch band stop filter
            BNRF2 = butterworthNRFilter(imgPad.shape, radius=9, uk=-60, vk=80, n=2)
            BNRF3 = butterworthNRFilter(imgPad.shape, radius=9, uk=60, vk=160, n=2)
            BNRF4 = butterworthNRFilter(imgPad.shape, radius=9, uk=-60, vk=160, n=2)
            freFilter = BNRF1 * BNRF2 * BNRF3 * BNRF4
        else:
            print("Error of unknown filter")
            return -1

        # (5) Modify Fourier transform in frequency domain: Fourier transform point multiplication filter transfer function
        freTrans = fft * freFilter
        # (6) Inverse Fourier transform
        ifft = np.fft.ifft2(freTrans)
        # (7) Decentralized inverse transform image
        M, N = img.shape[:2]
        mask2 = np.ones(imgPad.shape)
        mask2[1::2, ::2] = -1
        mask2[::2, 1::2] = -1
        ifftCenPad = ifft.real * mask2
        # (8) Intercept the upper left corner, the size is equal to the input image
        imgFilter = ifftCenPad[:M, :N]
        imgFilter = np.clip(imgFilter, 0, imgFilter.max())
        imgFilter = np.uint8(normalize(imgFilter) * 255)
        return imgFilter


    # Using notch filtering to delete moire patterns in digital printed images
    # (1) Read original image
    img = cv2.imread("../images/Fig0464a.tif", flags=0)  # flags=0 read as grayscale image
    fig = plt.figure(figsize=(10, 5))
    plt.subplot(141), plt.title("Original"), plt.axis('off'), plt.imshow(img, cmap='gray')

    # (2) Image Gaussian low pass filtering
    imgGLPF = imgFrequencyFilter(img, lpTyper="GaussLP", radius=30)  # Image Gaussian low pass filtering
    plt.subplot(142), plt.title("GaussLP filter"), plt.axis('off'), plt.imshow(imgGLPF, cmap='gray')

    # (3) Image Butterworth notch band stop filtering
    imgBNRF = imgFrequencyFilter(img, lpTyper="ButterworthNR", radius=9)
    plt.subplot(143), plt.title("ButterworthNR filter"), plt.axis('off'), plt.imshow(imgBNRF, cmap='gray')

    # (4) Superimposed Butterworth notch band stop filtering
    imgSBNRF = imgFrequencyFilter(img, lpTyper="MButterworthNR", radius=9)
    plt.subplot(144), plt.title("Superimposed BNRF"), plt.axis('off'), plt.imshow(imgSBNRF, cmap='gray')

    plt.tight_layout()
    plt.show()

(end of this section)

youcans@xupt Original works, reprint must be marked with the original link

Crated: 2022-2-1

Welcome to pay attention "100 complete OpenCV routines" Series, continuously updating
Welcome to pay attention "Python Xiaobai's OpenCV learning course" Series, continuously updating

Topics: Python OpenCV Algorithm Computer Vision image processing

Programmer Think

[OpenCV complete routine] 90 Frequency domain notch filter