python: bandpass filter of an image

Question

I have a data image with an imaging artifact that comes out as a sinusoidal background, which I want to remove. Since it is a single frequency sine wave, it seems natural to Fourier transform and either bandpass filter or "notch filter" (where I think I'd use a gaussian filter at +-omega).

In trying to do this, I notice two things:

1) simply by performing the fft and back, I have reduced the sine wave component, shown below. There seems to be some high-pass filtering of the data just by going there and back?

import numpy as np

f = np.fft.fft2(img)                  #do the fourier transform
fshift1 = np.fft.fftshift(f)          #shift the zero to the center
f_ishift = np.fft.ifftshift(fshift1)  #inverse shift
img_back = np.fft.ifft2(f_ishift)     #inverse fourier transform
img_back = np.abs(img_back)

This is an image of img_back:

Maybe the filtering here is good enough for me, but I'm not that confident in it since I don't have a good understanding of the background suppression.

2) To be more sure of the suppression at the unwanted frequencies, I made a boolean 'bandpass' mask and applied it to the data, but the fourier transform ignores the mask.

a = shape(fshift1)[0]
b = shape(fshift1)[1]

ro = 8
ri = 5
y,x = np.ogrid[-a/2:a/2, -b/2:b/2] 
m1 = x*x + y*y >= ro*ro 
m2 = x*x + y*y <= ri*ri
m3=np.dstack((m1,m2))       
maskcomb =[]
for r in m3:
    maskcomb.append([any(c) for c in r])  #probably not pythonic, sorry
newma = np.invert(maskcomb)
filtdat = ma.array(fshift1,mask=newma) 
imshow(abs(filtdat))
f_ishift = np.fft.ifftshift(filtdat) 
img_back2 = np.fft.ifft2(f_ishift) 
img_back2 = np.abs(img_back2)

Here the result is the same as before, because np.fft ignores masks. The fix to that was simple:

filtdat2 = filtdat.filled(filtdat.mean())

Unfortunately, (but upon reflection also unsurprisingly) the result is shown here:

The left plot is of the amplitude of the FFT, with the bandpass filter applied. It is the dark ring around the central (DC) component. The phase is not shown.

Clearly, the 'brickwall' filter is not the right solution. The phenomenon of making rings from this filter is well explained here: What happens when you apply a brick-wall filter to a 1D dataset.

So now I'm stuck. Perhaps it would be better to use one of the built in scipy methods, but they seem to be for 1d data, as in this implementation of a butterworth filter. Possibly the right thing to do involves using fftconvolve() as is done here to blur an image. My question about fftconvolve is this: Does it require both 'images' (the image and the filter) to be in real space? I think yes, but in the example they use a gaussian, so it's ambiguous (fft(gaussian)=gaussian). If so, then it seems wrong to try to make a real space bandpass filter. Maybe the right strategy uses convolve2d() with the fourier space image and a homemade filter. If so, do you know how to make a good 2d filter?

Answer 1

So, one problem here is that your background sinusoid has a period not terribly different from the signal components you are trying to preserve. i.e., the spacing of the signal peaks is about the same as the period of the background. This is going to make filtering difficult.

My first question is whether this background is truly constant from experiment to experiment, or does it depend on the sample and experimental setup? If it is constant, then background frame subtraction would work better than filtering.

Most of the standard scipy.signal filter functions (bessel, chebychev, etc.) are, as you say, designed for 1-D data. But you can easily extend them to isotropic filtering in 2-D. Each filter in frequency space is a rational function of f. The two representations are [a,b] which are the coefficiets of the numerator and denominator polynomial, or [z,p,k] which is the factored representation of the polynomial i.e.,: H(f) = k(f-z0)*(f-z1)/(f-p0)*(f-p1) You can just take the polynomial from one of the filter design algorithms, evaluate it as a function of sqrt(x^2+y^2) and apply it to your frequency domain data.

Can you post a link to the original image data?