Inverse and Pseudoinverse Filtering

The Problem

Given a known image function, f(x,y), and a blurring function, h(x,y), we need to recover f(x,y) from the convolution g(x,y)=f(x,y)*h(x,y).

(Note: * is used to represent convolution and x is used to represent multiplication).

A Possible Solution: Inverse Filtering

In the frequency domain we have: G(u,v)=F(u,v) x H(u,v)

F is the FFT2 of the image matrix, f.
H is the FFT2 of the complex distoring function, h.
G is the FFT2 of the "observed" image, g

The novice solution to this problem is to filter by dividing the 2-dimensional fft of the blurred function by the reciprocal of H(u,v):
G(u,v)/H(u,v)=F(u,v) x H(u,v)/H(u,v) = F(u,v).

This is commonly reffered to as the inverse filtering method where 1/H(u,v) is the inverse filter.

Difficulties with Inverse Filtering

The first problem in this formulation is that 1/H(u,v) does not necessairily exist. If H(u,v)=0 or is close to zero, it may not be computationally possible to compute 1/H(u,v).

If there are few values of H(u,v) which are close to zero then the ideal inverse filter can be approximated with a stabilized version of 1/H(u,v) given by :

Fapprox(u,v)=G(u,v) x Hinv(u,v)

where
  
Hinv(u,v)= 1/H(u,v)   if |H(u,v)| > threshold value
        = 0          otherwise

This works well if few elements of H have a magnitude below the threshold but if two many elements are lost, the frequency content of Fapprox will be much lower than F(u,v) and the image will appear distorted.

Here are some examples of how varying the threshold value effects the inverse filtering of an image:

This is the original image:

Image with a 31-Point horizontal blur.

Pseudoinverse filtered image
with 512 out of 65536 values of Hinv=0. Note the "ripples".

Pseudoinverse filtered image with 28682 out of 65536 values of Hinv=0. The image now appears to have horizontal "ghosts" but still provides more information than the blurred image.

Pseudoinverse filtered image
with 59648 out of 65536 values of Hinv=0.The image is no longer intelligible.

More Difficulties with Inverse Filtering

Another problem with inverse filtering is that it dosn't perform well when used on noisy images.
We revise the "observed image" model by including an error term:

gn(x,y)=g(x,y)+n(x,y)

the corresponding 2 dimensional FFT representation would be:

Gn(u,v)=G(u,v) + N(u,v).

Inverse filtering Gn:

Fn(u,v)=Gn(u,v)/H(u,v)=F(u,v)+N(u,v)/H(u,v).

Since N(u,v) can remain large when H(u,v) gets small, it is possible that the noise will dominate Fn and consequently also fn. In such instances, it may be necessary to use a different type of filter implementation.

The following images compare the results of the restoration of a blurred image with the addition of normally distributed random noise with a mean of 0.0 and a variance of 1.

Blurred image with added noise

Image of noise was added to the image.

Results of pseudoinverse filtering with threshold value set to minimize frequency content losses.

Results using a Wiener filter implementation.

More detail of the original image can be determined from the processed image using the Wiener filter implementation. The Wiener filter attempts to minimize the error due to noise.

Topics

Refrences

Lords of Infinity Mailbox