第 12 章：如何从光瞳函数计算 PSF

点源应该是最简单的物体。

它没有大小。它没有纹理。它没有边缘，没有图案，没有细节。在几何光学中，一个完美透镜会把来自该点的所有光线汇聚到像平面中的一个点。

但真实的光学系统不会形成数学上的点。

即使是完美的圆形孔径也会把一个点扩散成一个小衍射斑。像差会使它进一步扩散、变形、将能量移出中心，或产生不对称的裙边。

那个扩散后的点的像就是点扩散函数，或称PSF。

前一章构建了衍射计算所需的物体：

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

这是复光瞳函数。它存储孔径振幅和波前相位。

现在我们让那个光瞳函数形成一个像。

简化的计算链是：

复光瞳函数
→ 傅里叶变换
→ 复像平面场
→ 取幅度平方
→ PSF

用代码表示，骨架很短：

field = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(pupil)))
psf = np.abs(field) ** 2
psf = psf / np.sum(psf)

这三行就是本章的核心。

但它们不是魔法。我们需要拆解它们的含义，它们使用的假设，以及可能出错的地方。

1. PSF 代表什么

PSF 是点状物体在像平面产生的强度分布。

如果光学系统是完美的且没有衍射，PSF 将是一个无限小的点。但物理光有波长，孔径有有限尺寸。有限孔径不能将光聚焦成一个无限小的点。

对于一个圆形、无像差的光瞳，理想的衍射斑就是我们熟悉的艾里斑：

• 一个明亮的中央核；
• 环绕的圆环；
• 强度从中心向外递减。

这就是衍射受限的 PSF。

当存在像差时，光瞳相位不再平坦。该畸变复光瞳的傅里叶变换在像平面重新分配能量。PSF 的形状改变。

所以 PSF 不仅仅是一张漂亮的图片。它回答了一个具体的成像问题：

如果物体包含一个理想的点，系统在像中画出怎样的强度图样？

一旦知道了这一点，我们就非常接近理解模糊、分辨率和 MTF。

2. 从光瞳平面到像平面

光瞳函数存在于光瞳平面。它告诉我们孔径上各点的复场。

PSF 存在于像平面。它告诉我们聚焦点像周围的强度。

在通常的标量夫琅禾费衍射近似下，焦区复场正比于光瞳函数的傅里叶变换：

$$U(\xi, \eta) \propto \mathcal{F}{P(x, y)}$$

PSF 是那个复场的幅度平方：

$$\mathrm{PSF}(\xi, \eta) = |U(\xi, \eta)|^2$$

将两者结合：

$$\mathrm{PSF}(\xi, \eta) \propto \left|\mathcal{F}{P(x, y)}\right|^2$$

其中：

• $P(x,y)$ 是复光瞳函数；
• $U(\xi,\eta)$ 是复像平面场；
• $\xi,\eta$ 是像平面或角度坐标；
• $\mathcal{F}$ 表示傅里叶变换。

比例符号很重要。真实的物理 PSF 需要正确的坐标缩放和能量归一化。为了学习计算链，我们可以先计算一个归一化的 PSF：

$$\sum_{\xi,\eta}\mathrm{PSF}(\xi,\eta)=1$$

这意味着把 PSF 视为一个能量分布。总能量为 1，我们研究这些能量去了哪里。

这种归一化简捷而有用：

psf = psf / np.sum(psf)

它不能解决所有的物理缩放问题，但能让不同光瞳函数之间的对比更加清晰。

3. 为什么傅里叶变换会出现

傅里叶变换之所以出现，是因为像平面场是通过光瞳上所有点的相干叠加构建的。

光瞳上的每一点贡献一个小波。在远场或焦平面近似下，这些贡献之间的相位关系随像平面位置变化。傅里叶变换就是为每个输出位置将所有这些复贡献相加的数学操作。

这是光线图所隐藏的部分。

光线图说：

光瞳的这部分将光射向这个方向。

衍射计算说：

所有光瞳点都贡献复振幅，这些振幅发生干涉。

PSF 是由干涉产生的。

这就是为什么光瞳函数必须是复的。如果我们仅保留 OPD 作为一个实值映射，我们会知道波前误差，但还没有相干叠加所需的对象。

现在的计算链看起来是这样的：

OPD 映射
→ 相位映射
→ 复光瞳函数
→ 相干傅里叶叠加
→ 像平面强度

这就是波前误差变成可见模糊的时刻。

4. 一个完美圆形光瞳

从最简单的情况开始：一个清晰、无像差的圆形孔径。

在光瞳内部：

$$A(x,y)=1$$

$$W(x,y)=0$$

所以：

$$P(x,y)=1$$

在光瞳外部：

$$P(x,y)=0$$

这个平坦的圆形光瞳产生衍射受限的艾里斑。

让我们复用前一章的光瞳网格概念。

import numpy as np
import matplotlib.pyplot as plt


def make_pupil_grid(n=256):
    coord = np.linspace(-1.0, 1.0, n)
    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

现在创建一个完美光瞳：

wavelength = 550e-9  # meters

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

zero_opd = np.zeros_like(rho)

perfect_pupil, perfect_phase = pupil_function_from_opd(
    opd=zero_opd,
    mask=mask,
    wavelength=wavelength
)

光瞳是复的，但在这种完美情况下，其相位在孔径内为零。所以复值简单地为：

$$1 + 0i$$

在圆形内部。

现在我们准备好计算 PSF。

5. 用 FFT 计算 PSF

离散傅里叶变换通过 numpy.fft.fft2 计算。

因为我们通常希望中央峰在显示图像的中心，我们使用 fftshift。

一个谨慎的版本在变换前使用 ifftshift，变换后使用 fftshift：

def compute_psf(pupil, normalize=True):
    """
    Compute a normalized PSF from a complex pupil function.

    Parameters
    ----------
    pupil : 2D complex array
        Complex pupil function.
    normalize : bool
        If True, normalize total PSF energy to 1.

    Returns
    -------
    psf : 2D array
        Intensity PSF.
    field : 2D complex array
        Complex field in the image plane.
    """
    field = np.fft.fftshift(
        np.fft.fft2(
            np.fft.ifftshift(pupil)
        )
    )

    psf = np.abs(field) ** 2

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

    return psf, field

运行它：

perfect_psf, perfect_field = compute_psf(perfect_pupil)

这个 perfect_psf 数组就是我们采样的圆形光瞳的衍射受限 PSF。

它还没有以微米或角秒为单位标记。这是一个采样的傅里叶平面结果。我们稍后在本章回到坐标缩放。

首先，看一下图样。

6. 显示 PSF

PSF 通常有一个非常亮的中央峰和弱得多的周围结构。如果你用线性度显示它，中央峰会强烈地主导，以致于圆环很难看到。

所以我们经常同时显示：

1. 线性强度；
2. 对数强度。

def show_psf(psf, title="PSF"):
    plt.figure(figsize=(5, 4))
    plt.imshow(psf, origin="lower")
    plt.colorbar(label="Normalized intensity")
    plt.title(title + " - linear scale")
    plt.xlabel("Image sample x")
    plt.ylabel("Image sample y")
    plt.tight_layout()
    plt.show()

    plt.figure(figsize=(5, 4))
    plt.imshow(np.log10(psf + 1e-12), origin="lower")
    plt.colorbar(label="log10 intensity")
    plt.title(title + " - log scale")
    plt.xlabel("Image sample x")
    plt.ylabel("Image sample y")
    plt.tight_layout()
    plt.show()


show_psf(perfect_psf, title="Diffraction-limited circular pupil")

对数图并不是另一个 PSF。它是同一个 PSF 用了一种显示变换来展示。

这个区别很重要。对数尺度的 PSF 会让弱的光环看起来视觉上很重要。线性尺度的 PSF 会隐藏仍然影响 MTF 的结构。好的分析往往需要这两种视角。

一个合理的习惯是：

用线性尺度来判断能量集中度。
用对数尺度来检视微弱结构。

中央峰通常对清晰度最重要，但周围微弱的能量告诉你对比度去了哪里。

7. 零填充：让显示的 PSF 更平滑

FFT 输出的采样取决于输入数组的大小。如果光瞳数组很小，PSF 可能看起来像块状或欠采样。

一种常见技术是在 FFT 之前对光瞳进行零填充。

零填充并不添加新的光学信息。它更精细地插值傅里叶域的显示。

这是一个带填充的 PSF 函数：

def pad_array_centered(array, pad_factor=2):
    """
    Zero-pad a 2D array around its center.

    Parameters
    ----------
    array : 2D array
        Input array.
    pad_factor : int
        Output size will be pad_factor times the input size.

    Returns
    -------
    padded : 2D array
        Centered zero-padded array.
    """
    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=2, normalize=True):
    """
    Compute a PSF from a complex pupil function, with optional zero padding.
    """
    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

现在计算一个更平滑的显示：

perfect_psf, perfect_field = compute_psf(
    perfect_pupil,
    pad_factor=4
)

show_psf(perfect_psf, title="Diffraction-limited PSF with zero padding")

再次强调，零填充不是神奇的“分辨率”。它给了我们在相同物理特征之间更多的输出采样。

一个有用的说法是：

零填充让绘制的 PSF 更平滑。
它并不让光学系统更清晰。

这句话能避免很多错误解读。

8. 解读完美 PSF

完美圆形光瞳的 PSF 有三个主要特征。

第一，中央峰紧凑且对称。这就是艾里核。

第二，圆环围绕中央峰。这些圆环不是像差。它们是由有限圆形孔径的衍射产生的。

第三，由于光瞳是径向对称且相位平坦，整个图样是径向对称的。

这是一个重要的基准。

如果光学系统有一个圆形清晰孔径且无像差，PSF 并不是一个点。它已经被衍射扩散了。

所以当有人说“衍射受限”时，他们并不是指“无限清晰”。他们指的是系统主要受限于衍射而非像差。

在实践中：

完美波前 + 有限孔径
→ 衍射受限 PSF

畸变波前 + 有限孔径
→ 像差 PSF

PSF 只有相对于这个基准才能被正确判断。正确的比较对象不是数学上的点，而是在相同孔径和波长下的衍射受限 PSF。

9. 添加像差

现在让我们复用一个合成的 OPD 映射。我们将制作一个类似离焦和像散的波前，就像前一章那样。

def synthetic_opd(rho, theta, mask, wavelength):
    """
    Create a synthetic OPD map in meters.

    This is a teaching wavefront, not a full optical design result.
    """
    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, aberrated_phase = pupil_function_from_opd(
    opd=aberrated_opd,
    mask=mask,
    wavelength=wavelength
)

aberrated_psf, aberrated_field = compute_psf(
    aberrated_pupil,
    pad_factor=4
)

现在显示它：

show_psf(aberrated_psf, title="Aberrated PSF")

与完美 PSF 相比，像差 PSF 应该显示出能量从中心重新分配。根据像差的不同，中央峰可能变宽、变得不对称，或发展出侧翼结构。

这是本章核心的视觉课程：

像差是光瞳中的相位误差。
PSF 显示了该相位误差把像的能量送到了哪里。

波前图可能看起来很抽象。PSF 则展示了其后果。

10. 对比完美和像差 PSF

并列显示很有用。

def show_two_psfs(psf_a, psf_b, title_a, title_b, log_scale=False):
    if log_scale:
        image_a = np.log10(psf_a + 1e-12)
        image_b = np.log10(psf_b + 1e-12)
        label = "log10 intensity"
    else:
        image_a = psf_a
        image_b = psf_b
        label = "Normalized intensity"

    vmin = min(np.nanmin(image_a), np.nanmin(image_b))
    vmax = max(np.nanmax(image_a), np.nanmax(image_b))

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

    plt.subplot(1, 2, 1)
    plt.imshow(image_a, origin="lower", vmin=vmin, vmax=vmax)
    plt.colorbar(label=label)
    plt.title(title_a)
    plt.xlabel("Image sample x")
    plt.ylabel("Image sample y")

    plt.subplot(1, 2, 2)
    plt.imshow(image_b, origin="lower", vmin=vmin, vmax=vmax)
    plt.colorbar(label=label)
    plt.title(title_b)
    plt.xlabel("Image sample x")
    plt.ylabel("Image sample y")

    plt.tight_layout()
    plt.show()


show_two_psfs(
    perfect_psf,
    aberrated_psf,
    "Perfect PSF",
    "Aberrated PSF",
    log_scale=False
)

show_two_psfs(
    perfect_psf,
    aberrated_psf,
    "Perfect PSF",
    "Aberrated PSF",
    log_scale=True
)

线性对比告诉我们有多少能量仍集中在中心附近。

对数对比显示微弱的结构。

两者都很重要。如果中央能量下降，精细细节的对比度通常会受损。如果微弱能量扩散到远处，图像会获得光晕或低水平模糊。

这就是为什么 PSF 比单个数字更有信息。它显示了模糊的空间形状。

我们将在下一章计算的 MTF，将这种空间信息压缩成频率响应。那很强大，但也隐藏了一些可见的结构。PSF 让我们直接看到模糊。

11. 包围能量：一个简单的 PSF 摘要

在转到 MTF 之前，引入一个简单的基于 PSF 的摘要：包围能量。

包围能量问的是：

PSF 总能量中有多少位于给定半径内？

这与 MTF 不同，但常常有用。更锐利的 PSF 将能量集中到更小的半径内。

这是一个以像素为单位简单实现：

def encircled_energy(psf):
    """
    Compute encircled energy as a function of radius in pixel units.

    Parameters
    ----------
    psf : 2D array
        Normalized PSF.

    Returns
    -------
    radii : 1D array
        Radius values in pixels.
    energy : 1D array
        Encircled energy for each radius.
    """
    n, m = psf.shape
    cy = (n - 1) / 2.0
    cx = (m - 1) / 2.0

    yy, xx = np.indices(psf.shape)
    r = np.sqrt((xx - cx)**2 + (yy - cy)**2)

    r_flat = r.ravel()
    psf_flat = psf.ravel()

    order = np.argsort(r_flat)
    r_sorted = r_flat[order]
    psf_sorted = psf_flat[order]

    cumulative = np.cumsum(psf_sorted)

    return r_sorted, cumulative

绘制完美对像差：

r_perfect, e_perfect = encircled_energy(perfect_psf)
r_aberr, e_aberr = encircled_energy(aberrated_psf)

plt.figure(figsize=(6, 4))
plt.plot(r_perfect, e_perfect, label="Perfect")
plt.plot(r_aberr, e_aberr, label="Aberrated")
plt.xlabel("Radius [pixels]")
plt.ylabel("Encircled energy")
plt.title("Encircled energy comparison")
plt.legend()
plt.grid(True)
plt.tight_layout()
plt.show()

如果像差 PSF 把能量向外扩散，它的包围能量曲线上升得更慢。

这给我们一个简单的说法：

像差系统需要一个更大的半径来包含相同比例的点源能量。

再次强调，这不是 MTF 的替代品。这是观察成像形成的另一种方式。

PSF 是空间的。包围能量总结空间集中度。MTF 将描述空间频率响应。

它们是关联的，但并非完全相同。

12. 斯特利尔比：中央峰作为警示灯

另一个常见的基于 PSF 的量是斯特利尔比。

在简单的计算术语中，斯特利尔比比较了像差 PSF 的峰值强度与衍射受限 PSF 的峰值强度，在匹配的归一化和采样下：

$$S = \frac{\max(\mathrm{PSF}$$\text{aberrated})} {\max(\mathrm{PSF}\text{diffraction-limited})}

代码：

def strehl_ratio(aberrated_psf, perfect_psf):
    return np.max(aberrated_psf) / np.max(perfect_psf)


strehl = strehl_ratio(aberrated_psf, perfect_psf)
print(strehl)

较低的斯特利尔比意味着中央峰相对于理想衍射受限情况损失了能量。

这是一个有用的警示灯，但它不是一个完整的像质描述。

两个 PSF 可以有相似的峰值但形状不同。一个可能有对称展宽，另一个可能有彗星状拖尾。它们的视觉效果可能不同。

所以规则是：

斯特利尔比总结中央峰。
PSF 图像显示整个模糊结构。

PSF 仍然值得观察。

13. 为什么点列图和 PSF 不同

在这一点上，值得将 PSF 与前面章节的点列图进行比较。

点列图来自追迹光线。它显示光线在哪与像平面相交。

PSF 来自衍射计算。它使用复光瞳函数，包括相位，并计算干涉。

对于严重像差的系统，它们可能讲述相似的故事，但它们不是相同的计算。

点列图回答：

采样光线落在了哪里？

PSF 回答：

一个点源在波干涉后产生怎样的强度分布？

对于较大的像差，几何点列图可以给出一个有意义的模糊第一印象。但在接近衍射极限时，点列图可能产生误导。它可能显示一个紧密的光线束，却未能显示衍射环。它也不能直接告诉我们能量在干涉作用下是如何分布的。

这就是光学软件提供多种图的原因：

• 点列图；
• 光线扇形图；
• OPD；
• PSF；
• MTF。

它们不是冗余的。它们是同一光学系统的不同投影。

完整的计算路径很重要：

光线追迹
→ OPD
→ 光瞳函数
→ PSF

点列图从光线交点分支出来。PSF 从波前相位分支出来。

这一区别是本书第四部分的核心。

14. 艾里斑与衍射极限

对于一个完美圆形孔径，解析的 PSF 就是艾里斑。其强度轮廓通常用贝塞尔函数表示。我们这里不需要完整的解析推导，但我们应该知道实际意义。

第一暗环的角半径近似为：

$$\theta \approx 1.22\frac{\lambda}{D}$$

其中：

• $\theta$ 是角半径，单位为弧度；
• $\lambda$ 是波长；
• $D$ 是孔径直径。

在焦平面，对应的半径近似为：

$$r \approx 1.22\lambda N$$

其中 $N$ 是 f 数。

这告诉我们一些物理上重要的事情：

更长的波长产生更大的衍射斑。
更小的孔径直径产生更大的衍射斑。
更大的 f 数产生更大的衍射斑。

这就是为什么缩小光圈可以改善像差但会使衍射模糊更糟。PSF 是由像差和孔径尺寸共同塑造的。

我们的 FFT 计算已经根据光瞳函数再现了艾里斑的结构。为了给输出采样分配物理单位，我们需要孔径大小、焦距、波长以及采样关系。这种完整的缩放在工程应用中很重要，但核心的生成机制已经可见：

清晰圆形光瞳
→ 傅里叶变换
→ 类艾里 PSF

这是我们先需要的计算事实。

15. 坐标缩放：FFT 输出的含义

FFT 返回频率域的采样。如果光瞳坐标被归一化，输出坐标也是归一化的倒易坐标。这对于比较形状没问题，但物理解释需要缩放。

对于一个在光瞳坐标 $X,Y$ 上采样的物理孔径场，焦平面坐标与衍射角有关。大略地：

$$\theta_x \sim \lambda f_x$$

$$\theta_y \sim \lambda f_y$$

其中 $f_x,f_y$ 是与光瞳平面采样相关的空间频率坐标。

在焦距为 $F$ 的透镜焦平面上：

$$X_\text{image} \sim F\theta_x$$

$$Y_\text{image} \sim F\theta_y$$

确切的实现取决于你如何采样物理光瞳以及如何定义傅里叶变换归一化。

对于本章，最安全的实用划分是：

首先计算出正确的归一化 PSF 形状。
然后使用已知的孔径采样、波长和焦距附加物理坐标。

试图同时做两件事常常让初学者困惑。

这里是一个使用物理光瞳直径的简单坐标辅助函数：

def psf_angular_coordinates(n_output, pupil_sample_spacing, wavelength):
    """
    Estimate angular coordinates for FFT-based Fraunhofer PSF.

    Parameters
    ----------
    n_output : int
        Number of samples in the padded pupil / PSF array.
    pupil_sample_spacing : float
        Physical pupil-plane sample spacing, in meters.
    wavelength : float
        Wavelength, in meters.

    Returns
    -------
    theta : 1D array
        Angular coordinate samples, in radians.
    """
    freq = np.fft.fftshift(
        np.fft.fftfreq(n_output, d=pupil_sample_spacing)
    )

    theta = wavelength * freq

    return theta

如果你的原始光瞳直径是 $D$，而未填充的光瞳网格在整个正方形坐标范围内有 $n$ 个采样，大致的采样间距是：

pupil_sample_spacing = D / n

然后：

D = 5e-3  # 5 mm aperture diameter
n_input = perfect_pupil.shape[0]
pad_factor = 4
n_output = n_input * pad_factor

dx_pupil = D / n_input
theta = psf_angular_coordinates(
    n_output=n_output,
    pupil_sample_spacing=dx_pupil,
    wavelength=wavelength
)

这给出了以弧度为单位的角采样。如果焦距已知，转换为焦平面长度：

focal_length = 50e-3  # 50 mm
image_coord = focal_length * theta

这种坐标缩放很有用，但也需要谨慎。随着你的模型越接近实际，你必须更仔细地定义：

• 光瞳直径；
• 入瞳 vs 出瞳；
• 焦距；
• 像空间折射率；
• 采样约定；
• FFT 归一化；
• 场点；
• 波长。

目前，记住主要的警告：

PSF 数组很容易计算。
物理坐标标签需要仔细的簿记。

第 14 章会更积极地回到采样和归一化。

16. 一个完整的最小 PSF 脚本

这是一个构建完美和像差 PSF 的紧凑脚本。

import numpy as np
import matplotlib.pyplot as plt


def make_pupil_grid(n=256):
    coord = np.linspace(-1.0, 1.0, n)
    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 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


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=2):
    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


def show_two_psfs(psf_a, psf_b, title_a, title_b, log_scale=False):
    if log_scale:
        image_a = np.log10(psf_a + 1e-12)
        image_b = np.log10(psf_b + 1e-12)
        label = "log10 intensity"
    else:
        image_a = psf_a
        image_b = psf_b
        label = "Normalized intensity"

    vmin = min(np.nanmin(image_a), np.nanmin(image_b))
    vmax = max(np.nanmax(image_a), np.nanmax(image_b))

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

    plt.subplot(1, 2, 1)
    plt.imshow(image_a, origin="lower", vmin=vmin, vmax=vmax)
    plt.colorbar(label=label)
    plt.title(title_a)
    plt.xlabel("Image sample x")
    plt.ylabel("Image sample y")

    plt.subplot(1, 2, 2)
    plt.imshow(image_b, origin="lower", vmin=vmin, vmax=vmax)
    plt.colorbar(label=label)
    plt.title(title_b)
    plt.xlabel("Image sample x")
    plt.ylabel("Image sample y")

    plt.tight_layout()
    plt.show()


wavelength = 550e-9

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

zero_opd = np.zeros_like(rho)
perfect_pupil, _ = pupil_function_from_opd(
    opd=zero_opd,
    mask=mask,
    wavelength=wavelength
)

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(perfect_pupil, pad_factor=4)
aberrated_psf, _ = compute_psf(aberrated_pupil, pad_factor=4)

show_two_psfs(
    perfect_psf,
    aberrated_psf,
    "Perfect PSF",
    "Aberrated PSF",
    log_scale=False
)

show_two_psfs(
    perfect_psf,
    aberrated_psf,
    "Perfect PSF",
    "Aberrated PSF",
    log_scale=True
)

strehl = np.max(aberrated_psf) / np.max(perfect_psf)
print("Approximate Strehl ratio:", strehl)

这个脚本仍然是一个教学模型。它不能替代完整的光学设计软件包。但它清晰地暴露了计算骨架：

制作光瞳
→ 从 OPD 添加相位
→ FFT
→ 幅度平方
→ 归一化
→ 对比

这正是我们想要打开的那个隐藏步骤。

17. 将此与 Optiland 连接

在一个真正的 Optiland 工作流中，你通常不会手动根据合成多项式创建 OPD 映射。你会定义或加载一个光学系统，选择一个场点和波长，然后请求分析工具提供波前或 PSF 结果。

概念上，流程是：

Optiland 光学系统
→ 光线追迹和波前分析
→ 光瞳函数或等效内部表示
→ PSF 分析

一个验证式的对比工作流看起来像这样：

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

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

# 2. 选择场点和波长。
# field = ...
# wavelength = ...

# 3. 使用 Optiland 分析工具计算波前 / OPD。
# wavefront = ...
# opd = wavefront.opd_map(...)

# 4. 从该 OPD 构建一个光瞳函数。
# pupil, phase = pupil_function_from_opd(opd, mask, wavelength)

# 5. 使用 NumPy 计算一个教学 PSF。
# psf_numpy, field_numpy = compute_psf(pupil, pad_factor=4)

# 6. 使用 Optiland 自身的 PSF 分析计算 PSF。
# psf_optiland = ...

# 7. 对比形状、归一化、采样和坐标缩放。

二十行 NumPy 代码不能取代一个成熟的光学设计工具。但它们让隐藏的计算变得可检视。

当 Optiland 给出一个 PSF 时，你现在应该能够问：

• 采样的是哪个光瞳？
• 使用了哪个波长？
• 使用了哪个场点？
• 包含了什么 OPD 或相位？
• 应用了什么孔径掩模？
• 使用了什么 FFT 采样和归一化？
• 图是线性还是对数的？
• 坐标是像平面、角度还是归一化的？

这些问题比记忆一个按钮顺序更有价值。

工具可以产生 PSF。理解意味着知道为了让那个 PSF 存在，必须发生什么。

18. PSF 计算中的常见错误

错误 1：直接对 OPD 进行傅里叶变换

OPD 不是光瞳函数。

错误：

psf = np.abs(np.fft.fft2(opd)) ** 2

正确：

phase = 2.0 * np.pi * opd / wavelength
pupil = amplitude * np.exp(1j * phase)
psf = np.abs(np.fft.fft2(pupil)) ** 2

FFT 作用于复场，而不是作为高度图的 OPD。

错误 2：忘记取幅度平方

傅里叶变换给出一个复场。PSF 是强度。

错误：

psf = np.fft.fft2(pupil)

正确：

field = np.fft.fft2(pupil)
psf = np.abs(field) ** 2

强度是幅度平方。

错误 3：比较未归一化的 PSF

如果两个 PSF 的总和不同，峰值对比可能会产生误导。

更好的是：

psf = psf / np.sum(psf)

然后你可以更公平地比较能量分布。

错误 4：相信零填充能提高分辨率

零填充让 PSF 显示更平滑。它并不提高物理光学分辨率。

错误 5：忽略显示尺度

线性图和对数图会让同一个 PSF 感觉非常不同。

在解读圆环、光晕或微弱拖尾之前，始终检查正在使用哪个尺度。

错误 6：忘记 PSF 取决于波长和场点

一个 PSF 不是一个透镜的通用属性。

它取决于：

• 波长；
• 场点；
• 聚焦位置；
• 孔径；
• 光瞳采样；
• 如果包括的话，偏振假设；
• 像差状态。

一个单独的 PSF 图只是更大的光学行为中的一个切片。

19. 本章为完整链添加了什么

我们现在可以扩展计算链：

透镜处方
→ 表面和材料
→ 光线追迹
→ 光程长度
→ OPD 映射
→ 光瞳振幅和相位
→ 复光瞳函数
→ 傅里叶变换
→ PSF

新的步骤是：

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

为实际比较应用了归一化。

这是一个重要的里程碑。

我们已经从几何光线走向了一个点的波动光学图像。我们也看到了为什么在一个真实的光学系统中，一个点永远不会保持一个点。

一个清晰的圆形光瞳产生衍射受限的类艾里图样。一个有像差的光瞳改变了孔径上的相位，并重新分配了 PSF 中的能量。一个有遮挡或变迹的光瞳改变了振幅，也改变了 PSF。

PSF 是孔径、衍射和像差在像空间交汇的地方。

下一章将进一步迈出一步。与其问一个点是如何模糊的，我们将问不同的空间频率是如何传输的。

这将我们从 PSF 带到 OTF，然后到 MTF。

在第 1 章看起来像软件输出的 MTF 曲线现在几乎伸手可及了。

本章摘要

PSF 是一个理想点源在像平面产生的强度分布。

在通常的标量夫琅禾费衍射近似下，复像平面场正比于复光瞳函数的傅里叶变换：

$$U(\xi,\eta) \propto \mathcal{F}{P(x,y)}$$

PSF 是该场的幅度平方：

$$\mathrm{PSF}(\xi,\eta) = |U(\xi,\eta)|^2$$

用最简 NumPy 代码表示：

field = np.fft.fftshift(np.fft.fft2(np.fft.ifftshift(pupil)))
psf = np.abs(field) ** 2
psf = psf / np.sum(psf)

一个完美的清晰圆形光瞳产生衍射受限的类艾里 PSF。一个有像差的光瞳将能量从理想的中央峰重新分配出去。一个改变了孔径振幅也会改变 PSF，即使波前相位是平坦的。

PSF 应被仔细检视，通常同时在线性和对数显示下。零填充可以让显示的 PSF 更平滑，但它并不提高实际的光学分辨率。

最重要的是，FFT 必须应用于复光瞳函数，而不是直接应用于 OPD 映射。

现在的计算链达到：

OPD
→ 光瞳函数
→ PSF

接下来，我们将问这个 PSF 如何对不同空间频率的图像细节作用。那个问题将引向 OTF 和 MTF。

生成的验证图

来源与验证说明

本章中的 PSF 计算使用适合教授光瞳到像计算的标量傅里叶光学模型。它有意省略了矢量衍射、偏振效应、部分相干、探测器采样以及详细的物理单位。这些省略是本章的边界，而不是声称这些效应不重要。