CODEWITHSUNDEEP
pythonnumpyimage-processingfftgaussian
Elder
Elder

Reputation: 31

Gaussian notch filter in Python

I'm trying to design a Gaussian notch filter in Python to remove periodic noise.I tried implementing the following formula:

Gaussian Notch Filter

And here is the code:

import numpy as np

def gaussian_bandpass_filter(image):
    image_array = np.array(image)
    #Fourier Transform
    fourier_transform = np.fft.fftshift(np.fft.fft2(image_array)) 

    #Size of Image
    m = np.shape(fourier_transform)[0]
    n = np.shape(fourier_transform)[1]

    u = np.arange(m)
    v = np.arange(n)

    # Find the center
    u0 = int(m/2)
    v0 = int(n/2)

    # Bandwidth
    D0 = 10

    gaussian_filter = np.zeros(np.shape(fourier_transform))

    for x in u:
        for y in v:
            D1 = math.sqrt((x-m/2-u0)**2 + (y-n/2-v0)**2)
            D2 = math.sqrt((x-m/2+u0)**2 + (y-n/2+v0)**2)
            gaussian_filter[x][y] = 1 - math.exp(-0.5 * D1*D2/(D0**2))

    #Apply the filter
    fourier_transform = fourier_transform + gaussian_filter
    image_array = np.fft.ifft2(np.fft.ifftshift(fourier_transform))

    return image_array

this function is supposed to apply the Gaussian notch filter to an image and return the filtered image but it doesn't seem to work. I don't know where I went wrong with this (maybe I didn't understand the formula correctly?) so if anyone could help me I would really appreciate it.

Edit:

As an example, here is a noisy image. Using the existing gaussian_filter function in scipy.ndimage library, I get this, which is acceptable. But my function returns this. (I'm using PIL.Image.fromarray function to convert array to image)

Upvotes: 3

Views: 1961

Answers (1)

Sandipan Dey
Sandipan Dey

Reputation: 23129

Note that frequency response of a Gaussian Band-Reject Filter (notch) function is defined as:

enter image description here

where D0 is the cutoff frequency and W is the bandwidth of the filter. Implementing the same with the function gaussian_bandstop_filter() (it seems there is some issue with the frequency response formula you used) and using vectorization instead of loop for faster performance,

from skimage.io import imread
import numpy as np
import matplotlib.pylab as plt

def gaussian_bandstop_filter(image):
    #Fourier Transform
    image_freq = np.fft.fftshift(np.fft.fft2(image)) 
    #Size of Image
    m, n = np.shape(image_freq)
    # Find the center
    u0, v0 = int(n/2), int(m/2)
    u, v = np.meshgrid(range(n), range(m))
    D0, W = 4, 10 # cutoff and bandwidth for the BRF
    D = np.sqrt((u-u0)**2+(v-v0)**2)
    gaussian_filter = np.exp(-((D**2-D0**2)/(D*W))**2)
    #Apply the filter
    fourier_transform = image_freq * gaussian_filter
    rec_image = np.fft.ifft2(np.fft.ifftshift(fourier_transform))
    return rec_image, image_freq, gaussian_filter

def plot_spectrum(F, ax, method=None, cmap='coolwarm'):
    F_disp = (20*np.log10(0.01 + np.abs(F))).astype(int) if method == 'log' else np.minimum(np.abs(F).astype(int)+100, 255) 
    ax.imshow(F_disp, cmap=cmap, aspect='auto')

im = imread('images/USc6e.png', as_gray=True) 
im = im / im.max() # normalize 
im_rec, F, G = gaussian_bandstop_filter(im)
im_rec = im_rec.real
# plot images
fig, axes = plt.subplots(nrows=2, ncols=2, figsize=(12,10))
axes = axes.ravel()
plt.gray()
axes[0].imshow(im), axes[0].axis('off')
axes[0].set_title('image with periodic noise', size=15)
plot_spectrum(F.real, axes[1]), axes[1].set_title('noisy image spectrum', size=15)
axes[2].imshow(im_rec/im_rec.max()), axes[2].axis('off')
axes[2].set_title('image with band-reject filter', size=15)
axes[3].imshow(G), axes[3].set_title('band reject filter spectrum', size=15)
plt.tight_layout()
plt.show()

The next figure shows the result obtained, along with the frequency spectrum of the noisy image and the BRF:

enter image description here

You may need to fine-tune the parameters D0 and W (see the frequency spectrum of the noisy image).

Upvotes: 0

Related Questions