I have a problem with getting the right phase information out of the 2D fast-fourier transform (2D-FFT) with scipy.
I have already found similar questions/solutions, which have been insightful, but have not been yielding the right answer for me.
In the following I will showcase my code and my thought behind it and I hope somebody can help me.
My code is split up on multiple files, but here are the in my opinion relevant sections from it that can be combined in a file and executed to show the output, for readability's sake I left them split up
First of all I set up a 2D-radial grid with a pixel scaling that is unit/px:
import numpy as np
from typing import Any, Dict, List, Optional, Union
def set_size(mas_size: int, size: int,
incline_params: List[float] = None,
sampling: Optional[int] = None) -> np.array:
"""Sets the size of the model and returns a 2D-radial grid
Parameters
----------
mas_size: int
Sets the size of the images [mas] and implicitly the pixels
size: int
Sets the range of the model grid and implicitly the x-, y-axis.
incline_params: List[float]
A list of the inclination parameters [axis_ratio, pos_angle] [mas, rad]
sampling: int, optional
The pixel sampling
Returns
-------
radius: np.array
The 2D-radius grid
"""
with np.errstate(divide='ignore'):
fov_scale = mas_size/size
if sampling is None:
sampling = size
x = np.linspace(-size//2, size//2, sampling, endpoint=False)*fov_scale
y = x[:, np.newaxis]
axis_ratio, pos_angle = incline_params[0], np.radians(incline_params[1])
xr, yr = -x*np.cos(pos_angle)+y*np.sin(pos_angle), \
(x*np.sin(pos_angle)+y*np.cos(pos_angle))/axis_ratio
return np.sqrt(xr**2+yr**2), fov_scale
Then I plug this into a function for a ring, for example
def ring(axis_ratio: float, pos_angle: int, mas_size: int, px_size: int,
sampling: Optional[int] = None,
inner_radius: Optional[int] = None,
outer_radius: Optional[int] = None) -> np.array:
"""Evaluates the model
Parameters
----------
axis_ratio: float, optional
The ratio that determines ellipse (>1.) or ring (=1.)
pos_angle: int | float, optional
The angle of rotation of the model image
mas_size: int
Sets the size of the images [mas] and implicitly the pixels
px_size: int
The size of the model image
sampling: int, optional
The sampling of the object-plane
inner_radius: int, optional
A set inner radius overwriting the sublimation radius
Returns
--------
model: np.array
"""
radius, fov_scale = set_size(mas_size, px_size,[axis_ratio, pos_angle], sampling)
if inner_radius:
radius[radius < inner_radius] = 0.
else:
radius[radius < 1.] = 0.
if outer_radius:
radius[radius > outer_radius] = 0.
if inner_radius:
radius[np.where(radius != 0)] = 1/(2*np.pi*inner_radius)
else:
radius[np.where(radius != 0)] = 1/(2*np.pi)
return radius, fov_scale
And I then do the FFT of this in the following way
import matplotlib.pyplot as plt
from scipy.fft import ifftshift, fftshift, fftfreq, fft2
model, fov_scale = ring(1.5, 135, 10, 2**9, inner_radius=1., outer_radius=1.1)
model_centre = model.shape[0]//2
# Zero padding-order of 3
zero_padding = 2**int(np.log2(model.shape[0]) + 3)
freq_axis = fftshift(fftfreq(zero_padding, fov_scale))
# Scale by the wavelength of your model
freq_scaling_m = np.diff(freq_axis)[0]*8e-6
padded_image = np.zeros((zero_padding, zero_padding))
pad_centre = padded_image.shape[0]//2
mod_min, mod_max = pad_centre-model_centre, pad_centre+model_centre
padded_image[mod_min:mod_max, mod_min:mod_max] = model
ft = fftshift(fft2(ifftshift(padded_image)))
amp, phase = abs(ft), np.angle(ft, deg=True)
fig, (ax, bx, cx) = plt.subplots(1, 3)
ax.imshow(model)
bx.imshow(amp)
cx.imshow(phase)
plt.show()
And finally, here the plotted output images with the left being the object plane image (the model), the middle one being the amplitude of the FFT (normed) and the right most one being the phase space of the FFT:
2D-FFT of Ring
If I am not completely mistaken then the phase information should have circles that have either +180° or -180° and should not have mixed degrees circles...
Is there any glaring error here or do I have a wrong mindset about what the phase should look like for the FFT?
Thank you everyone for your help already!
PS: An additional question I have would be, does anyone know where there are good repositories or similar that show 2D-FFTs of various images (e.g., a Ring, Gaussian, etc.) or at least what the phases of a 2D-FFT of a 2D-Gaussian should look like?
