IFFT of FFT-image results in weird upside-down overlap of original image

Question

In the following code I have a function which returns an image cropped to a centered circle of a certain radius. Then I perform fourier-tranformation of an image, and reverse fourier-transformation again, to restore the image, which works fine.

Then I have calculated, that a centered circle with radius of 43 of the energy spectrum (excluded here) will yield 99% of the energy of the original image. But when I try to apply my function "imgRadius" to the amplitude-spectrum (the fourier-transformed image), and then perform inverse-fourier-tranformation on that image, I get this weird upside-down overlap of the original image.

def imgRadius(img, radius):
    result = np.zeros(img.shape,np.float64)
    centerX = (img.shape[0])/2
    centerY = (img.shape[1])/2
    for m in range(img.shape[0]):
        for n in range(img.shape[1]):
            if math.sqrt((m-centerX)**2+(n-centerY)**2) < radius:
                result[m,n] = img[m,n]
    return result

imgpath = directory+"NorbertWiener.JPG"
img = cv.imread(imgpath, cv.IMREAD_GRAYSCALE)
norm_image = cv.normalize(img, None, alpha=0, beta=1, norm_type=cv.NORM_MINMAX, dtype=cv.CV_32F)
amp = (np.fft.fftshift(np.fft.fft2(norm_image)))
amp_log = np.log(np.abs(amp))
norm_amp = cv.normalize(amp_log, None, alpha=0, beta=1, norm_type=cv.NORM_MINMAX, dtype=cv.CV_32F)
restoredAMP = np.abs(np.fft.ifft2(np.fft.ifftshift(amp)))
radamp = imgRadius(norm_amp,43)
rest99 = np.abs(np.fft.ifft2((imgRadius(amp,43))))

plt.subplot(231),plt.imshow(norm_image, "gray", vmin=0, vmax=1),plt.title('Image')
plt.xticks([]), plt.yticks([])
plt.subplot(232),plt.imshow(norm_amp, "gray", vmin=0, vmax=1),plt.title('Amplitude')
plt.xticks([]), plt.yticks([])
plt.subplot(233),plt.imshow(restoredAMP, "gray", vmin=0, vmax=1),plt.title('Restored from amplitude')
plt.xticks([]), plt.yticks([])
plt.subplot(235),plt.imshow(rest99, "gray", vmin=0, vmax=1),plt.title('Restored from r=43')
plt.xticks([]), plt.yticks([])
plt.subplot(234),plt.imshow(radamp, "gray", vmin=0, vmax=1),plt.title('Amplitude r=43')
plt.xticks([]), plt.yticks([])
plt.show()

Does it somehow have anything to do with the fact that I use absolute values? I don't see how I should be able to plot the image without getting rid of the imaginary parts of the images, but I can see how maybe some information gets lost in the inverse fft.

I just don't get why I don't get problems when performing IFFT on the original amplitude spectrum.

Answer 1

Discarding the imaginary part of the FFT (and also the sign of the real part) is exactly what the problem is leading to the 'backfolding' of the inverted image into itself. The FFT of a function that is symmetric about its origin is real (i.e. imaginary part 0). By discarding the imaginary part, the image has thus been somehow 'symmetrized'.

After performing operations on the complex FFT result in Fourier space, one can take the inverse FFT, and then only plot the real part np.fft.ifft2((imgRadius(amp,43))).real of it.

Answer 2

The problem happens here:

def imgRadius(img, radius):
    result = np.zeros(img.shape,np.float64)

You are creating a real-valued array, and copying over the complex values. Likely either the real component or the magnitude is written to the array. In any case, the complex-valued frequency domain data becomes real-valued. The inverse transform is a symmetric matrix.

To solve the problem, initialize result as a complex-valued array.

After this, make sure to use the real component of the inverse transform, not the magnitude, as Gianluca already suggested in their answer.

Answer 3

The following works for me in Python/OpenCV/Numpy and shows the difference between using a sharp boundary circle and one that has been smoothed by Gaussian blurring in order to mitigate ringing artifacts

Input:

import numpy as np
import cv2

# read input and convert to grayscale
img = cv2.imread('lena_gray.png', cv2.IMREAD_GRAYSCALE)

# do dft saving as complex output
dft = np.fft.fft2(img)

# apply shift of origin to center of image
dft_shift = np.fft.fftshift(dft)

# generate spectrum from magnitude image (for viewing only)
mag = np.abs(dft_shift)
spec = np.log(mag) / 20

# create circle mask
radius = 32
mask = np.zeros_like(img)
cy = mask.shape[0] // 2
cx = mask.shape[1] // 2
cv2.circle(mask, (cx,cy), radius, (255,255,255), -1)[0]

# blur the mask
mask2 = cv2.GaussianBlur(mask, (19,19), 0)

# apply mask to dft_shift
dft_shift_masked = np.multiply(dft_shift,mask) / 255
dft_shift_masked2 = np.multiply(dft_shift,mask2) / 255


# shift origin from center to upper left corner
back_ishift = np.fft.ifftshift(dft_shift)
back_ishift_masked = np.fft.ifftshift(dft_shift_masked)
back_ishift_masked2 = np.fft.ifftshift(dft_shift_masked2)


# do idft saving as complex output
img_back = np.fft.ifft2(back_ishift)
img_filtered = np.fft.ifft2(back_ishift_masked)
img_filtered2 = np.fft.ifft2(back_ishift_masked2)

# combine complex components to form original image again
img_back = np.abs(img_back).clip(0,255).astype(np.uint8)
img_filtered = np.abs(img_filtered).clip(0,255).astype(np.uint8)
img_filtered2 = np.abs(img_filtered2).clip(0,255).astype(np.uint8)


cv2.imshow("ORIGINAL", img)
cv2.imshow("SPECTRUM", spec)
cv2.imshow("MASK", mask)
cv2.imshow("MASK2", mask2)
cv2.imshow("ORIGINAL DFT/IFT ROUND TRIP", img_back)
cv2.imshow("FILTERED DFT/IFT ROUND TRIP", img_filtered)
cv2.imshow("FILTERED2 DFT/IFT ROUND TRIP", img_filtered2)
cv2.waitKey(0)
cv2.destroyAllWindows()

# write result to disk
cv2.imwrite("lena_dft_numpy_mask.png", mask)
cv2.imwrite("lena_dft_numpy_mask_blurred.png", mask2)
cv2.imwrite("lena_dft_numpy_roundtrip.png", img_back)
cv2.imwrite("lena_dft_numpy_lowpass_filtered1.png", img_filtered)
cv2.imwrite("lena_dft_numpy_lowpass_filtered2.png", img_filtered2)

Sharp Mask:

Blurred Mask:

Simple Round Trip:

Low Pass Filtered Result from sharp mask (ringing obvious):

Low Pass Filtered Result from blurred mask (ringing mitigated):