Inverse Filter of spatially convolved versus frequency convolved image

Question

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.

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.

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:

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

but get the same trash Inverse Filtered image result.

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?

Spatial Convolution:

Frequency Convolution: note the padding along top/left.

Simplest possible code to create the problem is 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" :-(

Many thanks for any help.

[EDIT] I'm running on Mac OSX 10.6.8 using Numpy 1.6.0 via Enthought's free 32-bit version. (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 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 to my code before the inverse 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 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.

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

The problem is clearly that F and F_HAT2 are not identical. The fact that you need to call nan_to_num is a clear indication that something is going wrong between the multiplication and division by K. A possible cause is integer overflow. Try converting f to a floating point type after loading.