How to inverse a DFT with magnitude with opencv python

Question

I'm new to all of this, I would like to get a magnitude spectrum from an image and then rebuild the image from a modified magnitude spectrum.. But for now i'am getting a very dark reconstitution.

import numpy as np
import cv2
from matplotlib import pyplot as plt

img = cv2.imread('IMG.jpg',0)

dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

m, a = np.log(cv2.cartToPolar(dft_shift[:,:,0],dft_shift[:,:,1]))

# do somthing with m

x, y = cv2.polarToCart(np.exp(m), a)


back = cv2.merge([x, y])


f_ishift = np.fft.ifftshift(back)
img_back = cv2.idft(f_ishift)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])

plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(m, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back, cmap = 'gray')
plt.title('result'), plt.xticks([]), plt.yticks([])
plt.show()

the result

Can you guys help me figure out why is this so dark.

Thank in advance :)

EDIT

I tryed to normalise the image, but it's not working. I'm still having a very dark image.

import numpy as np
import cv2
from matplotlib import pyplot as plt

img = cv2.imread('IMG.jpg',0)

dft = cv2.dft(np.float32(img),flags = cv2.DFT_COMPLEX_OUTPUT)
dft_shift = np.fft.fftshift(dft)

m, a = np.log1p(cv2.cartToPolar(dft_shift[:,:,0],dft_shift[:,:,1]))

# modify m, then use the modify m to reconstruct


x, y = cv2.polarToCart(np.expm1(m), a)


back = cv2.merge([x, y])


f_ishift = np.fft.ifftshift(back)
img_back = cv2.idft(f_ishift, flags=cv2.DFT_SCALE)
img_back = cv2.magnitude(img_back[:,:,0],img_back[:,:,1])


min, max = np.amin(img, (0,1)), np.amax(img, (0,1))
print(min,max)

# re-normalize to 8-bits
min, max = np.amin(img_back, (0,1)), np.amax(img_back, (0,1))
print(min,max)
img_back = cv2.normalize(img_back, None, alpha=0, beta=252, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)


min, max = np.amin(img_back, (0,1)), np.amax(img_back, (0,1))
print(min,max)


plt.subplot(131),plt.imshow(img, cmap = 'gray')
plt.title('Input Image'), plt.xticks([]), plt.yticks([])
plt.subplot(132),plt.imshow(m, cmap = 'gray')
plt.title('Magnitude Spectrum'), plt.xticks([]), plt.yticks([])
plt.subplot(133),plt.imshow(img_back, cmap = 'gray')
plt.title('result'), plt.xticks([]), plt.yticks([])
plt.show()
cv2.waitKey(0)
cv2.destroyAllWindows()

output:

0 252
0.36347726 5867.449
0 252

I would like to modify the magnitude spectrum and used the modify version to reconstruct the image.

Gorlin

Did you have a look at the Opencv DFT documentation? This might help you. You see that idft(dft(X))=X/size(X) so you have an attenuation factor – Ivanaivanah

You can try to narrow down the issue by removing all your magnitude modifications and just keep the dft and idft to check the output of these functions – Ivanaivanah

Hi, i have try the DFT_SCALE flag but with no sucess. If y try to pass the result of the dft to the idft, it's working fine. – Gorlin

You should extract the magnitude and phase images from cartToPolar() without applying the log. Then do the log separately only to view the spectrum, keeping the magnitude in its original form. The modify the original magnitude as desired before doing the inverse dft. What do you want to do in the way of modifications? Please post individual images, so others can easily start with your input. – Sublittoral

You might want to check the min, max and type of your input (img) and output (img_back) – Ivanaivanah

Answer 1

If you need to modify the magnitude by raising it to a power near 1 (called coefficient rooting or alpha rooting), then it is just a simple modification of my code above using Python/OpenCV. Simply add cv2.pow(mag, 1.1) before converting the magnitude and phase back to real and imaginary components.

Input:

import numpy as np
import cv2

# read input as grayscale
img = cv2.imread('lena.png', 0)

# convert image to floats and do dft saving as complex output
dft = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)

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

# extract magnitude and phase images
mag, phase = cv2.cartToPolar(dft_shift[:,:,0], dft_shift[:,:,1])

# get spectrum for viewing only
spec = np.log(mag) / 30

# NEW CODE HERE: raise mag to some power near 1
# values larger than 1 increase contrast; values smaller than 1 decrease contrast
mag = cv2.pow(mag, 1.1)

# convert magnitude and phase into cartesian real and imaginary components
real, imag = cv2.polarToCart(mag, phase)

# combine cartesian components into one complex image
back = cv2.merge([real, imag])

# shift origin from center to upper left corner
back_ishift = np.fft.ifftshift(back)

# do idft saving as complex output
img_back = cv2.idft(back_ishift)

# combine complex components into original image again
img_back = cv2.magnitude(img_back[:,:,0], img_back[:,:,1])

# re-normalize to 8-bits
min, max = np.amin(img_back, (0,1)), np.amax(img_back, (0,1))
print(min,max)
img_back = cv2.normalize(img_back, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)

cv2.imshow("ORIGINAL", img)
cv2.imshow("MAG", mag)
cv2.imshow("PHASE", phase)
cv2.imshow("SPECTRUM", spec)
cv2.imshow("REAL", real)
cv2.imshow("IMAGINARY", imag)
cv2.imshow("COEF ROOT", img_back)
cv2.waitKey(0)
cv2.destroyAllWindows()

# write result to disk
cv2.imwrite("lena_grayscale_opencv.png", img)
cv2.imwrite("lena_grayscale_coefroot_opencv.png", img_back)

Original Grayscale:

Coefficient Rooting Result:

Here is an animation showing the differences (created using ImageMagick):

Sublittoral

I find the solution!! I had np.log the magnitude and the phase but i did not np.exp the phase back to normal. Thanks for the help – Gorlin

Answer 2

You should extract the magnitude and phase images from cartToPolar() without applying the log. Then do the log separately only to view the spectrum, keeping the magnitude in its original form. Then modify the original magnitude as desired before doing the inverse dft.

One of the other issues is that the round trip image needs to be rescaled back to 8-bit range and datatype. I do that with cv2.normalize(). You can see that need from the printed min and max values.

Here is how to do the dft, get the spectrum and then do the inverse dft in Python/OpenCV. I start with a color image, but I convert that to grayscale when reading it in. The final returned round trip dft/idft will still be grayscale.

Input:

import numpy as np
import cv2

# read input as grayscale
img = cv2.imread('lena.png', 0)

# convert image to floats and do dft saving as complex output
dft = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)

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

# extract magnitude and phase images
mag, phase = cv2.cartToPolar(dft_shift[:,:,0], dft_shift[:,:,1])

# get spectrum for viewing only
spec = np.log(mag) / 30

# convert magnitude and phase into cartesian real and imaginary components
real, imag = cv2.polarToCart(mag, phase)

# combine cartesian components into one complex image
back = cv2.merge([real, imag])

# shift origin from center to upper left corner
back_ishift = np.fft.ifftshift(back)

# do idft saving as complex output
img_back = cv2.idft(back_ishift)

# combine complex components into original image again
img_back = cv2.magnitude(img_back[:,:,0], img_back[:,:,1])

# re-normalize to 8-bits
min, max = np.amin(img_back, (0,1)), np.amax(img_back, (0,1))
print(min,max)
img_back = cv2.normalize(img_back, None, alpha=0, beta=255, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)

cv2.imshow("ORIGINAL", img)
cv2.imshow("MAG", mag)
cv2.imshow("PHASE", phase)
cv2.imshow("SPECTRUM", spec)
cv2.imshow("REAL", real)
cv2.imshow("IMAGINARY", imag)
cv2.imshow("ORIGINAL DFT/IFT ROUND TRIP", img_back)
cv2.waitKey(0)
cv2.destroyAllWindows()

# write result to disk
cv2.imwrite("lena_dft_ift_opencv.png", img_back)

Result:

Sublittoral

Thanks for the normalisation tips, but I would like to modify the magnitude spectrum and used the modify version to reconstruct the image. I have edited my post. – Gorlin

Answer 3

Here is how to remove repetitive patterned noise from an image using notch filtering in Fourier domain using Python/OpenCV

• Read image
• Do DFT
• Generate the magnitude and phase components from the real and imaginary ones
• Create spectrum From magnitude
• Threshold the spectrum image to make a mask covering over the DC region in the center of the thresholded image with black
• Apply the mask to the magnitude
• Combine the new magnitude and the original phase
• Convert those to real and imaginary components
• Do IDFT
• Save the result

Input with repetitive patterned noise:

import numpy as np
import cv2

# read input as grayscale
img = cv2.imread('clown.jpg', 0)

# get min and max values of img
img_min, img_max = np.amin(img, (0,1)), np.amax(img, (0,1))
print(img_min,img_max)

# convert image to floats and do dft saving as complex output
dft = cv2.dft(np.float32(img), flags = cv2.DFT_COMPLEX_OUTPUT)

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

# extract magnitude and phase images
mag, phase = cv2.cartToPolar(dft_shift[:,:,0], dft_shift[:,:,1])

# get spectrum
spec = np.log(mag) / 20

# create mask from spectrum keeping only the brightest spots as the notches
mask = cv2.normalize(spec, None, alpha=0, beta=1, norm_type=cv2.NORM_MINMAX)
mask = cv2.threshold(mask, 0.65, 1, cv2.THRESH_BINARY)[1]

# dilate mask
kernel = cv2.getStructuringElement(cv2.MORPH_ELLIPSE, (3,3))
mask = cv2.morphologyEx(mask, cv2.MORPH_DILATE, kernel)

# cover center DC component by circle of black leaving only a few white spots on black background
xcenter = mask.shape[1] // 2
ycenter = mask.shape[0] // 2
mask = cv2.circle(mask, (xcenter,ycenter), radius=10, color=0, thickness=cv2.FILLED)

# apply mask to magnitude such that magnitude is made zero where mask is one, ie at spots
mag[mask!=0] = 0

# convert new magnitude and old phase into cartesian real and imaginary components
real, imag = cv2.polarToCart(mag, phase)

# combine cartesian components into one complex image
back = cv2.merge([real, imag])

# shift origin from center to upper left corner
back_ishift = np.fft.ifftshift(back)

# do idft saving as complex output
img_back = cv2.idft(back_ishift)

# combine complex components into original image again
img_back = cv2.magnitude(img_back[:,:,0], img_back[:,:,1])

# re-normalize to 8-bits in range of original
min, max = np.amin(img_back, (0,1)), np.amax(img_back, (0,1))
print(min,max)
notched = cv2.normalize(img_back, None, alpha=img_min, beta=img_max, norm_type=cv2.NORM_MINMAX, dtype=cv2.CV_8U)

cv2.imshow("ORIGINAL", img)
cv2.imshow("MAG", mag)
cv2.imshow("PHASE", phase)
cv2.imshow("SPECTRUM", spec)
cv2.imshow("MASK", mask)
cv2.imshow("NOTCHED", notched)
cv2.waitKey(0)
cv2.destroyAllWindows()

# write result to disk
cv2.imwrite("clown_mask.png", (255*mask).clip(0,255).astype(np.uint8))
cv2.imwrite("clown_notched.png", notched)

Spectrum:

Mask:

Notch Filtered Result (noise removed):

Animation (created separately using Imagemagick):

Sublittoral