Removing periodic noise from image using Fourier transform

Question

I'm trying to get rid of some periodic flicker noise through post-processing of the recorded images. The reason for these artifacts is that the electronic rolling shutter of the camera reads each line sequentially, and this mechanism results in brightness variations in the image, when the picture is taken under fluorescent lighting (due to 50 or 60 Hz AC power frequency).

I tried calculating the Fourier transform of the image and suppressing the components that might cause flickering. I inspected the magnitude spectrum (as you can see below for an example image) and set the vertical components around the DC component to a very small value (I preserved the DC component). In the end, I managed to eliminate the flicker because I suppressed the flicker frequency, but also the image was sort of deformed since I, probably, also suppressed some other frequencies that contributed to the signal's energy.

I would like to know what would be a better way to detect the flicker frequency and eliminate/suppress only that frequency? In my current implementation, I could be doing something fundamentally wrong, so I would appreciate any guidance and hints. Below, I show the resulting image that I got for the kind of filtering that I explained above. You can also find my Python code below.

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

img = cv2.imread('images/flicker2.jpg',0)
f = np.fft.fft2(img)
fshift = np.fft.fftshift(f)
# calculate amplitude spectrum
mag_spec = 20*np.log(np.abs(fshift))

r = f.shape[0]/2        # number of rows/2
c = f.shape[1]/2        # number of columns/2   
p = 3                         
n = 1                   # to suppress all except for the DC component       
fshift2 = np.copy(fshift)

# suppress upper part
fshift2[0:r-n , c-p:c+p] = 0.001
# suppress lower part
fshift2[r+n:r+r, c-p:c+p] = 0.001
# calculate new amplitude spectrum
mag_spec2 = 20*np.log(np.abs(fshift2))
inv_fshift = np.fft.ifftshift(fshift2)
# reconstruct image
img_recon = np.real(np.fft.ifft2(inv_fshift))

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(mag_spec, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(mag_spec2, cmap = 'gray')
plt.title('Magnitude Spectrum after suppression'), plt.xticks([]), plt.yticks([])
plt.show()

plt.figure()
plt.subplot(121),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(122),plt.imshow(img_recon, cmap = 'gray')
plt.title('Output Image'), plt.xticks([]), plt.yticks([])
plt.show()

Answer 1

Just eyeballing the first image tells me that you have a very simple band pattern. They're in fact mostly horizontal (in the spatial domain). You see a clear cross in the Fourier domain. That's because most of your spatial features in that hallway are oriented horizontally and vertically. Your horizontal band similarly show up as a pair of dots above and below the origin. The bands are slightly rotated to the right, and so will the pair of dots.

You'll need to erase a very small area around these dots. But how to locate them? The easiest solution may be to take your recorded image, open it in Paint or a similar program, and paint pure black and white bars over the actual bands. Now take the FFT of this bar pattern, and you get the exact location of the two points you're looking for. It in general won't be exactly a single FFT bin, so you could fit a gaussian curve to the neighborhood of those two points, and use the opposite 1-gauss(x-xc,y-yc) as the frequency domain filter.