2025-11-29T14:34:18.811816

Revisiting PSF models: unifying framework and high-performance implementation

Liu, Stergiopoulou, Chuah et al.
Localization microscopy often relies on detailed models of point spread functions. For applications such as deconvolution or PSF engineering, accurate models for light propagation in imaging systems with high numerical aperture are required. Different models have been proposed based on 2D Fourier transforms or 1D Bessel integrals. The most precise ones combine a vectorial description of the electric field and precise aberration models. However, it may be unclear which model to choose, as there is no comprehensive comparison between the Fourier and Bessel approaches yet. Moreover, many existing libraries are written in Java (e.g. our previous PSF generator software) or MATLAB, which hinders the integration into deep learning algorithms. In this work, we start from the original Richards-Wolf integral and revisit both approaches in a systematic way. We present a unifying framework in which we prove the equivalence between the Fourier and Bessel strategies and detail a variety of correction factors applicable to both of them. Then, we provide a high-performance implementation of our theoretical framework in the form of an open-source library that is built on top of PyTorch, a popular library for deep learning. It enables us to benchmark the accuracy and computational speed of different models, thus allowing for an in-depth comparison of the existing models for the first time. We show that the Bessel strategy is optimal for axisymmetric beams while the Fourier approach can be applied to more general scenarios. Our work enables efficient PSF computation on CPU or GPU, which can then be included in simulation and optimization pipelines.
academic

Revisiting PSF models: unifying framework and high-performance implementation

基本信息

  • 论文ID: 2502.03170
  • 标题: Revisiting PSF models: unifying framework and high-performance implementation
  • 作者: Yan Liu, Vasiliki Stergiopoulou, Jonathan Chuah, Eric Bezzam, Gert-Jan Both, Michael Unser, Daniel Sage, Jonathan Dong
  • 机构: École Polytechnique Fédérale de Lausanne (EPFL), HHMI Janelia Research Campus
  • 分类: physics.optics
  • 发表时间: 2025年10月28日 (arXiv v2)
  • 论文链接: https://arxiv.org/abs/2502.03170

摘要

点扩散函数(PSF)模型是定位显微镜的核心工具。本文针对高数值孔径(NA)成像系统中的光传播建模问题，提出了统一的理论框架，证明了基于2D傅里叶变换和1D贝塞尔积分两类方法的等价性。研究开发了基于PyTorch的高性能开源库，首次实现了不同模型的系统性基准测试。实验表明，贝塞尔策略对轴对称光束最优，而傅里叶方法适用于更通用场景。该工作支持CPU/GPU高效计算，可无缝集成到深度学习和优化流程中。

研究背景与动机

1. 核心问题

点扩散函数(PSF)是光学显微镜的基础概念，描述了成像系统的脉冲响应。在高数值孔径(NA)成像系统中，精确的PSF模型对以下应用至关重要：

  • 单分子定位显微镜(SMLM)：需要精确的PSF进行分子定位
  • 去卷积显微镜：需要准确的PSF模型恢复图像
  • PSF工程：通过设计特殊PSF实现超分辨率成像

2. 问题的重要性

  • PSF建模是计算成像的核心，直接影响超分辨率显微镜(如STED、MINFLUX)的性能
  • 高NA系统需要考虑矢量场特性、球差、折射率失配等复杂物理效应
  • 精确的PSF模型可以提高3D定位精度，实现突破光学衍射极限的空间分辨率

3. 现有方法的局限性

  • 理论层面：存在两类主流方法（傅里叶变换方法和贝塞尔函数方法），但它们之间的关系不清晰，缺乏系统性比较
  • 实现层面：现有库多用Java或MATLAB编写（如作者之前的PSF Generator），难以集成到现代深度学习框架
  • 应用层面：缺乏准确性和计算速度的系统性基准测试，用户难以选择合适的模型

4. 研究动机

  • 从Richards-Wolf积分出发，系统性地重新审视两类方法
  • 建立统一框架，证明两类方法的等价性
  • 提供基于PyTorch的高性能实现，支持自动微分和GPU加速
  • 首次进行全面的准确性和速度基准测试

核心贡献

  1. 理论统一框架：从Richards-Wolf积分出发，证明了傅里叶(Cartesian)和贝塞尔(Spherical)方法本质上是同一传播积分的不同参数化形式
  2. 通用校正因子：系统性地推导了多种物理校正因子（Gibson-Lanni球差、消光因子、Fresnel透射系数、任意相位畸变等），并将其统一应用于两类方法
  3. 高性能PyTorch实现：开发了开源库psf-generator，实现了四种传播器（标量/矢量 × 笛卡尔/球面），支持CPU/GPU计算和自动微分
  4. 系统性基准测试：首次对不同PSF模型进行准确性和计算速度的全面比较，为实际应用提供选择指导
  5. 生态系统集成：提供napari图形界面插件和chromatix光学模拟框架集成，促进开源社区采用

方法详解

任务定义

输入：

  • 物理参数：数值孔径(NA)、波长(λ)、折射率(n)、焦距(f)
  • 入射场(pupil function)：einc(s)，定义在瞳孔平面
  • 空间位置：焦点附近的三维坐标ρ = (x, y, z)

输出：

  • 焦点区域的电场分布E(ρ)，即PSF

约束条件：

  • 适用于高NA系统（需考虑矢量场效应）
  • 需包含各种物理校正（球差、消光、透射等）

理论框架

1. Richards-Wolf模型基础

所有精确PSF模型的起点是Richards-Wolf积分（Debye积分的矢量扩展）：

E(ρ)=−ifk2π∫∫Ωe∞(s)exp⁡(iks⋅ρ)dΩE(\rho) = -\frac{ifk}{2\pi} \int\int_{\Omega} e_{\infty}(s) \exp(iks \cdot \rho) d\Omega

其中：

  • s = (sx, sy, sz)：入射光线方向的单位矢量
  • k = 2πn/λ：波数
  • f：透镜焦距
  • Ω：由NA决定的立体角区域

2. 两类参数化方法

笛卡尔参数化（傅里叶方法）： 使用(sx, sy)坐标，dΩ = dsxdsy/sz

E(ρ)=−ifk2π∫∫sx2+sy2≤smax2e∞(sx,sy)szexp⁡(ikszz)exp⁡(ik(sxx+syy))dsxdsyE(\rho) = -\frac{ifk}{2\pi} \int\int_{s_x^2+s_y^2 \leq s_{max}^2} \frac{e_{\infty}(s_x,s_y)}{s_z} \exp(iks_z z) \exp(ik(s_x x + s_y y)) ds_x ds_y

这本质上是二维逆傅里叶变换，可利用FFT高效计算。

球面参数化（贝塞尔方法）： 使用球坐标(θ, φ)，dΩ = sinθdθdφ

E(ρ)=−ifk2π∫0θmaxdθ∫02πdϕ e∞(θ,ϕ)exp⁡(ikρsin⁡θcos⁡(ϕ−φ))exp⁡(ikzcos⁡θ)sin⁡θE(\rho) = -\frac{ifk}{2\pi} \int_0^{\theta_{max}} d\theta \int_0^{2\pi} d\phi \, e_{\infty}(\theta,\phi) \exp(ik\rho\sin\theta\cos(\phi-\varphi)) \exp(ikz\cos\theta) \sin\theta

对于轴对称入射场，φ积分可用贝塞尔函数J0显式计算：

E(ρ)=−ifk∫0θmaxe∞(θ)J0(kρsin⁡θ)exp⁡(ikzcos⁡θ)sin⁡θ dθE(\rho) = -ifk \int_0^{\theta_{max}} e_{\infty}(\theta) J_0(k\rho\sin\theta) \exp(ikz\cos\theta) \sin\theta \, d\theta

3. 标量vs矢量模型

标量模型：低NA近似，e∞(s) = einc(s)

矢量模型：考虑电场矢量性质，需进行柱坐标到球坐标的基变换：

e∞(θ,ϕ)=[复杂的3×2矩阵变换]⋅[eincx,eincy]Te_{\infty}(\theta,\phi) = \text{[复杂的3×2矩阵变换]} \cdot [e_{inc}^x, e_{inc}^y]^T

包含Fresnel透射系数qs和qp，描述界面处的偏振依赖透射。

矢量球面参数化可简化为：

E(ρ)=−ifk2[I0x−I2xcos⁡2φ−I2ysin⁡2φ−I2xsin⁡2φ+I0y+I2ycos⁡2φ−2iI1xcos⁡φ−2iI1ysin⁡φ]E(\rho) = -\frac{ifk}{2} \begin{bmatrix} I_0^x - I_2^x\cos2\varphi - I_2^y\sin2\varphi \\ -I_2^x\sin2\varphi + I_0^y + I_2^y\cos2\varphi \\ -2iI_1^x\cos\varphi - 2iI_1^y\sin\varphi \end{bmatrix}

其中Ina(ρ,z)是包含n阶贝塞尔函数Jn的一维积分。

校正因子

1. Gibson-Lanni球差

描述分层介质中折射率失配引起的球差：

W(s)=2πλ(tsns2−ni2sin⁡2θ+tini2−ni2sin⁡2θ−ti∗ni∗2−ni2sin⁡2θ+...)W(s) = \frac{2\pi}{\lambda}\left(t_s\sqrt{n_s^2-n_i^2\sin^2\theta} + t_i\sqrt{n_i^2-n_i^2\sin^2\theta} - t_i^*\sqrt{n_i^{*2}-n_i^2\sin^2\theta} + ...\right)

其中ts、tg、ti分别是样品、玻璃盖片、浸没介质的厚度。

2. 消光因子(Apodization)

确保从柱坐标到球坐标变换时的能量守恒：

A(s)=cos⁡θA(s) = \sqrt{\cos\theta}

3. 高斯包络

描述非理想平面波照明：

A(s)=exp⁡(−sin⁡2θsenv2)A(s) = \exp\left(-\frac{\sin^2\theta}{s_{env}^2}\right)

4. 任意相位畸变

用Zernike多项式参数化：

W(s)=∑k=0K−1ckZk(s)+W0(s)W(s) = \sum_{k=0}^{K-1} c_k Z_k(s) + W_0(s)

所有校正因子统一表示为：

E(ρ)=−ifk2π∫∫Ωa(s)exp⁡(iW(s))e∞(s)exp⁡(iks⋅ρ)dΩE(\rho) = -\frac{ifk}{2\pi} \int\int_{\Omega} a(s) \exp(iW(s)) e_{\infty}(s) \exp(iks \cdot \rho) d\Omega

技术创新点

  1. 等价性证明：首次严格证明笛卡尔和球面参数化是同一积分的不同表示，消除了领域内长期存在的理论模糊性
  2. 校正因子通用化：将以往仅应用于贝塞尔方法的校正因子（如Gibson-Lanni、消光因子）推广到傅里叶方法
  3. 自定义FFT实现：基于chirp Z变换实现任意像素尺寸的2D FFT，解决定位显微镜中极小像素尺寸的采样问题
  4. 高效数值积分：球面方法采用Simpson规则实现4阶精度，并通过torch.vmap向量化批处理
  5. 可微分贝塞尔函数：扩展PyTorch原生不支持的贝塞尔函数自动微分功能

实验设置

数据集

本研究主要进行理论验证和算法基准测试，不涉及真实数据集。测试场景包括：

  • 解析解验证：使用Airy盘（圆孔径傅里叶变换的解析解）
  • 参数设置：波长632nm，NA范围0.5-1.3，折射率ns=1.3, ni=1.5, ng=1.5

评价指标

1. 准确性评估

  • L2误差：δ = ||E - FAD||2，其中FAD是解析Airy盘
  • 收敛阶：分析数值方法的理论收敛速率

2. 计算效率

  • 运行时间：生成单个201×201 PSF图像的时间
  • 参数扫描：瞳孔尺寸从32到1024像素
  • 硬件平台：Intel i9-10900X CPU和NVIDIA RTX 3090 GPU

对比方法

四种传播器的系统性比较：

  1. ScalarCartesianPropagator：标量笛卡尔方法
  2. ScalarSphericalPropagator：标量球面方法
  3. VectorialCartesianPropagator：矢量笛卡尔方法
  4. VectorialSphericalPropagator：矢量球面方法

实现细节

  • 编程框架：PyTorch 2.3+
  • 数值积分：球面方法使用复合Simpson规则
  • FFT实现：基于chirp Z变换的自定义2D FFT
  • 并行化：利用PyTorch的原生并行能力和torch.vmap
  • 数据格式：张量形状(z, channel, x, y)，channel=1(标量)或3(矢量)

实验结果

主要结果

1. 准确性基准（图7）

与解析Airy盘的比较显示：

  • 球面方法（Riemann规则）：1阶收敛率（O(h)）
  • 球面方法（Simpson规则）：4阶收敛率（O(h⁴)）
  • 笛卡尔方法：介于1-2阶之间的收敛率
  • 结论：球面方法使用Simpson规则具有最高精度

2. 计算速度基准（图6）

CPU性能：

  • 小尺寸(<512像素)：笛卡尔方法更快
  • 大尺寸(>512像素)：球面方法更快
  • 标量vs矢量：标量快1.5倍(笛卡尔)或3倍(球面)

GPU性能：

  • 笛卡尔方法：与CPU相比有适度加速
  • 球面方法：显示平坦曲线，受益于GPU并行化显著
  • 标量vs矢量：速度差异很小，尤其是笛卡尔方法

关键发现：

  • 矢量模型的精度提升远超其计算开销
  • GPU上球面方法的可扩展性最佳
  • 对于中小尺寸图像，矢量笛卡尔方法是通用默认选择

PSF可视化（图4-5）

低NA vs 高NA对比（图4）

  • NA=0.5：标量和矢量模型结果相似，经典Airy盘形态
  • NA=1.3：矢量模型显示明显差异，能量分散到不同场分量，环形结构模糊

相位畸变效果（图5）

  • Gibson-Lanni校正：折射率失配导致球差，显著降低焦点质量
  • 甜甜圈PSF：涡旋相位掩模产生中心为零的环形PSF（用于STED）
  • 半月PSF：π相位跳变产生非对称PSF
  • 散光：Zernike多项式引入的散光在不同z平面呈现椭圆旋转

实验发现

  1. 方法选择指南：
    • 轴对称瞳孔函数→球面方法（高精度+GPU可扩展）
    • 非对称畸变/PSF工程→笛卡尔方法（更通用）
    • 默认推荐：矢量笛卡尔方法（通用性+合理开销）
  2. 矢量效应的重要性：高NA系统中，矢量模型必不可少，标量近似产生显著误差
  3. 校正因子的影响：Gibson-Lanni球差等校正对精确建模至关重要
  4. 计算复杂度：
    • 笛卡尔：O(n log n)，n为横向平面尺寸
    • 球面：O(n)，n为积分步数

相关工作

1. PSF建模理论

  • 经典理论：Richards-Wolf积分(1959)、Debye积分
  • 球差模型：Gibson-Lanni模型(1991)、Török扩展(1997)
  • 矢量模型：Aguet(2009)、Novotny & Hecht(2012)

2. 现有软件实现

  • 商业软件：Huygens PSF（广泛使用但闭源）
  • Java实现：作者之前的PSF Generator插件(2013)
  • MATLAB实现：Nasse & Woehl(2010)、Miora等(2024)
  • Python实现：PyFocus(2022)、SPITFIR(e)(2023)

3. 本文相对优势

  • 理论层面：首次证明两类方法等价性，建立统一框架
  • 实现层面：PyTorch原生支持，无缝集成深度学习
  • 性能层面：首次系统性基准测试，GPU加速显著
  • 生态系统：napari插件、chromatix集成，促进开源采用

4. 深度学习应用

相关工作已将PSF模型用于：

  • 网络训练：DeepSTORM(2020)、DECODE(2021)
  • 物理约束学习：Li等(2022)结合成像过程
  • 自监督学习：Kobayashi等(2020)的去卷积
  • 生成模型：FluoGAN(2023)用于荧光图像去卷积

本文的PyTorch实现为这些应用提供了高效、可微分的PSF生成工具。

结论与讨论

主要结论

  1. 理论统一：从Richards-Wolf积分出发，证明了傅里叶(Cartesian)和贝塞尔(Spherical)方法是同一传播积分的不同参数化，消除了长期存在的理论模糊性
  2. 方法选择：
    • 轴对称瞳孔：球面方法（高精度、GPU可扩展）
    • 非对称畸变：笛卡尔方法（通用性强）
    • 通用推荐：矢量笛卡尔方法（平衡性能与适用性）
  3. 性能优势：矢量模型的精度提升远超计算开销，在高NA系统中应优先使用
  4. 开源贡献：提供完整的PyTorch库、napari插件和chromatix集成，支持CPU/GPU计算和自动微分

局限性

  1. 理论假设：
    • Richards-Wolf积分有其适用条件（Wolf & Li, 1981讨论了有效性条件）
    • 未涵盖所有光传播模型（如完整的Maxwell方程求解）
  2. 数值精度：
    • 笛卡尔方法基于FFT，收敛率低于球面Simpson方法
    • chirp Z变换增加了计算开销（3倍FFT调用）
  3. 应用范围：
    • 主要针对荧光显微镜，其他成像模态可能需要调整
    • 未考虑散射介质等复杂样品特性
  4. 验证局限：
    • 准确性验证主要基于解析Airy盘
    • 缺乏与实验测量PSF的系统性比较

未来方向

  1. 深度学习集成：
    • 用于生成大规模训练数据集
    • 作为可微分层嵌入神经网络
    • 支持端到端PSF工程优化
  2. 扩展应用：
    • 虚拟SMLM显微镜仿真
    • 单粒子追踪算法开发
    • 自适应光学系统设计
  3. 性能优化：
    • 进一步优化GPU并行策略
    • 探索混合精度计算
    • 开发专用硬件加速
  4. 模型扩展：
    • 纳入散射效应
    • 支持时域脉冲传播
    • 多色PSF建模

深度评价

优点

  1. 理论贡献显著：
    • 首次严格证明两类主流方法的等价性，填补了领域空白
    • 统一框架使校正因子的应用更加系统和清晰
    • 数学推导严谨，从Richards-Wolf积分到具体实现逻辑清晰
  2. 实现质量高：
    • PyTorch原生实现，支持自动微分和GPU加速
    • 代码设计优雅，四种传播器结构统一
    • 自定义FFT和可微分贝塞尔函数体现深厚的工程能力
  3. 实验全面：
    • 首次系统性基准测试准确性和速度
    • 覆盖多种NA、偏振、畸变条件
    • 提供清晰的方法选择指南
  4. 开源生态完善：
    • 完整的文档和示例
    • napari插件降低使用门槛
    • chromatix集成促进跨平台应用
    • 遵循FAIR原则，强调可复现性
  5. 写作清晰：
    • 结构合理，从背景到理论到实现层次分明
    • 图表丰富（PSF可视化、基准测试曲线）
    • 技术细节充分但不冗余

不足

  1. 实验验证有限：
    • 准确性验证主要依赖解析Airy盘，缺乏与真实实验PSF的对比
    • 未提供与其他软件（如Huygens、PSF Generator Java版）的直接比较
    • 缺少在真实SMLM数据上的应用案例
  2. 理论深度可扩展：
    • 等价性证明较直观（坐标变换），缺少更深层的数学洞察
    • 未讨论数值稳定性和误差传播
    • 对Richards-Wolf积分的适用边界讨论不足
  3. 性能分析不够细致：
    • 未分析内存占用
    • 缺少不同硬件平台的系统性测试
    • 未提供与JAX、TensorFlow等其他框架的性能对比
  4. 应用场景单一：
    • 主要关注荧光显微镜
    • 未探讨在其他成像模态（如全息、定量相位成像）的应用
    • 缺少实际PSF工程优化的完整案例
  5. 技术细节待补充：
    • Simpson积分的步长选择策略未详细说明
    • FFT采样定理的满足条件讨论不足
    • 自动微分的数值稳定性未验证

影响力

  1. 学术价值：
    • 统一框架将成为PSF建模的标准参考
    • 为光学物理和计算成像的交叉研究提供理论基础
    • 预计将被超分辨率显微镜领域广泛引用
  2. 实用价值：
    • 直接服务于SMLM、STED、MINFLUX等先进显微技术
    • 降低深度学习在显微成像中的应用门槛
    • napari插件使非专业用户也能生成高质量PSF
  3. 可复现性：
    • 开源代码、详细文档、图形界面三位一体
    • PyTorch生态系统的广泛性保证长期可维护性
    • 遵循FAIR原则，促进科学透明度
  4. 领域推动：
    • 为物理信息神经网络(PINN)在光学中的应用铺平道路
    • 可能催生新的PSF工程优化算法
    • 促进开源显微镜软件生态的发展

适用场景

  1. 直接应用：
    • 单分子定位显微镜(SMLM)的PSF建模
    • 3D去卷积显微镜的系统设计
    • PSF工程（如双螺旋PSF、甜甜圈PSF设计）
  2. 研究工具：
    • 深度学习方法的训练数据生成
    • 虚拟显微镜仿真平台开发
    • 光学系统性能评估
  3. 教育用途：
    • 光学成像课程的可视化工具
    • 高NA光学效应的演示
    • 数值方法比较的教学案例
  4. 扩展场景：
    • 自适应光学的波前传感
    • 计算光学的端到端优化
    • 多模态成像的联合建模

参考文献（关键文献）

  1. Richards & Wolf (1959): 原始Richards-Wolf积分，矢量衍射理论基础
  2. Gibson & Lanni (1991): Gibson-Lanni球差模型，分层介质建模
  3. Leutenegger et al. (2006): 快速焦场计算，笛卡尔方法的现代实现
  4. Aguet (2009): 矢量PSF模型的系统性论文（博士论文）
  5. Kirshner et al. (2013): 作者之前的PSF Generator Java版
  6. Miora et al. (2024): 最新的PSF计算方法综述和MATLAB实现

总体评价：这是一篇理论与实践结合紧密的优秀论文。理论上首次证明了两类主流PSF建模方法的等价性，建立了统一框架；实践上提供了高质量的PyTorch实现和系统性基准测试。论文的开源精神和生态系统集成（napari、chromatix）体现了对科研社区的重要贡献。主要不足在于缺乏与实验数据和现有软件的直接对比。该工作将成为计算显微成像领域的重要工具，特别是在深度学习与物理模型结合的应用中具有广阔前景。