CODEWITHSUNDEEP
pythonimage-processingnumpy
David Poole
David Poole

Reputation: 3510

Am I using Numpy to calculate the Inverse Filter correctly?

As part of a digital image processing class, we have been assigned the Inverse Filter for image restoration. I'm using numpy. The variable names below try to follow the names in Digital Image Processing Gonzalez+Woods, 3e.

A zoom of the original image. A zoom of the original image.

Gaussian kernel "zz.tif" same size as original image. Gaussian kernel "zz.tif" same size as original image

Zoom of the gaussian smoothed image with no noise added Zoom of the gaussian smoothed image with no noise added

f = imtools.load_image( sys.argv[1], mode="L", dtype="float" )
zz = imtools.load_image( "zz.tif", mode="L", dtype="float" )

F = np.fft.fft2( f )
F2 = np.fft.fftshift( F )

# normalize to [0,1]
H = zz/255.

# calculate the damaged image
G = H * F2

# Inverse Filter 
F_hat = G / H

# cheat? replace division by zero (NaN) with zeroes
a = np.nan_to_num(F_hat)
f_hat = np.fft.ifft2( np.fft.ifftshift(a) )

imtools.save_image( np.abs(f_hat), "out.tif" )

imtools is just my wrapper using PIL+numpy to load/store images. (Can post that src, too.)

Zoom of the inverse filtered image. Zoom of the inverse filtered image.

Am I calculating the Inverse Filter correctly? Am I using numpy correctly?

Is the ringing in the final image expected or am I doing something wrong?

Upvotes: 1

Views: 4148

Answers (2)

user4337532
user4337532

Reputation: 31

I don't know much about Python but the 'ringing' is normal for the inverse filter. The Gibbs phenomenon lies at the basis of the ringing. Since the input is not entirely smooth but has some discontinuities, an infinite number of Fourier components is in principle needed to represent it completely. A finite number of components is sufficient here since the display resolution is finite, the image is pixelated. However, some information is lost in the recorded image because of the multiplication by zeros in H, by consequence the restored image approximates the input image with components covering a finite bandwidth, lower than that of the display, revealing the Gibbs oscillations.

To mitigate this use proper regularization as with a 2D Wiener filter: F_hat=G * H.conjugate()/(abs(H)2+NSR2) where NSR is an estimate of the noise to signal ratio, e.g. linearly increasing from 0 to 10 at the highest spatial frequency. This will account for the finite signal to noise ratio and when the NSR estimate is close enough you should see little 'ringing' after restoration.

Upvotes: 3

Joe Kington
Joe Kington

Reputation: 284652

Generally, yes you seem to be doing things correctly, as far as I know.

The ringing is due to an overly "sharp" high pass filter, but that's what the method you're using does.

However, you might consider using numpy.fft.rfft2 ("real fft") and numpy.fft.irfft2 instead of numpy.fft.fft2 and numpy.fft.ifft2 because you're dealing purely with real values. It should be slightly faster.

Upvotes: 3

Related Questions