第14章：采样、归一化以及那些破坏PSF与MTF的小错误

至此，我们已有足够的代码来绘制令人信服的图像。

这是好消息。

但也很危险。

当采样出现错误时，PSF图看起来依然可以很平滑且科学。MTF曲线可以从1开始，优雅地下降，但其频率轴却是错的。傅里叶变换在数学上可能正确，却代表着错误的物理孔径。一个平移的数组可能产生看起来几乎正确的结果，但其中心却位于错误的位置。

本章是一章“补救”性质的章节。它不会引入新的光学概念，而是让前三个章节更为安全。

我们已经建立起这样一条链路：

光程差 (OPD)
→ 光瞳函数
→ PSF
→ OTF
→ MTF

现在，我们放慢脚步并追问：

这条链路何时才能产生物理上有意义的结果？

答案取决于那些看似细小，却会彻底破坏一切的细节：

• FFT平移（shift）；
• 阵列中心化；
• 孔径掩膜（mask）；
• 补零（zero padding）；
• 光瞳采样；
• 像平面采样；
• 频率轴单位；
• PSF归一化；
• OTF归一化；
• 波长单位；
• 焦距单位；
• 子午与弧矢方向。

这些并非装饰性的实现细节，它们决定了图像的实际含义。

本章有一条尖锐但有用的准则：

一个正确的公式救不了一个采样错误的数组。

让我们逐一梳理常见的陷阱。

1. 问题所在：漂亮但错误的图

傅里叶变换极为“顺从”。

它不会询问你的波长是以米还是纳米为单位；它不知道你的光瞳直径是5毫米还是5米；它不知道零频是否应位于阵列中心；它也不知道你的PSF是否已经归一化。

它只处理数字。

这正是为什么基于FFT的光学计算可能静悄悄地出错。

错误的光线–表面交点通常会产生明显的错误：光线飞散、交点消失，或者程序崩溃。

错误的MTF计算却可以产生一条完美平滑的曲线。

这更为危险，因为它看起来像模像样。

实践中的教训是：

在FFT光学中，调试的第一个问题不是：
“这张图看起来像光学图吗？”

而是：
“它的单位、采样间隔、归一化方式和坐标约定是什么？”

本章将为你提供一个回答此问题的核查清单。

2. 三个数组，三种坐标系

光瞳数组、PSF数组和MTF数组并不生活在同一个坐标系中。

这是混淆的首要来源。

光瞳数组

光瞳数组描述的是孔径上的场分布。

其坐标可以是：

归一化的光瞳坐标

也可以是物理坐标，例如：

光瞳上的米或毫米

光瞳函数为：

$$P(x,y)$$

其中 $x,y$ 是光瞳平面坐标。

PSF数组

PSF数组描述的是像平面上的强度分布。

其坐标可以是：

角坐标

也可以是：

焦平面上的长度坐标
$$

PSF为：

$$
\mathrm{PSF}(X,Y)
$$

其中 $X,Y$ 是像平面坐标。

### MTF数组

MTF数组描述的是空间频率响应。

其坐标是频率坐标：

```text
线对/毫米

或者：

线对/弧度

有时是归一化的频率。

MTF为：

$$ \mathrm{MTF}(f_x,f_y) $$

其中 $f_x,f_y$ 是空间频率。

因此，这条链路不仅是：

数组 → 数组 → 数组

而是：

光瞳平面场
→ 像平面强度
→ 像平面空间频率响应

这种区分可以节省大量调试时间。

3. DFT的采样关系

离散傅里叶变换将采样间隔联系在一起。

若以间隔：

$$ \Delta x $$

对一个函数进行 $N$ 点采样，则对应频率样本的间隔为：

$$ \Delta f = \frac{1}{N\Delta x} $$

对于偶数长度的数组，经过 fftshift 之后，平移后的频率坐标为：

$$ f_k = \frac{k - N/2}{N\Delta x} $$

在NumPy中，这表示为：

freq = np.fft.fftshift(np.fft.fftfreq(N, d=dx))

这一行代码是构建频率轴最安全的方式。

重要之处在于：

• N 是被变换数组的样本点数；
• dx 是输入域中的样本间隔；
• 输出单位是每 dx 单位的周期数（线对数）。

因此，若：

dx = 0.01  # 毫米

则：

freq 的单位是 线对/毫米

若：

dx = 10e-6  # 米

则：

freq 的单位是 线对/米

FFT本身不携带单位。你通过使用正确的 dx 来赋予其单位。

这是本章最重要的采样概念。

4. 从光瞳采样到PSF采样

当我们从光瞳函数计算PSF时，实际上是在对光瞳场进行傅里叶变换。

对于夫琅禾费衍射，光瞳的空间频率映射的是衍射角度。

一个有用的实践关系是：

$$ \theta = \lambda f_\text{pupil} $$

其中：

• $\theta$ 是衍射角度（弧度）；
• $\lambda$ 是波长；
• $f_\text{pupil}$ 是与光瞳平面采样相关的空间频率。

若光瞳采样间隔为 $\Delta u$，并且补零后的光瞳数组有 $N$ 个样本点，则角度采样间隔近似为：

$$ \Delta \theta = \frac{\lambda}{N\Delta u} $$

在焦距 $F$ 下，焦平面上的采样间隔为：

$$ \Delta X = F\Delta\theta $$

因此：

$$ \Delta X = \frac{F\lambda}{N\Delta u} $$

这个公式值得慢慢阅读。

它表明PSF的采样间隔取决于：

• 焦距；
• 波长；
• FFT补零后的尺寸 $N$；
• 物理光瞳的采样间隔。

补零会改变 $N$，因此会改变显示的PSF采样间隔，为PSF的显示提供更多的采样点。

但它不会改变物理孔径。

以下是一个辅助函数：

import numpy as np


def psf_sample_spacing_from_pupil(
    wavelength,
    focal_length,
    pupil_sample_spacing,
    n_fft
):
    """
    估计夫琅禾费传播下焦平面PSF的采样间隔。

    参数
    ----------
    wavelength : float
        波长，单位：米。
    focal_length : float
        焦距，单位：米。
    pupil_sample_spacing : float
        光瞳采样点之间的物理间距，单位：米。
    n_fft : int
        补零后的FFT数组尺寸。

    返回
    -------
    dx_image : float
        近似的像平面采样间隔，单位：米。
    """
    dtheta = wavelength / (n_fft * pupil_sample_spacing)
    dx_image = focal_length * dtheta

    return dx_image

如果你希望结果以毫米为单位：

dx_image_mm = dx_image * 1e3

这是从光瞳数组通往物理PSF坐标的桥梁。

5. 示例：具有物理坐标标记的PSF

让我们创建一个简单的物理设置。

假设：

wavelength = 550 nm
focal length = 50 mm
pupil diameter = 10 mm

F数为：

$$ N = \frac{F}{D} = \frac{50}{10} = 5 $$

现在构建光瞳。

import numpy as np
import matplotlib.pyplot as plt


def make_pupil_grid(n=256):
    coord = np.linspace(-1.0, 1.0, n, endpoint=False)
    x, y = np.meshgrid(coord, coord)

    rho = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    mask = rho <= 1.0

    return x, y, rho, theta, mask


def pupil_function_from_opd(opd, mask, wavelength, amplitude=None):
    if amplitude is None:
        amplitude = mask.astype(float)

    phase = 2.0 * np.pi * opd / wavelength
    pupil = amplitude * np.exp(1j * phase)
    pupil[~mask] = 0.0 + 0.0j

    return pupil, phase


def pad_array_centered(array, pad_factor=4):
    if pad_factor == 1:
        return array

    n, m = array.shape
    new_n = n * pad_factor
    new_m = m * pad_factor

    padded = np.zeros((new_n, new_m), dtype=array.dtype)

    start_n = (new_n - n) // 2
    start_m = (new_m - m) // 2

    padded[start_n:start_n + n, start_m:start_m + m] = array

    return padded


def compute_psf(pupil, pad_factor=4, normalize=True):
    padded_pupil = pad_array_centered(pupil, pad_factor=pad_factor)

    field = np.fft.fftshift(
        np.fft.fft2(
            np.fft.ifftshift(padded_pupil)
        )
    )

    psf = np.abs(field) ** 2

    if normalize:
        total = np.sum(psf)
        if total > 0:
            psf = psf / total

    return psf, field

现在定义物理参数：

wavelength = 550e-9      # 550 nm
focal_length = 50e-3     # 50 mm
pupil_diameter = 10e-3   # 10 mm

n_pupil = 256
pad_factor = 4
n_fft = n_pupil * pad_factor

x, y, rho, theta, mask = make_pupil_grid(n=n_pupil)

zero_opd = np.zeros_like(rho)

pupil, _ = pupil_function_from_opd(
    opd=zero_opd,
    mask=mask,
    wavelength=wavelength
)

psf, _ = compute_psf(
    pupil,
    pad_factor=pad_factor
)

归一化坐标网格从 $-1$ 到 $1$，代表整个光瞳直径。因此粗略的光瞳采样间隔为：

$$ \Delta u = \frac{D}{N_\text{pupil}} $$

代码中：

pupil_sample_spacing = pupil_diameter / n_pupil

dx_image = psf_sample_spacing_from_pupil(
    wavelength=wavelength,
    focal_length=focal_length,
    pupil_sample_spacing=pupil_sample_spacing,
    n_fft=n_fft
)

dx_image_mm = dx_image * 1e3

print("Image-plane sample spacing [mm]:", dx_image_mm)
print("Image-plane sample spacing [um]:", dx_image * 1e6)

现在我们可以为PSF坐标轴赋予物理意义。

def image_coordinate_axis(n, dx):
    """
    创建中心化的像平面坐标轴。
    """
    center = n // 2
    index = np.arange(n) - center
    return index * dx


coord_image_m = image_coordinate_axis(psf.shape[0], dx_image)
coord_image_um = coord_image_m * 1e6

extent_um = [
    coord_image_um[0],
    coord_image_um[-1],
    coord_image_um[0],
    coord_image_um[-1],
]

plt.figure(figsize=(5, 4))
plt.imshow(psf, origin="lower", extent=extent_um)
plt.colorbar(label="Normalized intensity")
plt.xlabel("Image x [um]")
plt.ylabel("Image y [um]")
plt.title("Diffraction-limited PSF with physical axis")
plt.tight_layout()
plt.show()

现在，这张PSF图不再仅仅是一个数组图像。其坐标轴具有物理含义。

这是一个巨大的进步。

但它也带来了责任——如果光瞳直径、焦距、波长或补零尺寸出错，坐标轴的标签就会随之出错。

6. 检验艾里斑半径

一个很好的合理性检查是艾里斑的第一暗环半径。

对于圆形孔径，焦平面上的第一暗环半径近似为：

$$ r_\text{Airy} = 1.22\lambda N $$

其中 $N$ 为F数。

代入：

λ = 550 nm
N = 5

得到：

$$ r_\text

1.22 \times 550\text{ nm} \times 5 $$

$$ r_\text{Airy} \approx 3.36\ \mu\text{m} $$

代码：

f_number = focal_length / pupil_diameter
airy_radius = 1.22 * wavelength * f_number

print("F-number:", f_number)
print("Airy first dark radius [um]:", airy_radius * 1e6)

这为我们提供了检验PSF坐标轴的一个粗略方法。

所绘PSF中的第一暗环应该接近该半径。

由于采样、有限阵列效应、显示插值以及我们定位最小值的方式，它或许不会完全匹配。但如果该环出现在3米或0.003纳米处，那就一定有问题。

在计算光学中，这种合理性检验不是可选的。

一个好习惯是：

在计算衍射受限PSF之后，
检查艾里斑的尺度在物理上是否合理。

仅此一项检验就能抓住许多单位上的错误。

7. 错误采样示例：光瞳直径错误下的PSF

让我们故意犯一个错误。

假设实际光瞳直径为10 mm，但我们错误地将其标记为1 mm。

计算出的PSF数组并不会改变。为什么？因为数组值来自归一化的光瞳掩膜。

但物理坐标轴却会剧烈变化。

wrong_pupil_diameter = 1e-3  # 错误：1 mm 而非 10 mm

wrong_pupil_sample_spacing = wrong_pupil_diameter / n_pupil

wrong_dx_image = psf_sample_spacing_from_pupil(
    wavelength=wavelength,
    focal_length=focal_length,
    pupil_sample_spacing=wrong_pupil_sample_spacing,
    n_fft=n_fft
)

print("Correct dx [um]:", dx_image * 1e6)
print("Wrong dx [um]:", wrong_dx_image * 1e6)

由于光瞳采样间隔变为了原来的十分之一，像平面采样间隔变为原来的十倍。

PSF图像看起来依然像艾里斑图样，但其物理尺寸已被错误标记。

这就是那种静默失效的模式。

数值阵列本身是不够的。坐标模型至关重要。

错误的图可能看起来很“光学”，甚至很漂亮，但它描述的是一个完全不同的物理系统。

8. 从PSF采样到MTF频率轴

现在我们向前一步。

MTF是由PSF计算得到的。其频率轴由PSF的采样间隔决定。

若PSF采样间隔为 $\Delta X$，则MTF的频率样本为：

$$ f_k = \frac{k - N/2}{N\Delta X} $$

NumPy中：

freq = np.fft.fftshift(np.fft.fftfreq(N, d=dx_image))

若 dx_image 以米为单位， freq 的单位即为 线对/米。

要转换为 线对/毫米：

freq_cyc_per_mm = freq / 1000.0

因为1米包含1000毫米，线对/米除以1000得到线对/毫米。

以下是一个包含坐标轴生成的整洁MTF计算：

def compute_otf_mtf(psf):
    otf = np.fft.fftshift(
        np.fft.fft2(
            np.fft.ifftshift(psf)
        )
    )

    center = (otf.shape[0] // 2, otf.shape[1] // 2)

    if np.abs(otf[center]) > 0:
        otf = otf / otf[center]

    mtf = np.abs(otf)

    return otf, mtf


def mtf_frequency_axis(n, dx_image):
    """
    生成MTF的频率轴。

    参数
    ----------
    n : int
        PSF / MTF 的采样点数。
    dx_image : float
        PSF 采样间隔，单位：米。

    返回
    -------
    freq_cyc_per_mm : 1维数组
        空间频率，单位：线对/毫米。
    """
    freq_cyc_per_m = np.fft.fftshift(
        np.fft.fftfreq(n, d=dx_image)
    )

    freq_cyc_per_mm = freq_cyc_per_m / 1000.0

    return freq_cyc_per_mm

现在进行计算并绘图：

otf, mtf = compute_otf_mtf(psf)

freq = mtf_frequency_axis(
    n=mtf.shape[0],
    dx_image=dx_image
)

center = mtf.shape[0] // 2
mtf_slice = mtf[center, :]

plt.figure(figsize=(6, 4))
plt.plot(freq[center:], mtf_slice[center:])
plt.xlabel("Spatial frequency [cycles/mm]")
plt.ylabel("MTF")
plt.ylim(0, 1.05)
plt.title("Diffraction-limited MTF")
plt.grid(True)
plt.tight_layout()
plt.show()

现在x轴具有了物理意义。

9. 检验衍射截止频率

对于非相干、衍射受限、圆形孔径的成像系统，截止频率近似为：

$$ f_c = \frac{1}{\lambda N} $$

其中：

• $f_c$ 的单位是 线对 / 长度；
• $\lambda$ 使用相同的长度单位；
• $N$ 为F数。

若 $\lambda$ 以毫米为单位，则 $f_c$ 的单位是 线对/毫米。

wavelength_mm = wavelength * 1e3
cutoff_cyc_per_mm = 1.0 / (wavelength_mm * f_number)

print("Diffraction cutoff [cycles/mm]:", cutoff_cyc_per_mm)

对于550 nm和F/5，截止频率近似为：

$$ f_c \approx 364\ \text{cycles/mm} $$

MTF应在该频率附近趋近于零。

再次强调，在小型教学脚本中结果可能并不完美，但应在合理区间内。

如果你的MTF图在3.64线对/毫米或3640000线对/毫米处趋零，请怀疑有单位错误。

截止频率提供了一个强力的合理性检查：

艾里半径检验PSF轴。
衍射截止频率检验MTF轴。

这两项检查应当成为自动化的步骤。

10. 错误的MTF轴示例

我们故意用错误的PSF采样间隔来生成错误的MTF频率轴。

MTF数值相同，但频率标签却变了。

wrong_freq = mtf_frequency_axis(
    n=mtf.shape[0],
    dx_image=wrong_dx_image
)

plt.figure(figsize=(7, 4))
plt.plot(freq[center:], mtf_slice[center:], label="Correct axis")
plt.plot(wrong_freq[center:], mtf_slice[center:], "--", label="Wrong axis")
plt.xlabel("Spatial frequency [cycles/mm]")
plt.ylabel("MTF")
plt.ylim(0, 1.05)
plt.title("Same MTF values, different frequency labels")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

这正是那种可以逃过肉眼检查的错误。

曲线形状完全相同，只有x轴发生了变化。

但对于工程解读来说，这绝非小问题。它会改变以下问题的答案：

在50线对/毫米处的MTF是多少？
截止频率是多少？
它与传感器像素间距相比如何？

一个频率轴被错误标记的系统可能将一个性能孱弱的系统误判为强大，或者相反。

仅靠曲线形状是不够的。

11. FFT平移（shift）：中心在哪？

FFT平移错误是另一个经典陷阱。

NumPy的FFT函数使用这样一种数组约定：在进行平移之前，零频率位于数组的起始位置。

为了显示和光学解读，我们通常希望将零频率或PSF峰值置于中心。

因此我们使用：

np.fft.fftshift(...)

有时还需：

np.fft.ifftshift(...)

由光瞳计算PSF的常见模式是：

field = np.fft.fftshift(
    np.fft.fft2(
        np.fft.ifftshift(pupil)
    )
)

由PSF计算OTF的常见模式是：

otf = np.fft.fftshift(
    np.fft.fft2(
        np.fft.ifftshift(psf)
    )
)

在变换过程中应保持中心化数组的中心化。

存在几种有效的约定，但混合使用不同约定是危险的。

一个实用的测试方法是：

peak_index = np.unravel_index(np.argmax(psf), psf.shape)
print("PSF peak index:", peak_index)
print("Expected center:", (psf.shape[0] // 2, psf.shape[1] // 2))

对于一个轴上理想光瞳，PSF的峰值应位于或非常接近中心。

如果它出现在角落，那意味着你用于显示的平移约定是错误的。

对于MTF：

otf_center = (otf.shape[0] // 2, otf.shape[1] // 2)
print("OTF center value:", otf[otf_center])
print("MTF center value:", mtf[otf_center])

归一化后，MTF的中心值应接近1。

这些检查很简单，请使用它们。

12. 奇数与偶数数组尺寸

偶数尺寸和奇数尺寸的数组在中心行为上略有不同。

对于一个 $N=1024$ 的偶数数组，fftshift 后的中心索引为：

N // 2

对于一个 $N=1025$ 的奇数数组，它也是：

N // 2

但中心周围的对称性略有不同。

这并不意味着奇数尺寸被禁止，而是意味着你应该保持一致性，并检查峰值和零频率出现的位置。

对于教学代码，2的幂次或偶数尺寸很方便：

256, 512, 1024

现代FFT库并不强制要求这些尺寸，但它们能使示例显得干净。

调试时，避免一次改变太多东西。不要同时更改：

• 数组尺寸；
• 补零因子；
• 孔径直径；
• 波长；
• FFT平移约定。

一次只改一件事，观察发生了什么。

这样做速度虽慢，但行之有效。

13. 补零：是插值，而非增加物理信息

补零很有用，但常被误解。

当你在计算PSF之前对光瞳进行补零时，你增加了傅里叶域输出的样本点数。

显示的PSF变得更平滑。

但物理孔径并未改变，光学分辨率也未提高。

正确的表述是：

补零是对采样后的PSF进行插值。
它不会创造新的光学信息。

让我们比较不同的补零因子。

pad_factors = [1, 2, 4, 8]

plt.figure(figsize=(7, 4))

for pf in pad_factors:
    psf_pf, _ = compute_psf(pupil, pad_factor=pf)

    n_fft_pf = n_pupil * pf
    dx_pf = psf_sample_spacing_from_pupil(
        wavelength=wavelength,
        focal_length=focal_length,
        pupil_sample_spacing=pupil_sample_spacing,
        n_fft=n_fft_pf
    )

    axis_pf = image_coordinate_axis(psf_pf.shape[0], dx_pf) * 1e6
    center_pf = psf_pf.shape[0] // 2
    slice_pf = psf_pf[center_pf, :]

    plt.plot(axis_pf, slice_pf / np.max(slice_pf), label=f"pad {pf}")

plt.xlim(-15, 15)
plt.xlabel("Image coordinate [um]")
plt.ylabel("Normalized peak intensity")
plt.title("Zero padding changes display sampling")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

这些曲线应描述同一个衍射图样，但更高的补零因子给出了更平滑的采样轮廓。

如果补零显著改变了物理图像，那意味着还有其他地方出了问题。

一个有用的检查是：

当补零改变时，艾里斑半径不应改变。
横跨艾里斑的采样点数量应该改变。

这句话是一个很好的调试工具。

14. PSF归一化：按峰值还是按能量？

PSF的显示通常有两种归一化方式。

能量归一化

PSF总和等于1：

$$ \sum \mathrm{PSF} = 1 $$

代码：

psf_energy = psf / np.sum(psf)

这在比较能量分布时很有用。

峰值归一化

PSF峰值等于1：

$$ \max(\mathrm{PSF}) = 1 $$

代码：

psf_peak = psf / np.max(psf)

这在绘制剖面图或视觉上比较形状时很有用。

它们回答的问题不同。

能量归一化保留了总能量的可比性；峰值归一化使中心峰值在视觉上具有可比性。

对于斯特列尔比，你需要匹配的能量归一化和采样：

strehl = np.max(aberrated_psf_energy) / np.max(perfect_psf_energy)

如果一个PSF是峰值归一化，另一个是能量归一化，这个比值就毫无意义。

这是那种可以静悄悄地潜伏在笔记本中的错误。

一个实用的习惯是：

清晰地命名归一化后的数组：
psf_energy_normalized
psf_peak_normalized

不要把一切都命名为 psf。

15. OTF归一化：除以零频分量

OTF是PSF的傅里叶变换。

OTF的零频值即为PSF的总能量。若PSF之和为1，平移后OTF的中心应为1。

但数值上的缩放和归一化选择仍可能有所变化。

因此我们通常这样做归一化：

center = (otf.shape[0] // 2, otf.shape[1] // 2)
otf = otf / otf[center]
mtf = np.abs(otf)

这使得：

零频率处的MTF = 1

一个基本检查：

print(mtf[center])

它应该非常接近：

1.0

如果不是，请检查：

• PSF中心化；
• OTF平移；
• OTF中心索引；
• PSF归一化；
• 偶发的复数或负PSF值。

零频频点是你的锚点，切勿忽视。

16. 错误归一化示例

让我们创建两个PSF：理想的与有像差的。

def synthetic_opd(rho, theta, mask, wavelength):
    opd = np.zeros_like(rho)

    defocus = 0.50 * wavelength
    astigmatism = 0.25 * wavelength

    opd += defocus * (2.0 * rho**2 - 1.0)
    opd += astigmatism * rho**2 * np.cos(2.0 * theta)

    opd[~mask] = 0.0

    return opd


aberrated_opd = synthetic_opd(
    rho=rho,
    theta=theta,
    mask=mask,
    wavelength=wavelength
)

aberrated_pupil, _ = pupil_function_from_opd(
    opd=aberrated_opd,
    mask=mask,
    wavelength=wavelength
)

perfect_psf, _ = compute_psf(
    pupil,
    pad_factor=4,
    normalize=True
)

aberrated_psf, _ = compute_psf(
    aberrated_pupil,
    pad_factor=4,
    normalize=True
)

现在正确地计算斯特列尔比：

strehl_correct = np.max(aberrated_psf) / np.max(perfect_psf)
print("Correct Strehl:", strehl_correct)

现在，我们通过先将两者峰值归一化来进行一次错误比较：

perfect_peak_norm = perfect_psf / np.max(perfect_psf)
aberrated_peak_norm = aberrated_psf / np.max(aberrated_psf)

strehl_wrong = np.max(aberrated_peak_norm) / np.max(perfect_peak_norm)
print("Wrong Strehl after peak normalization:", strehl_wrong)

错误的结果将会是1。

这显然毫无用处——在计算峰值比之前，我们亲手抹去了峰值的差异。

这个例子虽然简单，但其教训具有普适性：

在决定了你要比较什么量之后再进行归一化。

而非在此之前。

17. 孔径掩膜：方形数组陷阱

一个方形的NumPy数组并非圆形孔径。

这种错误很容易犯：

pupil_wrong = np.ones((256, 256), dtype=complex)

这并不代表一个清晰的圆形透镜孔径，它代表的是一个方形孔径。

方形孔径的衍射图样和MTF均与圆形不同。

正确的圆形光瞳需要一个掩膜：

pupil_correct = mask.astype(complex)

让我们比较一下：

square_pupil = np.ones_like(pupil, dtype=complex)

square_psf, _ = compute_psf(square_pupil, pad_factor=4)
circular_psf, _ = compute_psf(pupil, pad_factor=4)

square_otf, square_mtf = compute_otf_mtf(square_psf)
circular_otf, circular_mtf = compute_otf_mtf(circular_psf)

center = circular_mtf.shape[0] // 2
freq = mtf_frequency_axis(circular_mtf.shape[0], dx_image)

plt.figure(figsize=(7, 4))
plt.plot(freq[center:], circular_mtf[center, center:], label="Circular pupil")
plt.plot(freq[center:], square_mtf[center, center:], "--", label="Square pupil")
plt.xlabel("Spatial frequency [cycles/mm]")
plt.ylabel("MTF")
plt.ylim(0, 1.05)
plt.title("Aperture shape changes MTF")
plt.grid(True)
plt.legend()
plt.tight_layout()
plt.show()

这种差异并非数值上的偶然，而是物理上的必然。

孔径形状会影响衍射。

因此，在计算PSF之前，第一个检查就是：

这个光瞳数组代表的是我所认为的那个孔径吗？

在对其进行傅里叶变换之前，先绘制光瞳的振幅图。

每次都请这样做。

18. 端点采样：linspace 可能悄然改变你的网格

在之前的章节中，我们经常使用：

coord = np.linspace(-1.0, 1.0, n)

默认情况下，这会包含两个端点。

对于许多教学绘图来说，这是可以接受的。但对于FFT采样，同时包含 $-1$ 和 $1$ 会略显尴尬，因为这些端点代表了周期性采样区间的重复边界。

一个对FFT更友好的版本是：

coord = np.linspace(-1.0, 1.0, n, endpoint=False)

这会在不重复端点的情况下对区间进行采样。

使用默认设置会毁掉你的结果吗？对于初学者示例，通常不会。但它可能引起微小的非对称性和采样差异。

对于细致的FFT工作，推荐使用显式采样：

coord = (np.arange(n) - n / 2) / (n / 2)

或：

coord = np.linspace(-1.0, 1.0, n, endpoint=False)

更深层的规则是：

了解你的网格是否包含端点。
不要让默认的采样选择悄无声息地定义你的光学模型。

这也是为什么生产级光学软件拥有精细的内部采样约定的原因之一。

19. 单位：最枯燥也最强大的调试工具

单位并非文书装饰。

它们是计算的一部分。

以下是本书此部分最常见的单位配对：

OPD: 米 (meters)
wavelength: 米
phase: 弧度

pupil diameter: 米
focal length: 米
PSF coordinate: 米

PSF sample spacing: 毫米
MTF frequency: 线对/毫米

危险的情况是混合了隐藏单位：

wavelength = 550       # 本意是纳米
focal_length = 50      # 本意是毫米
pupil_diameter = 10    # 本意是毫米

这段代码可能可以运行，但它完全没有可靠的物理意义。

更安全的编写方式是：

wavelength = 550e-9
focal_length = 50e-3
pupil_diameter = 10e-3

然后仅在显示时进行转换：

print(dx_image * 1e6, "um")
print(freq / 1000.0, "cycles/mm")

一条简单的纪律很有帮助：

在计算内部使用国际单位制 (SI)。
仅在输入和输出的边界处转换单位。

这个习惯能预防大量错误。

20. 符号约定与镜像结果

在第11章中我们曾提及，光瞳相位的符号约定可以不同：

$$ P(x,y)=A(x,y)\exp\left(i\frac{2\pi W(x,y)}{\lambda}\right) $$

或者：

$$ P(x,y)=A(x,y)\exp\left(-i\frac{2\pi W(x,y)}{\lambda}\right) $$

改变符号可能会根据像差和傅里叶约定的不同，镜像或共轭部分场行为。

对于对称像差，强度PSF看起来可能相似。对于非对称像差（如彗差类相位），PSF的方向可能会翻转。

这意味着两种实现可能在视觉上不一致，但它们各自内部可能是一致的。

一个实用的对比检查清单：

相同的 OPD 定义
相同的相位符号
相同的光瞳坐标取向
相同的 FFT 约定
相同的显示取向

不要只对比最终的图像，要对中间数组进行比对：

• OPD图；
• 相位图；
• 光瞳振幅；
• 光瞳相位；
• PSF中心；
• MTF中心值。

当最终的PSF出现镜像时，错误更可能存在于符号约定或坐标取向上，而非衍射理论中。

21. 子午与弧矢方向：不要猜测

对于具有旋转对称性的轴上系统，水平和垂直方向的MTF剖面可能几乎相同。

对于轴外视场，子午和弧矢方向至关重要。

但它们并非在所有图像中都是简单的“x”和“y”。

它们的定义与以下因素相关：

• 光轴；
• 视场点；
• 主光线或子午面。

因此，如下做法是草率的：

horizontal = tangential
vertical = sagittal

在特定的坐标设置下，这或许成立，但并非普遍成立。

更安全的表述是：

水平和垂直剖面是数组方向。
子午和弧矢剖面是光学视场方向。

要与Optiland或任何光学设计软件进行对比，应检查该软件如何为所选视场定义这些方向。

一个好的调试问题是：

我数组中哪个方向对应于视场的子午面？

没有这个答案，诸如“子午”和“弧矢”的标签就只是装饰，不可信赖。

22. 一个小型诊断工具箱

现在让我们收集几个能让调试更轻松的函数。

它们并不花哨，但很实用。

def check_array_center(name, array):
    """
    打印实值数组的峰值位置。
    对 PSF 和 MTF 检查非常有用。
    """
    peak = np.unravel_index(np.argmax(np.abs(array)), array.shape)
    center = (array.shape[0] // 2, array.shape[1] // 2)

    print(f"{name}")
    print("  shape: ", array.shape)
    print("  peak:  ", peak)
    print("  center:", center)
    print()


def check_psf_normalization(psf):
    """
    打印基本的 PSF 归一化信息。
    """
    print("PSF checks")
    print("  sum: ", np.sum(psf))
    print("  max: ", np.max(psf))
    print("  min: ", np.min(psf))
    print()


def check_mtf_normalization(mtf):
    """
    打印 MTF 的中心值。
    """
    center = (mtf.shape[0] // 2, mtf.shape[1] // 2)

    print("MTF checks")
    print("  center value:", mtf[center])
    print("  max value:   ", np.max(mtf))
    print("  min value:   ", np.min(mtf))
    print()

使用它们：

check_array_center("PSF", psf)
check_psf_normalization(psf)

otf, mtf = compute_otf_mtf(psf)
check_array_center("MTF", mtf)
check_mtf_normalization(mtf)

对于一个中心化的衍射受限PSF：

• PSF峰值应靠近中心；
• 若为能量归一化，PSF之和应为1；
• MTF中心值应为1；
• MTF最大值通常应接近1。

这些检查并不能证明一切正确，但它们能迅速捕捉到常见的失败。

23. 一个完整的修正后PSF/MTF工作流

这是一个将簿记工作保持显式的紧凑工作流。

import numpy as np
import matplotlib.pyplot as plt


def make_pupil_grid(n=256):
    coord = np.linspace(-1.0, 1.0, n, endpoint=False)
    x, y = np.meshgrid(coord, coord)

    rho = np.sqrt(x**2 + y**2)
    theta = np.arctan2(y, x)
    mask = rho <= 1.0

    return x, y, rho, theta, mask


def pupil_function_from_opd(opd, mask, wavelength, amplitude=None):
    if amplitude is None:
        amplitude = mask.astype(float)

    phase = 2.0 * np.pi * opd / wavelength
    pupil = amplitude * np.exp(1j * phase)
    pupil[~mask] = 0.0 + 0.0j

    return pupil, phase


def pad_array_centered(array, pad_factor=4):
    if pad_factor == 1:
        return array

    n, m = array.shape
    new_n = n * pad_factor
    new_m = m * pad_factor

    padded = np.zeros((new_n, new_m), dtype=array.dtype)

    start_n = (new_n - n) // 2
    start_m = (new_m - m) // 2

    padded[start_n:start_n + n, start_m:start_m + m] = array

    return padded


def compute_psf(pupil, pad_factor=4):
    padded_pupil = pad_array_centered(pupil, pad_factor=pad_factor)

    field = np.fft.fftshift(
        np.fft.fft2(
            np.fft.ifftshift(padded_pupil)
        )
    )

    psf = np.abs(field) ** 2
    psf = psf / np.sum(psf)

    return psf, field


def compute_otf_mtf(psf):
    otf = np.fft.fftshift(
        np.fft.fft2(
            np.fft.ifftshift(psf)
        )
    )

    center = (otf.shape[0] // 2, otf.shape[1] // 2)
    otf = otf / otf[center]

    mtf = np.abs(otf)

    return otf, mtf


def psf_sample_spacing_from_pupil(
    wavelength,
    focal_length,
    pupil_sample_spacing,
    n_fft
):
    dtheta = wavelength / (n_fft * pupil_sample_spacing)
    return focal_length * dtheta


def image_coordinate_axis(n, dx):
    center = n // 2
    return (np.arange(n) - center) * dx


def mtf_frequency_axis(n, dx_image):
    freq_cyc_per_m = np.fft.fftshift(
        np.fft.fftfreq(n, d=dx_image)
    )
    return freq_cyc_per_m / 1000.0


# 物理参数，计算内部使用国际单位制 (SI)
wavelength = 550e-9
focal_length = 50e-3
pupil_diameter = 10e-3

# 采样参数
n_pupil = 256
pad_factor = 4
n_fft = n_pupil * pad_factor

# 构建光瞳
x, y, rho, theta, mask = make_pupil_grid(n=n_pupil)
opd = np.zeros_like(rho)

pupil, phase = pupil_function_from_opd(
    opd=opd,
    mask=mask,
    wavelength=wavelength
)

# 计算 PSF
psf, field = compute_psf(
    pupil=pupil,
    pad_factor=pad_factor
)

# 坐标簿记
pupil_sample_spacing = pupil_diameter / n_pupil

dx_image = psf_sample_spacing_from_pupil(
    wavelength=wavelength,
    focal_length=focal_length,
    pupil_sample_spacing=pupil_sample_spacing,
    n_fft=n_fft
)

image_axis_um = image_coordinate_axis(psf.shape[0], dx_image) * 1e6

# 计算 MTF
otf, mtf = compute_otf_mtf(psf)
freq_cyc_per_mm = mtf_frequency_axis(mtf.shape[0], dx_image)

# 合理性检查
f_number = focal_length / pupil_diameter
airy_radius_um = 1.22 * wavelength * f_number * 1e6
cutoff_cyc_per_mm = 1.0 / ((wavelength * 1e3) * f_number)

print("F-number:", f_number)
print("Airy first dark radius [um]:", airy_radius_um)
print("PSF sample spacing [um]:", dx_image * 1e6)
print("Diffraction cutoff [cycles/mm]:", cutoff_cyc_per_mm)

# 绘制 PSF
extent_um = [
    image_axis_um[0],
    image_axis_um[-1],
    image_axis_um[0],
    image_axis_um[-1],
]

plt.figure(figsize=(5, 4))
plt.imshow(psf, origin="lower", extent=extent_um)
plt.colorbar(label="Normalized intensity")
plt.xlabel("Image x [um]")
plt.ylabel("Image y [um]")
plt.title("PSF with physical coordinates")
plt.tight_layout()
plt.show()

# 绘制 MTF 剖面
center = mtf.shape[0] // 2

plt.figure(figsize=(6, 4))
plt.plot(freq_cyc_per_mm[center:], mtf[center, center:])
plt.xlabel("Spatial frequency [cycles/mm]")
plt.ylabel("MTF")
plt.ylim(0, 1.05)
plt.title("MTF with physical frequency axis")
plt.grid(True)
plt.tight_layout()
plt.show()

这并非可能的最短代码。

这是刻意为之。

此版本更青睐于可追溯性而非紧凑性。

每一个物理参数都清晰可见。每一个采样步骤都明确写出。每一个输出坐标轴都源自已知量。

这种代码更容易让人信赖。

24. Optiland 对比：在指责库之前要检查什么

将本教学代码与Optiland进行对比时，不要期待能偶然匹配。

一个真正的光学设计库可能在以下方面包含更精细的处理：

• 出瞳几何；
• 参考球面；
• 视场相关的光瞳映射；
• 波长权重；
• 像空间折射率；
• 焦点位置；
• 渐晕；
• 偏振；
• 采样策略；
• 内部归一化；
• 子午/弧矢定义。

因此，不匹配并不自动意味着是错误。

在指责库或你的代码之前，请检查：

1. 同一光学系统？
2. 同一视场点？
3. 同一波长？
4. 同一孔径光阑和光瞳直径？
5. 同一焦平面？
6. 同一 OPD 参考？
7. 同一光瞳采样？
8. 同一 FFT 补零？
9. 同一 PSF 归一化？
10. 同一 MTF 频率单位？
11. 同一子午/弧矢定义？
12. 同一显示比例？

一个用于验证的对比工作流大致如下：

# 版本检查的工作流程概要。
# 检查你安装的 Optiland 版本以获取准确的 API 名称。

# 1. 构建或加载光学系统。
# optic = ...

# 2. 设置波长、视场点、孔径和焦点。
# wavelength = ...
# field = ...

# 3. 使用 Optiland 的分析工具生成 PSF 和 MTF。
# psf_optiland = ...
# mtf_optiland = ...

# 4. 如果可能，提取或重建可比较的数组。
# psf_array = ...
# dx_image = ...

# 5. 从采样后的 PSF 重新计算 OTF/MTF。
# otf_numpy, mtf_numpy = compute_otf_mtf(psf_array)

# 6. 比较：
# - PSF 峰值位置
# - PSF 总和
# - 零频率处的 MTF
# - 频率轴
# - 截止频率
# - 子午/弧矢方向定义

这里教育性的目的，并非要证明你简短的脚本比Optiland更正确——那是不可能的。

其目的是理解，在两种计算能够相互比较之前，哪些方面必须对齐。

这是一种非常具有实践意义的理解。

25. 可疑PSF的调试序列

当一个PSF看起来不对劲时，不要随意改动代码。

请遵循一个序列。

步骤1：检查光瞳振幅

绘制：

plt.imshow(np.abs(pupil), origin="lower")
plt.colorbar(label="Amplitude")
plt.title("Pupil amplitude")
plt.show()

追问：

孔径形状正确吗？
中心遮拦是按预期包含还是排除了？
孔径之外是否有非零值？

步骤2：检查光瞳相位

绘制：

plt.imshow(np.angle(pupil), origin="lower")
plt.colorbar(label="Phase [rad]")
plt.title("Wrapped pupil phase")
plt.show()

追问：

相位图案是否与预期的 OPD 匹配？
相位包裹是否在预期之内？
单位是否一致？

步骤3：检查PSF中心化

check_array_center("PSF", psf)

追问：

对于轴上点，峰值是否在中心附近？

步骤4：检查PSF能量

check_psf_normalization(psf)

追问：

归一化后，PSF 总和是否等于 1？

步骤5：检查物理尺度

计算：

airy_radius = 1.22 * wavelength * f_number

追问：

艾里斑的尺度是否合理？

这个序列能捕捉到大多数初学者的错误。

26. 可疑MTF的调试序列

当MTF曲线看起来不对劲时，请使用另一个序列。

步骤1：首先检查PSF

不要还没检查PSF就去调试MTF。

一个错误的PSF不可能产生正确的MTF。

步骤2：检查OTF中心

center = (otf.shape[0] // 2, otf.shape[1] // 2)
print(otf[center])

归一化后，中心应该接近：

1 + 0j

步骤3：检查MTF中心

print(mtf[center])

它应该接近：

1

步骤4：检查频率轴

计算衍射截止频率：

cutoff_cyc_per_mm = 1.0 / ((wavelength * 1e3) * f_number)

追问：

MTF 是否在预期截止频率附近趋于零？

步骤5：检查方向标签

追问：

这些是数组切片，还是真正的子午/弧矢切片？

步骤6：检查显示范围

确保图像没有隐藏负值、被裁剪或超出范围的值。

对于 MTF：

plt.ylim(0, 1.05)

通常对初步检查很有用。

27. 最常见的失败模式

以下是一张关于症状及其可能原因的紧凑表格。

这张表并不能代替理解，它是一个分流工具。

当出现问题时，它告诉了你首先应该去检查哪些地方。

28. 为什么本章要放在优化前面

人们或许会忍不住从MTF直接跳到优化。

毕竟，一旦我们能计算MTF，为什么不立刻开始改进它呢？

因为优化会放大错误。

如果PSF的尺度错误，优化器可能去改进一个错误的指标。

如果MTF频率轴错误，优化器可能去青睐错误的空间频率。

如果归一化不一致，优化器可能通过改变能量尺度而非图像质量来表面上提高性能。

如果子午和弧矢方向被标错，优化器可能修正的是错误的方向行为。

优化并不会让一个糟糕的指标变得诚实。

它只会忠实地遵循你给它的那个指标。

这就是本章坐落于此的原因。

在要求软件去改善一个镜头之前，我们需要清楚我们的数字意味着什么。

现在，这条链路应该带着安全措施来解读：

光瞳函数
→ 经过检查的振幅与相位
→ 具有已知采样和归一化的 PSF
→ 在零频处归一化的 OTF
→ 具有正确频率轴和方向标签的 MTF

只有如此，才适宜在评价函数中使用这些量。

29. 可实践的检查清单

每当你计算PSF或MTF时，都使用这份检查清单。

光瞳检查

[ ] 孔径掩膜正确吗？
[ ] 振幅场是振幅而非强度吗？
[ ] OPD 是否与波长使用相同的长度单位？
[ ] 相位符号约定是否已知？
[ ] 光瞳坐标取向是否已知？
[ ] 光瞳是否采样得足够密集？

PSF检查

[ ] FFT平移约定是否一致？
[ ] 对于轴上点，PSF峰值是否居中？
[ ] 需要时是否进行了PSF能量归一化？
[ ] 补零是否被理解为插值而非更好的光学？
[ ] PSF坐标轴是否由光瞳采样、波长和焦距推导得出？
[ ] 对于一个完美的圆形光瞳，艾里斑半径看起来合理吗？

MTF检查

[ ] OTF 是否从 PSF 计算，而非直接从 OPD 或光瞳计算？
[ ] OTF 是否由其零频值归一化？
[ ] MTF 是否从 1 开始？
[ ] 频率轴是否使用了预期的单位？
[ ] 衍射截止频率看起来合理吗？
[ ] 子午和弧矢方向被正确定义了吗？

对比检查

[ ] 相同波长？
[ ] 相同视场点？
[ ] 相同焦平面？
[ ] 相同孔径和光瞳定义？
[ ] 相同采样密度？
[ ] 相同归一化？
[ ] 相同显示比例？
[ ] 相同坐标取向？

这份检查清单看起来或许平淡无奇。但在实际工作中，平淡的清单能节省数小时的时间。

30. 本章为完整链路增添了什么

前面的章节给了我们这些公式：

$$ P(x,y)=A(x,y)\exp\left(i\frac{2\pi W(x,y)}{\lambda}\right) $$

$$ \mathrm{PSF}=|\mathcal{F}{P(x,y)}|^2 $$

$$ \mathrm{OTF}=\mathcal{F}{\mathrm{PSF}} $$

$$ \mathrm{MTF}=|\mathrm{OTF}| $$

本章增添了使这些公式在代码中有意义的条件。

它告诉我们，每一次变换都需要采样，每一张图都需要单位，每一次比较都需要归一化。

因此，真正的计算链路不仅仅是：

OPD
→ pupil
→ PSF
→ MTF

它更接近于：

具有正确单位的 OPD
→ 具有正确掩膜和相位约定的 pupil
→ 具有正确 FFT 平移、采样和归一化的 PSF
→ 在零频处归一化的 OTF
→ 具有正确频率轴和方向标签的 MTF

这个更长的版本不那么优雅。

但它也更接近真相。

下一章将进入优化。我们终于要问，光学软件在“改进”一个设计时究竟在做什么。

但现在，我们可以安全地提出这个问题。

优化并非魔法，它是由评价函数引导的重复计算。

经过本章，我们明白了一件重要的事情：

在优化一个光学指标之前，先确保这个指标的含义你已了然于胸。

本章小结

基于FFT的PSF和MTF计算功能强大，但对采样、单位、中心化和归一化非常敏感。

DFT的采样关系为：

$$ \Delta f = \frac{1}{N\Delta x} $$

由物理光瞳计算PSF时，近似的焦平面采样间隔为：

$$ \Delta X = \frac{F\lambda}{N\Delta u} $$

其中 $F$ 是焦距，$\lambda$ 是波长，$N$ 是补零后的FFT尺寸，$\Delta u$ 是光瞳平面的采样间隔。

对于MTF，频率轴由PSF采样间隔得到：

freq = np.fft.fftshift(np.fft.fftfreq(n, d=dx_image))

PSF应通过艾里斑半径进行检验：

$$ r_\text{Airy}=1.22\lambda N $$

MTF应通过衍射截止频率进行检验：

$$ f_c=\frac{1}{\lambda N} $$

其中 $N$ 是F数。

补零能改善显示的采样，但并不会提高物理光学分辨率。PSF能量归一化、PSF峰值归一化以及OTF零频归一化服务于不同的目的，不应随意混用。

最重要的实践习惯是检查中间结果：

光瞳振幅
光瞳相位
PSF中心与能量
OTF中心值
MTF频率轴

只有当采样和单位都值得信赖时，图像才值得信赖。

生成的验证图像