Inverse filter of spatially minimized and frequency minimized images

Question

My image processing class was assigned an image restoration project. I am currently working on an inverse filter. image -> deprade -> inverse filter -> restore image. I use a simple 5x5 filter for my degradation.

If I drill the image in the spatial domain, go to the frequency domain, then invert Filter thumbnail with fft core, I get a mess. If I drill the image in the frequency domain, then invert this image, I get a good image.

The frequency domain and the convolution of the spatial domain must be the same. My only thought is - am I doing something wrong with the core? I am using a 5x5 filter. Spatial convolution divides the final result into np.sum (box). I tried to normalize the window through:

box = np.ones( 25 ).reshape( 5,5 ) / 25.0 

but get the same inverse filter processing result.

I also noticed that the frequency minimized image ("g_freq.png" from the code below) is shifted, probably due to the fact that the FFT fills the upper and left lower / right parts of the image. Could this be a problem?

Spatial convolution:

Frequency convolution: pay attention to filling along the top / left side.

The simplest possible code to create the problem is given below. 100% numpy / scipy / matplotlib.

 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" ) 

My "f_hat_frequency"

My "f_hat_spatial" :-(

Thanks so much for any help.

[EDIT] I work on Mac OSX 10.6.8 using Numpy 1.6.0 through the free 32-bit version of Enthought. ( http://www.enthought.com/products/epd_free.php ) Python 2.7.2 | EPD_free 7.1-1 (32-bit)

EDIT 31-Oct-2011. I think that what I'm trying to do has deeper mathematical roots than I understand. http://www.owlnet.rice.edu/~elec539/Projects99/BACH/proj2/inverse.html helped a bit. Adding the following code to the code before the reverse filter:

 H_HAT = np.copy(K) np.putmask( H_HAT, H_HAT>0.0001, 0.0001 ) 

gives me an image, but with a lot of calls (perhaps due to my filter box, you need to switch to Gaussian). In addition, the offset of the frequency-filtered image is likely to cause a problem. My professor looked at my code, can not find the problem. Her suggestion is to continue to use a frequency-filtered image rather than a spatially filtered image.

I have a similar question on dsp.stackexchange.com: https://dsp.stackexchange.com/questions/538/using-the-inverse-filter-to-correct-a-spatially-convolved-image

Answer 1

Obviously, the problem is not that F and F_HAT2 not identical. The fact that you need to call nan_to_num is a clear sign that something is going wrong between multiplication and division by K A possible reason is an integer overflow. Try converting F to floating point type after loading.