title: Inverse Filter of spatially convolved versus frequency convolved image
body: <p>My image processing class has been assigned a project on image restoration. I'm currently working on the Inverse Filter. image -> degrade -> inverse filter -> restore image.  I'm using a simple 5x5 box filter for my degradation.</p>
  
  <p>If I convolve the image in the spatial domain, move to frequency domain, then Inverse Filter the convolved image with the kernel's fft, I get a mess. If I convolve the image in the frequency domain, then Inverse Filter that image, I get a good image. </p>
  
  <p>Frequency domain and spatial domain convolution should be identical. My only thought is I'm doing something wrong with the kernel? I'm using a 5x5 box filter. The spatial convolution divides the final result by np.sum(box). I've tried normalizing the box via:</p>
  
  <pre><code>box = np.ones( 25 ).reshape( 5,5 ) / 25.0
  </code></pre>
  
  <p>but get the same trash Inverse Filtered image result.  </p>
  
  <p>I've also noticed the frequency convolved image ("g_freq.png" from code below) is shifted, probably due to the FFT padding the top and left with the bottom/right of the image. Could this be causing a problem?</p>
  
  <p>Spatial Convolution:
  <img src="https://i.imgur.com/AmcFA.png" alt="spatial convolition"></p>
  
  <p>Frequency Convolution: note the padding along top/left.
  <img src="https://i.imgur.com/44Ldf.png" alt="frequency convolution"></p>
  
  <p>Simplest possible code to create the problem is below. 100% numpy/scipy/matplotlib.</p>
  
  <pre><code>import sys
  import matplotlib
  matplotlib.use( 'Agg' )
  import matplotlib.pyplot as plt
  import numpy as np
  import scipy
  from scipy import ndimage
  
  def save_image( data, filename ) : 
      print "saving",filename
      plt.cla()
      fig = plt.figure()
      ax = fig.add_subplot( 111 )
      ax.imshow( data, interpolation="nearest", cmap=matplotlib.cm.gray )
      fig.savefig( filename )
  
  f = scipy.misc.lena()
  save_image( f, "scipylena.png" )
  
  # create a simple box filter
  kernel = np.ones( 25 ).reshape( 5, 5 ) 
  kernel_padded = np.zeros_like(f,dtype="float")
  # put kernel into upper left
  kernel_padded[:5,:5] = kernel 
  
  # FFT kernel, save as image
  K = np.fft.fftshift( np.fft.fft2( kernel_padded ) )  
  save_image( np.abs(K), "K.png" )
  
  
  # degrade image via spatial convolution
  g = ndimage.convolve( f, kernel )
  if np.sum(kernel) != 0 :
      g /= np.sum(kernel)
  # save spatial image
  save_image( g, "g_spatial.png" )
  
  # take convolved image into frequency domain
  G = np.fft.fftshift( np.fft.fft2( g ) )
  
  # inverse filter the spatially convolved image
  F_HAT = G / K
  
  # back to spatial, save the reconstructed image 
  a = np.nan_to_num( F_HAT )
  f_hat = np.fft.ifft2( np.fft.ifftshift( F_HAT ) )  
  save_image( np.abs( f_hat ), "f_hat_spatial.png" )
  
  # 
  # now the same path but entirely in frequency domain
  #
  
  # create a frequency domain convolved image
  F = np.fft.fftshift( np.fft.fft2( f ) )
  G2 = F * K
  
  # back to spatial, save frequency convolved image  
  g2 = np.fft.ifft2( np.fft.ifftshift( G2 ) )
  save_image( np.abs(g2), "g_freq.png" )
  
  # inverse filter the frequency convolved image
  F_HAT2 = G2 / K
  a = np.nan_to_num( F_HAT2 )
  f_hat2 = np.fft.ifft2( np.fft.ifftshift( a ) ) 
  save_image( np.abs( f_hat2 ), "f_hat_freq.png" )
  </code></pre>
  
  <p>My "f_hat_frequency"
  <img src="https://i.imgur.com/Lg2Sz.png" alt="my f_hat_frequency"></p>
  
  <p>My "f_hat_spatial"  :-(
  <img src="https://i.imgur.com/0eBu2.png" alt="my f_hat_spatial"></p>
  
  <p>Many thanks for any help.</p>
  
  <p>[EDIT] I'm running on Mac OSX 10.6.8 using Numpy 1.6.0 via Enthought's free 32-bit version. (<a href="http://www.enthought.com/products/epd_free.php" rel="nofollow noreferrer">http://www.enthought.com/products/epd_free.php</a>) Python 2.7.2 |EPD_free 7.1-1 (32-bit)</p>
  
  <p>EDIT 31-Oct-2011.
  I think what I'm trying to do has deeper mathematical roots than I understand. 
  <a href="http://www.owlnet.rice.edu/~elec539/Projects99/BACH/proj2/inverse.html" rel="nofollow noreferrer">http://www.owlnet.rice.edu/~elec539/Projects99/BACH/proj2/inverse.html</a> helped a bit. Adding the following to my code before the inverse filter:</p>
  
  <pre><code>H_HAT = np.copy(K)
  np.putmask( H_HAT, H_HAT&gt;0.0001, 0.0001 )
  </code></pre>
  
  <p>gives me an image but with a lot of ringing (probably because of my box filter; need to switch to a Gaussian). Also, the offset of the frequency filtered image is quite likely causing a problem. My prof has looked over my code, can't find a problem. Her suggestion is to continue to use the frequency filtered image rather than the spatially filtered image.</p>
  
  <p>I have a similar question on dsp.stackexchange.com: <a href="https://dsp.stackexchange.com/questions/538/using-the-inverse-filter-to-correct-a-spatially-convolved-image">https://dsp.stackexchange.com/questions/538/using-the-inverse-filter-to-correct-a-spatially-convolved-image</a></p>