Scipy 2D-FFT计算二维环形模型相位结果异常问题

问题概述 #

使用Scipy库开展二维快速傅里叶变换（2D-FFT）计算时，无法提取到符合预期的相位信息，已查阅的同类问题解决方案虽具备参考价值，但未能适配当前场景得到正确结果。
当前认知为：环形模型的2D-FFT相位应当呈现整环统一为+180°或-180°的分布，不应出现度数混杂的环形结构。需要排查代码错误、确认相位认知是否存在偏差，同时获取典型图像2D-FFT结果的参考资料。

实现代码 #

1. 二维径向网格生成函数set_size #

功能为生成像素缩放为unit/px的二维径向网格，支持配置mas级视场尺寸、网格范围、倾角参数（轴比、位置角，单位分别为mas、rad）、像素采样率，返回二维半径网格与视场缩放系数。

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

2. 环形/椭圆环模型生成函数ring #

功能为生成环形/椭圆环模型，支持配置轴比（=1为正圆环，>1为椭圆环）、位置角、mas级视场尺寸、像素尺寸、物平面采样率、内外半径参数，返回模型数组与视场缩放系数。

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

3. FFT计算与绘图流程 #

执行逻辑如下：

• 调用接口生成轴比1.5、位置角135°、内半径1.0mas、外半径1.1mas的环形模型
• 执行3阶零填充
• 基于scipy.fft模块的ifftshift、fft2、fftshift接口完成FFT计算
• 分别提取振幅、相位结果后通过matplotlib绘图

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()

运行结果 #

代码运行后输出三幅结果图：左图为物平面环形模型，中图为FFT振幅结果，右图为FFT相位结果。

2D-FFT of Ring

待解答疑问 #

• 现有代码是否存在明显错误，导致相位结果不符合预期
• 环形模型2D-FFT相位应当呈现整环统一为+180°或-180°分布的认知是否存在偏差
• 是否有优质公开代码仓库或参考资料展示各类典型图像（环形、二维高斯等）的2D-FFT结果，尤其是二维高斯模型的2D-FFT相位分布特征

内容的提问来源于stack exchange，提问作者MBSck