title: python: bandpass filter of an image
body: <p>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). </p>
  
  <p><a href="https://i.stack.imgur.com/AacAP.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/AacAP.png" alt="My data.  The red spots are what I want, the background sine wave in kx+ky is unwanted."></a></p>
  
  <p>In trying to do this, I notice two things: </p>
  
  <p>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?</p>
  
  <pre><code>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)
  </code></pre>
  
  <p>This is an image of img_back:</p>
  
  <p><a href="https://i.stack.imgur.com/DKCEi.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/DKCEi.png" alt="The inverse fourier transform, no filter applied."></a></p>
  
  <p>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.</p>
  
  <p>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.  </p>
  
  <pre><code>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 &gt;= ro*ro 
  m2 = x*x + y*y &lt;= 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)
  </code></pre>
  
  <p>Here the result is the same as before, because np.fft ignores masks.  The fix to that was simple: </p>
  
  <p><code>filtdat2 = filtdat.filled(filtdat.mean())</code></p>
  
  <p>Unfortunately, (but upon reflection also unsurprisingly) the result is shown here: </p>
  
  <p><a href="https://i.stack.imgur.com/9Klkd.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/9Klkd.png" alt="The result of a brickwall bandpass filter."></a></p>
  
  <p>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.</p>
  
  <p>Clearly, the 'brickwall' filter is not the right solution. The phenomenon of making rings from this filter is well explained here: <a href="https://dsp.stackexchange.com/questions/724/low-pass-filter-and-fft-for-beginners-with-python/725#725">What happens when you apply a brick-wall filter to a 1D dataset.</a></p>
  
  <p>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 <a href="http://wiki.scipy.org/Cookbook/ButterworthBandpass" rel="nofollow noreferrer">in this implementation of a butterworth filter</a>.  Possibly the right thing to do involves using fftconvolve() as is done <a href="http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.fftconvolve.html#scipy.signal.fftconvolve" rel="nofollow noreferrer">here to blur an image.</a>  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? </p>