Reproducing the phase spectrum while using np.fft.fft2 and cv2.dft. Why are the results not similar?

106 views Asked by At

Another question was asking about the correct way of getting magnitude and phase spectra while using cv2.dft.

My answer was limited to the numpy approach and then I thought that using OpenCV for this would be even nicer. I am currently trying to reproduce the same results but I am seeing significant differences in the phase spectrum.

Here are my imports:

%matplotlib notebook
import matplotlib.pyplot as plt
import numpy as np
import cv2
im = np.zeros((50, 50), dtype = np.float32) # create empty array
im[2:10, 2:10] = 255 # draw a rectangle

The numpy example and results:


imFFTNumpy = np.fft.fft2(im)
imFFTNumpyShifted = np.fft.fftshift(imFFTNumpy)
magSpectrumNumpy = np.abs(imFFTNumpyShifted)
phaseSpectrumNumpy = np.angle(imFFTNumpyShifted)
fig, ax = plt.subplots(nrows = 1, ncols = 3)
ax[0].imshow(im)
ax[1].imshow(magSpectrumNumpy)
ax[2].imshow(phaseSpectrumNumpy)
plt.suptitle("Using Numpy np.fft.fft2 and np.abs/ np.angle")

Numpy results

The OpenCV example and results:

imFFTOpenCV = cv2.dft(im, flags=cv2.DFT_COMPLEX_OUTPUT)
imFFTOpenCVShifted = np.fft.fftshift(imFFTOpenCV)
magSpectrumOpenCV, phaseSpectrumOpenCV = cv2.cartToPolar(imFFTOpenCVShifted[:,:,0], imFFTOpenCVShifted[:,:,1])
fig, ax = plt.subplots(nrows = 1, ncols = 3)
ax[0].imshow(im)
ax[1].imshow(magSpectrumOpenCV)
ax[2].imshow(phaseSpectrumOpenCV)
plt.suptitle("Using OpenCV cv2.dft and cv2.cartToPolar")

OpenCV results

As you can see, while the magnitude spectrum looks the same (it has some expected deviations due to floating-point arithmetic), the phase spectrum looks significantly different. I dug around a bit and found out that OpenCV usually returns phase from 0 to 2π, whereas np.angle returns the phase from -π to +π. Subtracting π from the OpenCV phase does not correct difference though.

What could be the reason for this? Is it possible to get almost identical phase using both approaches, just like with magnitude?

2

There are 2 answers

4
Cris Luengo On BEST ANSWER

Given that imFFTOpenCV is a 3D array because OpenCV doesn’t understand complex numbers, np.fft.fftshift(imFFTOpenCV) will swap the real and complex planes. That is, the shift happens in all 3 dimensions of the array.

So when computing the phase and magnitude, you need to take this swap into account:

magSpectrumOpenCV, phaseSpectrumOpenCV = cv2.cartToPolar(imFFTOpenCVShifted[:,:,1], imFFTOpenCVShifted[:,:,0])

Alternatively, and this is probably more readable, you could tell NumPy which axes you want shifted:

imFFTOpenCVShifted = np.fft.fftshift(imFFTOpenCV, axes=(0, 1))
1
cards On

The main issue is that numpy.angle and cv2.phase are not comparable because their underline implementation is different.

  • numpy.angle works with a fix domain, (-π, π]
  • cv2.phase uses arctan2 which is a piecewise defined function and its output depends on the possible sign of the rise and the run. That's the reason why in OP is written "OpenCV usually returns phase from 0 to 2π"

To be able to compare the phase of a complex number with the two libraries is much easier to forget about the numpy.angle and just use numpy.arctan2:

im = np.array([[1, 0, 1], [1, 0, 1], [1, 0, 1]], dtype=np.uint8)

# numpy
fft_np = np.fft.fft2(im)
phase_np = np.mod(np.arctan2(np.real(fft_np), np.imag(fft_np)) * 180/np.pi, 360)

# opencv
fft_cv = cv2.dft(np.float32(im), flags=cv2.DFT_COMPLEX_OUTPUT)
phase_cv = cv2.phase(fft_cv[:,:,1], fft_cv[:,:,0], angleInDegrees=True)

# display results
print(phase_np)
print(phase_cv)

Output

[[ 90.  30. 150.]
 [  0.   0.   0.]
 [  0.   0.   0.]]

[[ 90.        29.997692 150.0023  ]
 [  0.         0.         0.      ]
 [  0.         0.         0.      ]]

But not think to be finished yet! numpy.arctan2 seems to return also negative angles but cv2.phase is always positive... so you need to use a further modulo operation of a full period, so either 360 or 2π, on its result. This can be done with numpy.mod.

You can check it out by using [[1, 0, 1, 1], [0, 1, 1, 0], [1, 0, 1, 1], [0, 0, 1, 0]] in the above example.


Here, the final numpy-counter part of cv2.phase:

def phase_np(x, y, angleInDegrees=False): # angle_np2cv
    if angleInDegrees:
        return np.mod(np.arctan2(x, y) * 180/np.pi, 360)
    return np.mod(np.arctan2(x, y), 2*np.pi)

Also another difference could be the datatypes: np.fft2 uses complex128 and computations (np.angle, np.arctan2, ...) on it are usually float64.


Conclusion: cv2.cartToPolar is arctan2-based as cv2.phase so the same holds.