Revisiting PSF models: unifying framework and high-performance implementation

Basic Information

• Paper ID: 2502.03170
• Title: Revisiting PSF models: unifying framework and high-performance implementation
• Authors: Yan Liu, Vasiliki Stergiopoulou, Jonathan Chuah, Eric Bezzam, Gert-Jan Both, Michael Unser, Daniel Sage, Jonathan Dong
• Institutions: École Polytechnique Fédérale de Lausanne (EPFL), HHMI Janelia Research Campus
• Classification: physics.optics
• Publication Date: October 28, 2025 (arXiv v2)
• Paper Link: https://arxiv.org/abs/2502.03170

Abstract

Point spread function (PSF) models are core tools for localization microscopy. This paper proposes a unified theoretical framework for light propagation modeling in high numerical aperture (NA) imaging systems, proving the equivalence between methods based on 2D Fourier transforms and 1D Bessel integrals. The authors developed a high-performance open-source library based on PyTorch and conducted the first systematic benchmarking of different models. Experiments demonstrate that the Bessel strategy is optimal for axially symmetric beams, while Fourier methods are suitable for more general scenarios. The work supports efficient CPU/GPU computation and seamlessly integrates into deep learning and optimization pipelines.

Research Background and Motivation

1. Core Problem

The point spread function (PSF) is a fundamental concept in optical microscopy, describing the impulse response of an imaging system. In high numerical aperture (NA) imaging systems, accurate PSF models are critical for:

• Single-molecule localization microscopy (SMLM): Requires precise PSF for molecular localization
• Deconvolution microscopy: Requires accurate PSF models for image restoration
• PSF engineering: Achieves super-resolution imaging through specially designed PSFs

2. Problem Significance

• PSF modeling is central to computational imaging, directly affecting the performance of super-resolution microscopy (e.g., STED, MINFLUX)
• High-NA systems must account for complex physical effects including vectorial field characteristics, spherical aberration, and refractive index mismatch
• Accurate PSF models can improve 3D localization accuracy and achieve spatial resolution beyond the optical diffraction limit

3. Limitations of Existing Methods

• Theoretical level: Two mainstream methods exist (Fourier transform and Bessel function methods), but their relationship is unclear, lacking systematic comparison
• Implementation level: Existing libraries are mostly written in Java or MATLAB (e.g., the authors' previous PSF Generator), making integration with modern deep learning frameworks difficult
• Application level: Lack of systematic benchmarking of accuracy and computational speed, making it difficult for users to select appropriate models

4. Research Motivation

• Systematically revisit both methods starting from the Richards-Wolf integral
• Establish a unified framework proving the equivalence of both methods
• Provide high-performance PyTorch implementation supporting automatic differentiation and GPU acceleration
• Conduct comprehensive accuracy and speed benchmarking for the first time

Core Contributions

1. Unified Theoretical Framework: Starting from the Richards-Wolf integral, proves that Fourier (Cartesian) and Bessel (Spherical) methods are essentially different parameterizations of the same propagation integral
2. Universal Correction Factors: Systematically derives multiple physical correction factors (Gibson-Lanni spherical aberration, apodization, Fresnel transmission coefficients, arbitrary phase distortions, etc.) and applies them uniformly to both methods
3. High-Performance PyTorch Implementation: Develops the open-source library psf-generator implementing four propagators (scalar/vectorial × Cartesian/spherical), supporting CPU/GPU computation and automatic differentiation
4. Systematic Benchmarking: First comprehensive comparison of different PSF models for accuracy and computational speed, providing guidance for practical applications
5. Ecosystem Integration: Provides napari graphical interface plugin and chromatix optical simulation framework integration, promoting adoption by the open-source community

Methodology Details

Task Definition

Input:

• Physical parameters: numerical aperture (NA), wavelength (λ), refractive index (n), focal length (f)
• Incident field (pupil function): einc(s), defined in the pupil plane
• Spatial position: three-dimensional coordinates ρ = (x, y, z) near the focal point

Output:

• Electric field distribution E(ρ) in the focal region, i.e., the PSF

Constraints:

• Applicable to high-NA systems (must consider vectorial field effects)
• Must include various physical corrections (spherical aberration, apodization, transmission, etc.)

Theoretical Framework

1. Richards-Wolf Model Foundation

The starting point for all accurate PSF models is the Richards-Wolf integral (vectorial extension of the Debye integral):

$E(\rho) = -\frac{ifk}{2\pi} \int\int_{\Omega} e_{\infty}(s) \exp(iks \cdot \rho) d\Omega$

where:

• s = (sx, sy, sz): unit vector of incident ray direction
• k = 2πn/λ: wave number
• f: lens focal length
• Ω: solid angle region determined by NA

2. Two Parameterization Methods

Cartesian Parameterization (Fourier Method): Using (sx, sy) coordinates, dΩ = dsxdsy/sz

$E(\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$

This is essentially a 2D inverse Fourier transform, enabling efficient computation via FFT.

Spherical Parameterization (Bessel Method): Using spherical coordinates (θ, φ), dΩ = sinθdθdφ

$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$

For axially symmetric incident fields, the φ integral can be explicitly computed using Bessel function J0:

$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. Scalar vs. Vectorial Models

Scalar Model: Low-NA approximation, e∞(s) = einc(s)

Vectorial Model: Accounts for electric field vectorial properties, requiring basis transformation from cylindrical to spherical coordinates:

$e_{\infty}(\theta,\phi) = \text{[complex 3×2 matrix transformation]} \cdot [e_{inc}^x, e_{inc}^y]^T$

Including Fresnel transmission coefficients qs and qp, describing polarization-dependent transmission at interfaces.

Vectorial spherical parameterization simplifies to:

$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}$

where Ina(ρ,z) are one-dimensional integrals containing nth-order Bessel functions Jn.

Correction Factors

1. Gibson-Lanni Spherical Aberration

Describes spherical aberration caused by refractive index mismatch in stratified media:

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

where ts, tg, ti are the thicknesses of sample, glass coverslip, and immersion medium, respectively.

2. Apodization Factor

Ensures energy conservation during transformation from cylindrical to spherical coordinates:

$A(s) = \sqrt{\cos\theta}$

3. Gaussian Envelope

Describes non-ideal plane wave illumination:

$A(s) = \exp\left(-\frac{\sin^2\theta}{s_{env}^2}\right)$

4. Arbitrary Phase Distortion

Parameterized using Zernike polynomials:

$W(s) = \sum_{k=0}^{K-1} c_k Z_k(s) + W_0(s)$

All correction factors are uniformly expressed as:

$E(\rho) = -\frac{ifk}{2\pi} \int\int_{\Omega} a(s) \exp(iW(s)) e_{\infty}(s) \exp(iks \cdot \rho) d\Omega$

Technical Innovations

1. Equivalence Proof: First rigorous proof that Cartesian and spherical parameterizations are different representations of the same integral, eliminating long-standing theoretical ambiguity in the field
2. Universal Correction Factors: Extends correction factors previously applied only to Bessel methods (e.g., Gibson-Lanni, apodization) to Fourier methods
3. Custom FFT Implementation: Implements 2D FFT for arbitrary pixel sizes using chirp Z transform, addressing sampling issues with extremely small pixel sizes in localization microscopy
4. Efficient Numerical Integration: Spherical method employs Simpson's rule for 4th-order accuracy, with batch processing vectorization via torch.vmap
5. Differentiable Bessel Functions: Extends automatic differentiation capabilities for Bessel functions not natively supported by PyTorch

Experimental Setup

Datasets

This study primarily involves theoretical validation and algorithm benchmarking without real datasets. Test scenarios include:

• Analytical solution verification: Using Airy disk (analytical solution of circular aperture Fourier transform)
• Parameter settings: Wavelength 632 nm, NA range 0.5-1.3, refractive indices ns=1.3, ni=1.5, ng=1.5

Evaluation Metrics

1. Accuracy Assessment

• L2 error: δ = ||E - FAD||2, where FAD is the analytical Airy disk
• Convergence order: Analysis of theoretical convergence rate of numerical methods

2. Computational Efficiency

• Runtime: Time to generate a single 201×201 PSF image
• Parameter sweep: Pupil sizes from 32 to 1024 pixels
• Hardware platforms: Intel i9-10900X CPU and NVIDIA RTX 3090 GPU

Comparison Methods

Systematic comparison of four propagators:

1. ScalarCartesianPropagator: Scalar Cartesian method
2. ScalarSphericalPropagator: Scalar spherical method
3. VectorialCartesianPropagator: Vectorial Cartesian method
4. VectorialSphericalPropagator: Vectorial spherical method

Implementation Details

• Programming Framework: PyTorch 2.3+
• Numerical Integration: Composite Simpson's rule for spherical method
• FFT Implementation: Custom 2D FFT based on chirp Z transform
• Parallelization: Leveraging PyTorch's native parallelization and torch.vmap
• Data Format: Tensor shape (z, channel, x, y), channel=1 (scalar) or 3 (vectorial)

Experimental Results

Main Results

1. Accuracy Benchmarking (Figure 7)

Comparison with analytical Airy disk shows:

• Spherical method (Riemann rule): 1st-order convergence (O(h))
• Spherical method (Simpson rule): 4th-order convergence (O(h⁴))
• Cartesian method: Convergence rate between 1st and 2nd order
• Conclusion: Spherical method with Simpson rule achieves highest accuracy

2. Computational Speed Benchmarking (Figure 6)

CPU Performance:

• Small sizes (<512 pixels): Cartesian method faster
• Large sizes (>512 pixels): Spherical method faster
• Scalar vs. vectorial: Scalar 1.5× faster (Cartesian) or 3× faster (spherical)

GPU Performance:

• Cartesian method: Modest acceleration compared to CPU
• Spherical method: Flat curve showing significant GPU parallelization benefits
• Scalar vs. vectorial: Speed difference minimal, especially for Cartesian method

Key Findings:

• Vectorial model accuracy improvement far exceeds computational overhead
• Spherical method shows best scalability on GPU
• For medium-small image sizes, vectorial Cartesian method is the universal default choice

PSF Visualization (Figures 4-5)

Low-NA vs. High-NA Comparison (Figure 4)

• NA=0.5: Scalar and vectorial models produce similar results, classical Airy disk morphology
• NA=1.3: Vectorial model shows significant differences, energy dispersed to different field components, ring structure blurred

Phase Distortion Effects (Figure 5)

• Gibson-Lanni Correction: Refractive index mismatch causes spherical aberration, significantly degrading focal quality
• Donut PSF: Vortex phase mask produces ring-shaped PSF with zero center (used for STED)
• Half-moon PSF: π phase jump produces asymmetric PSF
• Astigmatism: Zernike polynomial-induced astigmatism shows elliptical rotation across different z-planes

Experimental Findings

1. Method Selection Guide:
   • Axially symmetric pupil function → Spherical method (high accuracy + GPU scalability)
   • Non-symmetric distortion/PSF engineering → Cartesian method (greater generality)
   • Default recommendation: Vectorial Cartesian method (generality + reasonable overhead)
2. Importance of Vectorial Effects: In high-NA systems, vectorial models are essential; scalar approximation produces significant errors
3. Impact of Correction Factors: Corrections like Gibson-Lanni spherical aberration are crucial for accurate modeling
4. Computational Complexity:
   • Cartesian: O(n log n), where n is lateral plane size
   • Spherical: O(n), where n is number of integration steps

Related Work

1. PSF Modeling Theory

• Classical Theory: Richards-Wolf integral (1959), Debye integral
• Aberration Models: Gibson-Lanni model (1991), Török extension (1997)
• Vectorial Models: Aguet (2009), Novotny & Hecht (2012)

2. Existing Software Implementations

• Commercial Software: Huygens PSF (widely used but closed-source)
• Java Implementation: Authors' previous PSF Generator plugin (2013)
• MATLAB Implementation: Nasse & Woehl (2010), Miora et al. (2024)
• Python Implementation: PyFocus (2022), SPITFIR(e) (2023)

3. Advantages Relative to This Work

• Theoretical level: First proof of equivalence between two methods, establishing unified framework
• Implementation level: Native PyTorch support, seamless deep learning integration
• Performance level: First systematic benchmarking, significant GPU acceleration
• Ecosystem: napari plugin, chromatix integration, promoting open-source adoption

4. Deep Learning Applications

Related work has applied PSF models to:

• Network Training: DeepSTORM (2020), DECODE (2021)
• Physics-Constrained Learning: Li et al. (2022) incorporating imaging process
• Self-Supervised Learning: Kobayashi et al. (2020) deconvolution
• Generative Models: FluoGAN (2023) for fluorescence image deconvolution

This paper's PyTorch implementation provides efficient, differentiable PSF generation tools for these applications.

Conclusions and Discussion

Main Conclusions

1. Theoretical Unification: Starting from the Richards-Wolf integral, proves that Fourier (Cartesian) and Bessel (Spherical) methods are different parameterizations of the same propagation integral, eliminating long-standing theoretical ambiguity
2. Method Selection:
   • Axially symmetric pupil: Spherical method (high accuracy, GPU scalability)
   • Non-symmetric distortion: Cartesian method (strong generality)
   • Universal recommendation: Vectorial Cartesian method (balanced performance and applicability)
3. Performance Advantages: Vectorial model accuracy improvement far exceeds computational overhead; should be prioritized in high-NA systems
4. Open-Source Contribution: Provides complete PyTorch library, napari plugin, and chromatix integration, supporting CPU/GPU computation and automatic differentiation

Limitations

1. Theoretical Assumptions:
   • Richards-Wolf integral has applicability conditions (Wolf & Li, 1981 discuss validity conditions)
   • Does not cover all light propagation models (e.g., complete Maxwell equation solutions)
2. Numerical Accuracy:
   • Cartesian method based on FFT has lower convergence rate than spherical Simpson method
   • Chirp Z transform increases computational overhead (3× FFT calls)
3. Application Scope:
   • Primarily targets fluorescence microscopy; other imaging modalities may require adjustments
   • Does not consider complex sample characteristics like scattering media
4. Validation Limitations:
   • Accuracy validation primarily based on analytical Airy disk
   • Lacks systematic comparison with experimentally measured PSFs

Future Directions

1. Deep Learning Integration:
   • Generate large-scale training datasets
   • Embed as differentiable layer in neural networks
   • Support end-to-end PSF engineering optimization
2. Extended Applications:
   • Virtual SMLM microscopy simulation
   • Single-particle tracking algorithm development
   • Adaptive optics system design
3. Performance Optimization:
   • Further optimize GPU parallelization strategies
   • Explore mixed-precision computation
   • Develop specialized hardware acceleration
4. Model Extensions:
   • Incorporate scattering effects
   • Support time-domain pulse propagation
   • Multi-color PSF modeling

In-Depth Evaluation

Strengths

1. Significant Theoretical Contribution:
   • First rigorous proof of equivalence between two mainstream PSF modeling methods, filling a field gap
   • Unified framework makes correction factor application more systematic and clear
   • Rigorous mathematical derivation with clear logic from Richards-Wolf integral to concrete implementation
2. High-Quality Implementation:
   • Native PyTorch implementation supporting automatic differentiation and GPU acceleration
   • Elegant code design with unified structure for four propagators
   • Custom FFT and differentiable Bessel functions demonstrate deep engineering expertise
3. Comprehensive Experiments:
   • First systematic benchmarking of accuracy and speed across different PSF models
   • Covers multiple NA, polarization, and distortion conditions
   • Provides clear method selection guidelines
4. Complete Open-Source Ecosystem:
   • Complete documentation and examples
   • napari plugin lowers usage barriers
   • chromatix integration promotes cross-platform applications
   • Adheres to FAIR principles, emphasizing reproducibility
5. Clear Writing:
   • Well-structured with clear hierarchical organization from background to theory to implementation
   • Rich figures (PSF visualization, benchmarking curves)
   • Sufficient but non-redundant technical details

Weaknesses

1. Limited Experimental Validation:
   • Accuracy validation primarily relies on analytical Airy disk, lacking comparison with real experimental PSFs
   • No direct comparison with other software (e.g., Huygens, Java PSF Generator)
   • Lacks application cases on real SMLM data
2. Expandable Theoretical Depth:
   • Equivalence proof is relatively intuitive (coordinate transformation), lacking deeper mathematical insights
   • Insufficient discussion of numerical stability and error propagation
   • Limited discussion of Richards-Wolf integral applicability boundaries
3. Insufficient Performance Analysis:
   • No memory consumption analysis
   • Lacks systematic testing across different hardware platforms
   • No performance comparison with other frameworks (JAX, TensorFlow)
4. Single Application Scenario:
   • Primarily focuses on fluorescence microscopy
   • Unexplored applications in other imaging modalities (holography, quantitative phase imaging)
   • Lacks complete cases of practical PSF engineering optimization
5. Incomplete Technical Details:
   • Simpson integration step size selection strategy not detailed
   • Insufficient discussion of FFT sampling theorem satisfaction conditions
   • Numerical stability of automatic differentiation not verified

Impact

1. Academic Value:
   • Unified framework will become standard reference for PSF modeling
   • Provides theoretical foundation for cross-disciplinary research between optical physics and computational imaging
   • Expected to be widely cited in super-resolution microscopy field
2. Practical Value:
   • Directly serves advanced microscopy techniques (SMLM, STED, MINFLUX)
   • Lowers barriers to deep learning applications in microscopy imaging
   • napari plugin enables non-specialist users to generate high-quality PSFs
3. Reproducibility:
   • Open-source code, detailed documentation, graphical interface form complete package
   • PyTorch ecosystem breadth ensures long-term maintainability
   • Adheres to FAIR principles, promoting scientific transparency
4. Field Advancement:
   • Paves way for physics-informed neural networks (PINN) applications in optics
   • May catalyze new PSF engineering optimization algorithms
   • Promotes open-source microscopy software ecosystem development

Applicable Scenarios

1. Direct Applications:
   • PSF modeling for single-molecule localization microscopy (SMLM)
   • System design for 3D deconvolution microscopy
   • PSF engineering (e.g., double helix PSF, donut PSF design)
2. Research Tools:
   • Training data generation for deep learning methods
   • Virtual microscopy simulation platform development
   • Optical system performance evaluation
3. Educational Use:
   • Visualization tool for optical imaging courses
   • Demonstration of high-NA optical effects
   • Teaching case for numerical method comparison
4. Extended Scenarios:
   • Wavefront sensing for adaptive optics
   • End-to-end optimization for computational optics
   • Joint modeling for multimodal imaging

Key References

1. Richards & Wolf (1959): Original Richards-Wolf integral, foundation of vectorial diffraction theory
2. Gibson & Lanni (1991): Gibson-Lanni spherical aberration model, stratified media modeling
3. Leutenegger et al. (2006): Fast focal field computation, modern Cartesian method implementation
4. Aguet (2009): Systematic vectorial PSF model paper (doctoral thesis)
5. Kirshner et al. (2013): Authors' previous PSF Generator Java version
6. Miora et al. (2024): Latest PSF computation method review and MATLAB implementation

Overall Assessment: This is an excellent paper combining theory and practice closely. Theoretically, it first proves the equivalence of two mainstream PSF modeling methods, establishing a unified framework; practically, it provides high-quality PyTorch implementation and systematic benchmarking. The paper's open-source spirit and ecosystem integration (napari, chromatix) reflect important contributions to the research community. Main limitations lie in lacking direct comparison with experimental data and existing software. This work will become an important tool in computational microscopy imaging, particularly with broad prospects in applications combining deep learning with physical models.