Demo 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.
>>> import Numeric, FFT
>>> f = iaread('keyb.pgm').astype(Numeric.Float)
>>> F = FFT.fft2d(f) # Discrete Fourier Transform of f
>>> H = iabwlp(f.shape, 16, 4) # Butterworth filter, cutoff period 16, order 4
>>> G = F*H # Filtering in frequency domain
>>> g = FFT.inverse_fft2d(G).real # inverse DFT
>>> iashow(ianormalize(g, [0,255])) # display distorted image
(316, 295) Min= 0.0 Max= 255.0 Mean=127.038 Std=57.84
 |
|
| ianormalize(g, [0,255]) |
Spectrum of the blurred image and its inverse filter
>>> iashow(iadftview(G)) # Spectrum of the corrupted image
(316, 295) Min= 0 Max= 254 Mean=9.676 Std=29.71
>>> IH = 1.0/H # The inverse filter
>>> iashow(iadftview(IH)) # Spectrum of the inverse filter
(316, 295) Min= 9 Max= 255 Mean=177.053 Std=51.93
 |
 |
|
| iadftview(G) | iadftview(IH) |
Applying the inverse filter
>>> FR = G*IH # Inverse filtering
>>> iashow(iadftview(FR)) # Spectrum of the restored image
(316, 295) Min= 35 Max= 254 Mean=118.542 Std=17.79
>>> fr = FFT.inverse_fft2d(FR).real
>>> iashow(fr) # display the restored image
(316, 295) Min= -4.1157513403e-013 Max= 255.0 Mean=123.284 Std=61.65
 |
 |
|
| iadftview(FR) | fr |
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.
>>> gfix = Numeric.floor(g) # truncate pixel values
>>> Gfix = FFT.fft2d(gfix) # DFT of truncated filted image
>>> FRfix = Gfix * IH # applying the inverse filtering
>>> frfix = FFT.inverse_fft2d(FRfix).real
>>> iashow(gfix)
(316, 295) Min= 4.0 Max= 242.0 Mean=122.785 Std=54.08
>>> iashow(ianormalize(frfix, [0,255])) # display the restored image
(316, 295) Min= 0.0 Max= 255.0 Mean=129.299 Std=30.88
>>> iashow(iadftview(FRfix)) # spectrum of the restored image
(316, 295) Min= 58 Max= 255 Mean=186.088 Std=36.25
 |
 |
|
| gfix | ianormalize(frfix, [0,255]) |
 |
|
| iadftview(FRfix) |
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.
>>> import MLab
>>> fn = g - gfix # noise that was added
>>> print [MLab.mean(Numeric.ravel(fn)), MLab.min(Numeric.ravel(fn)), MLab.max(Numeric.ravel(fn))] # mean, minumum and maximum values
[0.4994099940305825, 9.2963516067356977e-006, 0.99996062920030226]
>>> iashow(ianormalize(fn, [0,255])) # display noise
(316, 295) Min= 0.0 Max= 255.0 Mean=127.353 Std=73.42
>>> FN = FFT.fft2d(fn)
>>> iashow(iadftview(FN)) # spectrum of the noise
(316, 295) Min= 3 Max= 255 Mean=99.354 Std=14.72
>>> FNI = FN*IH # inverse filtering in the noise
>>> fni = FFT.inverse_fft2d(FNI).real
>>> print [MLab.min(Numeric.ravel(fni)), MLab.max(Numeric.ravel(fni))] # min and max of restored noise
[-15348654.784230139, 15787958.430655491]
>>> iashow(ianormalize(fni, [0,255])) # display restoration of noise
(316, 295) Min= 0.0 Max= 255.0 Mean=125.701 Std=30.88
 |
 |
|
| ianormalize(fn, [0,255]) | iadftview(FN) |
 |
|
| ianormalize(fni, [0,255]) |
