Demo iainversefiltering

Illustrate the inverse filtering for restoration.

Description

The Inverse Filtering is a restoration method that tries to recover an image that has been corrupted by a known linear translation invariant (LTI) filter. The experiment shown here is only dicdatical and consists of corrupt an image by a known LTI filter and then apply the inverse of this filter to recover the image back to its original.

Script

Original and distorted images

The original image is corrupted using a low pass Butterworth filter with cutoff period of 8 pixels and order 4. This has the a similar effect of an out of focus image.

1. from ia636 import iabwlp, ianormalize
2. import numpy as np
3. 
4. f = iaread('keyb.pgm').astype(float)
5. F = np.fft.fft2(f)                # Discrete Fourier Transform of f
6. H = iabwlp(f.shape, 16, 4)      # Butterworth filter, cutoff period 16, order 4
7. G = F*H                         # Filtering in frequency domain
8. g = np.fft.ifft2(G).real   # inv DFT
9. adshow(ianormalize(g, [0,255])) # display distorted image

Spectrum of the blurred image and its inverse filter

1. from ia636 import iadftview
2. 
3. adshow(iadftview(G))  # Spectrum of the corrupted image
4. IH = 1.0/H            # The inv filter
5. adshow(iadftview(IH)) # Spectrum of the inv filter

Applying the inverse filter

1. FR = G*IH                              # Inverse filtering
2. adshow(iadftview(FR))                  # Spectrum of the restored image
3. fr = np.fft.ifft2(FR).real
4. adshow(ianormalize(fr).astype(uint8))  # display the restored image

Adding a little of noise

The previous example is rather didactical. Just as an experience, instead of using the corrupted image with pixels values represented in double, we will use a integer truncation to the pixel values.

1. gfix = floor(g)                        # truncate pixel values
2. Gfix = np.fft.fft2(gfix)               # DFT of truncated filted image
3. FRfix = Gfix * IH                      # applying the inv filtering
4. frfix = np.fft.ifft2(FRfix).real
5. adshow(gfix)
6. adshow(ianormalize(frfix, [0,255]))    # display the restored image
7. adshow(iadftview(FRfix))               # spectrum of the restored image

Why is the result so noisy?

When the distorted image is rounded, it is equivalent of subtracting a noise from the distorted image. This noise has uniform distribution with mean 0.5 pixel intensity. The inverted filter is very high pass frequency, and at high frequency, the noise outcomes the power of the signal. The result is a magnification of high frequency noise.

1. fn = g - gfix                         # noise that was added
2. print [fn.mean(), fn.min(), fn.max()] # mean, minimum and maximum values
3. adshow(ianormalize(fn, [0,255]))      # display noise
4. FN = np.fft.fft2(fn)
5. adshow(iadftview(FN))                 # spectrum of the noise
6. FNI = FN*IH                           # inv filtering in the noise
7. fni = np.fft.ifft2(FNI).real
8. print [fni.min(), fni.max()]          # min and max of restored noise
9. adshow(ianormalize(fni, [0,255]))     # display restoration of noise
[0.49940999403058939, 9.2963516067356977e-06, 0.99996062920030226]
[-15348654.784229999, 15787958.430654926]