# Iterative phase retrieval algorithm for space-variant PSF in optical systems with aberrations

Dylan Jonathan Brault, Corinne Fournier, Tatiana Latychevskaia

► To cite this version:

Dylan Jonathan Brault, Corinne Fournier, Tatiana Latychevskaia. Iterative phase retrieval algorithm for space-variant PSF in optical systems with aberrations. Optics Letters, 2025, 50 (6), pp.1767-1770. 10.1364/OL.550272 . hal-05002576

HAL Id: hal-05002576

https://hal.science/hal-05002576v1

Submitted on 24 Mar 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# ITERATIVE PHASE RETRIEVAL ALGORITHM FOR SPACE-VARIANT PSF IN OPTICAL SYSTEMS WITH ABERRATIONS

Dylan Brault ID 1,*, Corinne Fournier ID 1, and Tatiana Latychevskaia ID 2,3

¹Université Jean Monnet Saint-Etienne, CNRS, Institut d'Optique Graduate School, Laboratoire Hubert Curien UMR 5516, 42023, Saint-Etienne, France

²Center for Life Sciences, Paul Scherrer Institute, Forschungsstrasse 111, 5232 Villigen, Switzerland.

³Department of Physics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland

*dylan.brault@telecom-st-etienne.fr

## ABSTRACT

Iterative phase retrieval algorithms are widely used in digital optics for their efficiency and simplicity. Conventionally, these algorithms do not consider aberrations as they assume an ideal, aberration-free optical system. Here, we propose modified iterative phase retrieval algorithms that take into account the space-invariant and space-variant point spread function of the optical system.

## 1 INTRODUCTION

Iterative Phase Retrieval (IPR) algorithms [1, 2] allow reconstruction of the complete complex-valued wavefront from one or more intensity measurements by using back-and-forth wavefront propagation between the detector and sample's planes and applying constraints $\gamma$ based on the object's properties (most commonly, finite size and positive absorption). Originally, IPR algorithms were applied for the reconstruction of objects from X-ray diffraction patterns, for no lenses are available for X-rays [3]. Later, the IPR approach was generalized to digital in-line holography, where it allowed to solve the twin image problem [4]. In the case of in-line holography, the wavefront propagation is often calculated through a convolution with a free-space propagation function. A hologram can be obtained using either a lens-free system, as in the case of in-line, Gabor type or low-energy electron holography [5, 6], or by using an optical system with lenses, as in the case of digital holographic microscopy imaging [7]. Since there are no perfect lenses, such optical systems always have aberrations. In this study, we show how an iterative phase retrieval algorithm can be modified by including aberration correction directly into the algorithm. The proposed method allows considering both space-invariant and space-variant Point Spread Functions (PSF). We describe the principle and the algorithms, and as an example of practical application, we demonstrate the results obtained by applying a modified IPR algorithm for light optical holograms. The proposed here approach can be applied for optical systems with aberrations and/or space-variant PSF. It could be potentially useful for application in electron microscopy, where the resolution is mainly limited by the lens aberrations.

## 2 HOLOGRAM FORMATION MODELS

In this section, the conventional model of wavefront propagation through an ideal optical system is first recalled. Then, a wavefront propagation model that accounts for space-invariant aberration is introduced. Finally, we present a method to propagate a wavefront numerically with space-variant aberrations at a low computational cost.

### 2.1 Hologram formation model without aberration

Reconstruction algorithms aim to retrieve the complex-valued transmission function that describes how the incident wave and the sample interact. The transmission function is typically modeled as a 2D discrete complex-valued distribution, with its modulus, $a$, describing the absorption properties of the sample and phase $\phi$ describing the phase delay added by the sample to the incident wavefront:

$$ \underline{t}(x, y) = a(x, y) \exp(i\phi(x, y)) \quad (1) $$

where $(x, y)$ are the coordinates in the sample plane and the underlying values are complex-valued distributions. When a plane wave of unitary amplitude illuminates the sample, in the absence of noise, the data $\underline{d}$ recorded on the sensor plane, the so-called hologram, corresponds to the intensity of the interference pattern of the sample at a distance $z$ and is modeled by:

$$ d(x', y') = |\underline{U}_{\text{diff}}(x', y')|^2 = |\underline{t}_{x', y'} \underline{h}_{z}|^2 \quad (2) $$

where $\underline{U}_{\text{diff}}$ is the discrete complex amplitude of the diffracted wavefront, $\underline{t}$ is the transmission function, $\underline{h}$ is a discrete propagation kernel for a distance $z$, $(x', y')$ are the coordinates in the hologram plane and $*$ is the discrete convolution operator. For light optical holograms, Fresnel and Rayleigh-Sommerfeld kernels are commonly used to simulate the wavefront propagation [8]. This hologram formation model describes the propagation in a perfect optical system and is the core of most reconstruction algorithms. However, slight misalignment, tilts of the sample, nonstandard uses of microscope objectives, or poorly corrected objectives lead to aberrations that must be considered to perform morphologically reliable and quantitative sample reconstructions.

### 2.2 Model of propagation with a Space-Invariant PSF (SI-PSF)

For PSF that do not vary in the field of view (FOV), wavefront propagation can still be modeled as a convolution between the transmission function $\underline{t}$ and a new aberrated propagation kernel $\underline{h}^{\text{SI}}$. This aberrated propagation kernel is given by the convolution between the aberration-free propagation kernel and the
complex-valued PSF of the optical system $\underline{\underline{p}}$ for $z=0$:

$$ \underline{\underline{h}}^{\text{SI}}(x', y') = \underline{\underline{h}}_{z} * \underline{\underline{p}}_{x', y'} \quad (3) $$

$\underline{\underline{p}}$ for $z=0$ can, for example, be estimated with model fitting approaches like proposed in [9]. Computing the propagation of a wavefront with a space-invariant aberration has the same computational cost as conventional aberration-free propagation, because the wavefront propagation model remains a simple convolution.

## 2.3 Model of propagation with a Space-Variant PSF (SV-PSF)

For PSF that vary in the FOV, one has to interpolate the complex-valued PSF to compute the image formation model. Naive interpolation of the PSF is prohibitive for computational time reasons as it would require a convolution with a local propagation kernel for each image pixel. Several techniques have been developed to reduce the computational time of PSF interpolation [10, 11, 12, 13]. Most of these techniques have been applied to compute space-variant aberration models with real-valued PSF and real-valued objects. Also, contrary to holography applications, which require extended propagation kernels, these methods have been mostly used for PSF located on a small spatial support. We chose to adapt the method proposed by Flicker and Rigaut [11] to interpolate the PSF.

This method relies on a modal decomposition of the space-varying PSF $\underline{\underline{h}}^{\text{SV}}$:

$$ \underline{\underline{h}}^{\text{SV}}(x, y, x', y') \approx \sum_{p=1}^{P} \underline{\underline{w}}_{p}(x, y) \underline{\underline{m}}_{p}(x', y') \quad (4) $$

where $P$ is the number of spatial locations where the PSF has been estimated, $\underline{\underline{m}} = \{\underline{\underline{m}}_p\}_{p=1,...,P}$ is a set of space-invariant orthogonal components of the PSF (modes), and $\underline{\underline{w}} = \{\underline{\underline{w}}_p\}_{p=1,...,P}$ is a complex-valued weighting map that describes how the mode $\underline{\underline{m}}_p$ should be weighted in the FOV. This separable approximation simplifies the image formation model such that the complex amplitude in the detector plane $\underline{\underline{U}}_{\text{diff}}$ can be expressed as a sum of convolutions:

$$ \underline{\underline{U}}_{\text{diff}}(x', y') = \underline{\underline{h}}^{\text{SV}} * \sum_{(x', y')}^{\text{SV}} \underline{\underline{t}} = \sum_{p=1}^{P} \left[ \underline{\underline{w}}_{p} \odot \underline{\underline{t}} \right] * \underline{\underline{m}}_{p} \quad (5) $$

where $*^{\text{SV}}$ symbolizes the space-variant convolution and $\odot$ is the pixel-wise product.

In practice, only a few modes are required when $\underline{\underline{h}}^{\text{SV}}$ slowly varies in the FOV. Thus, the separable approximation of Eq. 4 can be performed using a low-rank approximation (see Fig. S2) as proposed in [11]. The matrix $\underline{\underline{K}}$, whose columns contain each estimated PSF, is decomposed using an Singular Value Decomposition (SVD):

$$ \underline{\underline{K}} = \underline{\underline{U}} \cdot \underline{\underline{\Sigma}} \cdot \underline{\underline{V}}^{\text{T}} = \sum_{p}^{P} \underline{\underline{u}}_{p} \cdot \underline{\underline{\sigma}}_{p} \underline{\underline{v}}_{p}^{\text{T}} \approx \sum_{p=1}^{Q} \underline{\underline{u}}_{p} \cdot \underline{\underline{\sigma}}_{p} \underline{\underline{v}}_{p}^{\text{T}} \quad (6) $$

where $\underline{\underline{u}}_p$ and $\underline{\underline{v}}_p$ are the $p$-th left and right singular vectors, $\underline{\underline{\sigma}}_p$ the corresponding singular value and $\cdot$ is the matrix multiplication. By analogy with Eq. 4, $\underline{\underline{U}}$ is a matrix that contains the

orthogonal modes $\underline{\underline{m}}$ on which the PSF are decomposed and $\underline{\underline{\Sigma}} \cdot \underline{\underline{V}}^{\text{T}}$ contains the corresponding weights $\underline{\underline{w}}$. As the singular value $\underline{\underline{\Sigma}}$ decreases quickly, a low-rank approximation can be performed with $Q < P$. Moreover, considering the complex-valued PSF is slowly varying in the FOV, instead of estimating the PSF on each pixel, only an estimation on a coarse regular grid is required. Then, a smooth interpolation of weights $\underline{\underline{w}}$ on the pixel grid is performed such that each $\underline{\underline{w}}_p$ has the same size as $\underline{\underline{t}}$. We chose to interpolate these weights with cubic functions $\varphi$. Note that in the following sections, the back-propagation kernel associated with $\underline{\underline{h}}$ is noted $\underline{\underline{h}}^{\text{T}}$.

# 3 GENERALIZATION OF THE IPR ALGORITHM FOR ABERRATED SYSTEMS

In the following, we describe how to adapt the IPR algorithm for optical systems with SI-PSF and SV-PSF. Alg. S1 reminds the pseudo code of the IPR algorithm [4], and Fig. S1 depicts a general scheme of this approach.

## 3.1 Iterative phase retrieval algorithm for SI-PSF

In the case of SI-PSF, the IPR algorithm is modified by replacing the aberration-free propagator with $\underline{\underline{h}}^{\text{SI}}$. The backpropagation step is performed by convolving the complex amplitude $\underline{\underline{U}}_{\text{diff}}$ with $\underline{\underline{h}}^{\text{SI}\dagger}$ where $\underline{\underline{h}}^{\text{SI}\dagger}$ is obtained by conjugating $\underline{\underline{h}}^{\text{SI}}$ in the Fourier domain such that convolving $\underline{\underline{h}}^{\text{SI}}$ with $\underline{\underline{h}}^{\text{SI}\dagger}$ is approximately an identity operation (sampling and truncation of the FOV aside). The IPR algorithm in this configuration is very similar to the classical IPR algorithm [4]. Thus, after calibrating the PSF, the computational time needed for the reconstruction is the same as that required for the classical IPR algorithm.

## 3.2 Iterative phase retrieval algorithm for SV-PSF

To account for SV-PSF in the IPR algorithm, we propose to exploit the image formation model described in Sec. 2.2.3 to perform the propagation steps. The backpropagation step is performed by conjugating the modes $\underline{\underline{m}}$ in the Fourier domain and by conjugating the weights $\underline{\underline{w}}$. Alg. 1 provides a pseudo code for the IPR algorithm considering SV-PSF. Following Eq. 5, the propagation steps are performed as the sum of weighted convolution instead of a simple convolution. The computational time of the IPR algorithm is then proportional to the rank $Q$ of the low-rank approximation described in Eq. 6.

# 4 RESULTS

The proposed approach has been applied to simulated and experimental light optical holograms. In the results presented below, the propagation was simulated using Rayleigh-Sommerfeld approximation [8].

## 4.1 Reconstruction results on simulated holograms

To test the robustness of our approach, two holograms where simulated with the model for SV-PSF described in section 2.2.3. The first hologram was simulated with a slowly spatially varying
**Algorithm 1 IPR algorithm for SV-PSF**

1: **Input:** Data **d** (L × C image)

2: **Input:** Propagation kernels matrix **K** (LC × P matrix)

3: **Input:** Number of iterations $N_{iter}$ (integer)

4: **Input:** Rank of the low rank approximation **Q** (integer)

5: **Output:** Transmission function **t** (L × C image)

**Computation of modes and weights**

6: $[\underline{U}, \underline{\Sigma}, \underline{V}] \leftarrow \text{SVD}(\underline{K}) \triangleright \text{Low rank approximation of } \underline{K}$

7: $\underline{w} \leftarrow \varphi(\underline{\Sigma} \cdot \underline{V}^T) \triangleright \text{Interpolating weights on the sensor grid}$

8: $\underline{m} \leftarrow \text{reshape}(\underline{U}) \triangleright \text{Reshaping } \underline{U}(\rightarrow L \times C \times P \text{ 3D matrix})$

**IPR algorithm for SI-PSF**

9: $\underline{U}_{\text{diff}} \leftarrow \sqrt{\underline{d}} \triangleright \text{Initialization of the wavefront on the sensor}$

10:

11: **for** $n = 1$ **to** $N_{iter}$ **do**

12: $\underline{t} \leftarrow \sum_{p=1}^{Q} (\underline{w}_{p}^{*} \odot \underline{U}_{\text{diff}}) * \underline{m}_{p}^{\dagger} \triangleright \text{Back-propagation of } \underline{U}_{\text{diff}}$

13:

14: **for** $k = 1$ **to** $L \times C$ **do**

15: **if** $\underline{t}_k \notin \gamma$ **then** ▶ Applying constraints

16: $\underline{t}_k = 1$

17:

18: $\underline{U}_{\text{diff}} \leftarrow \sum_{p=1}^{Q} (\underline{w}_p \odot \underline{t}) * \underline{m}_p \triangleright \text{Propagation of } \underline{t}$

19:

20: $\underline{U}_{\text{diff}} \leftarrow \sqrt{\underline{d}} \times \exp(i \times \arg(\underline{U}_{\text{diff}})) \triangleright \text{Updating modulus}$

21:

22: **return** $\underline{t}$

PSF and the second one with a rapidly varying PSF. These two SV-PSF holograms are reconstructed under three described model hypothesis.

A pure phase binary object with maximal phase shift of 1 radian was used in the simulations (see Fig.1). To perform realistic hologram simulations, the experimental PSF of the hologram described in section 4.4.2 was used. The PSF (defocus and aberrations) was estimated at $P = 49$ positions (7 × 7 regular grid) in the FOV using the method described in [9]. Typical Zernike coefficients are provided in Supplementary Information (Table S2). As this method estimates the PSF on an irregular grid, interpolating the estimated Zernike coefficients on the regular grid has been performed using LOWESS interpolation [14]. The simulation parameters (matching the experimental ones) are provided in Table S1 and the pupil functions are displayed in Fig. S3. In this case, aberrations are slowly varying, furthermore, the sample slide is slightly tilted (0.5° tilt). Simulations with these experimental PSFs are referred to as 'Slowly varying SV-PSF' simulations in the following. To study the efficiency of our approach, we performed reconstruction with a more rapidly varying SV-PSF that was achieved by tilting the sample slide of 3.5°. Even though this tilt introduces a linear spatial variation in $z$, the PSF, which is a Rayleigh-Sommerfeld kernel, does not evolve linearly. Simulations with these PSFs are referred to as 'Rapidly varying SV-PSF' simulations in the following. The classical IPR algorithm reconstructions have been performed using a distance measured with a sharpness autofocusing criterion [15], so that the evaluated reconstruction errors do not come

from a defocus error due to spherical aberrations. As the phase object introduces a positive phase shift, the sample constraint $\gamma$ has been chosen such that the reconstructed object should be a pure phase with a positive phase shift. To achieve convergence of all the proposed reconstruction methods, we set the number of iterations $N_{iter} = 200$. Fig. 1.(a) shows the reconstructions of the simulated SV-PSF holograms with the classical IPR algorithm and accounting for SI-PSF and SV-PSF presented in section 3. To provide a quantitative comparison the root mean square error (RMSE) and the structural similarity index measure (SSIM) were estimated between the reconstructed distributions and the ground truth. Fig. 1.(a) shows that the classical IPR

<figure><img src="image_2.png" /><figcaption>Figure 1: Phase distributions reconstructed from the simulated hologram (a) with classical IPR algorithm, (b) with SI-PSF IPR algorithm, and (c) with SV-PSF IPR algorithm ($P = 49$, $Q = 15$)</figcaption></figure>

algorithm fails to reduce the twin image in both simulated cases, leading to high error values and strong reconstruction artifacts. Compared to the other methods, the classical IPR algorithm provides the worst RMSE and SSIM. As expected, considering an SI-PSF (Fig. 1.(b)) with aberrations estimated at the center of the FOV improves the reconstruction quality for slowly varying PSF. However, when the PSF variation increases, the SI-PSF hypothesis is insufficient to reduce the twin image artifacts. As expected, considering the SV-PSF (Fig. 1.(c)), the quality of the reconstructions regarding RMSE and SSIM has increased significantly. The twin image artifacts are almost entirely removed in that case, leading to more accurate reconstructions. Only $Q = 15$ modes have been used to model the SV-PSF in these reconstructions. However, let us notice that the reconstruction quality could be further improved by increasing $Q$. As computational time is the bottleneck of our SV-PSF method, we proposed reconstructing the slowly varying SV-PSF simulated hologram with a limited convolution budget. Thus, we chose $n_{iter}$ and $Q$ such that the algorithm should stop after 400 convolutions. We compared the reconstruction results with the classical IPR algorithm and SI-PSF IPR algorithm with 400 convolutions (see Table 1). As reported in Table 1, the SV-PSF IPR algorithm provides more quantitatively correct reconstructions, even with a reduced number of iterations, than classical IPR and SI-PSF algorithms. Thus, for similar computational time, and if the PSF is calibrated in the FOV, SV-PSF algorithm provides the most quantitative reconstruction of the three methods.

## 4.2 Reconstruction results on experimental hologram

A similar reconstruction comparison has been performed on experimental hologram. These data have been acquired with the
<table><caption>Table 1: Reconstruction quality criterion with the 3 hologram formation models with limited convolution budget (400 convolutions). "conv." stands for "convolution".</caption><thead><tr><th></th><th>IPR</th><th>SI-PSF IPR</th><th colspan="3">SV-PSF IPR</th></tr></thead><tbody><tr><th scope="row">N<sub>iter</sub></th><td>200</td><td>200</td><td>40</td><td>20</td><td>10</td></tr><tr><th scope="row">Q</th><td></td><td></td><td>5</td><td>10</td><td>20</td></tr><tr><th scope="row"># of conv.</th><td>400</td><td>400</td><td>400</td><td>400</td><td>400</td></tr><tr><th scope="row">RMSE</th><td>0.0755</td><td>0.0267</td><td>0.0126</td><td>0.0114</td><td>0.0206</td></tr><tr><th scope="row">SSIM</th><td>0.1132</td><td>0.2943</td><td>0.9853</td><td>0.9827</td><td>0.9457</td></tr></tbody></table>

<figure><img src="image_3.png" /><figcaption>Figure 2: Phase of the reconstruction on three regions in the FOV obtained (a) with classical IPR algorithm, with IPR algorithm (b) with SI-PSF, and (c) with SV-PSF ($P = 28, Q = 15$)</figcaption></figure>

setup described in [9]. The sample consists of calibration latex beads on a microscope slide. The beads are monodispersed of diameter of $1.0 \pm 0.06 \ \mu\text{m}$ and refractive index of $1.58$. The reconstruction of the beads are expected to be similar, regardless of their position in the FOV. Space-variant aberrations have been added to the hologram recording by choosing a wrong correction collar tuning (0.21 mm instead of 0.17 mm). Experimental parameters for the acquisition are provided in Table S1 and pupil functions are displayed in Fig.S3. The holograms have been reconstructed on the $204 \times 273 \ \mu\text{m}^2$ FOV to highlight the effect of space-varying PSFs. Indeed, for the selected large FOV, both aberrations and defocus exhibit significant variations (Fig. S3). Fig. 2 shows the reconstruction results for the three described methods. The classical IPR algorithm (Fig. 2.(a)) poorly succeeds in removing twin-image artifacts for such an aberrated system. The radial symmetry of the reconstructed distribution, expected for these spherical objects, is broken. As seen in region 3, as the PSF varies, the morphological and refractive properties of the object are modified. This can lead to misinterpretation of the reconstruction in other imaging contexts like microbiological ones. Using a SI-PSF hypothesis (Fig. 2.(b)) slightly improves the reconstruction results. This can be seen in region 2, which is almost at the center of the FOV, where the SI-PSF has been calibrated. Indeed, this reconstruction recovers the expected radial symmetry of the beads. However, this method fails to reconstruct the beads reproducibly when they are located far from the center of the FOV. It should be noted that strong twin image artifacts are still visible on these images as because the selected SI-PSF hypothesis does not compensate for all the space-variant aberrations. Using an SV-PSF hypothesis (Fig. 2.(c)) highly improves the quality and reproducibility of the reconstruction, as was already demonstrated in the simulations. Most of the

twin-image artifacts are removed from the reconstructions and the beads reconstructions are repeatable regardless of position in the FOV. It should be noted that the reconstruction quality obtained with the SV-PSF models depends on the PSF grid sampling ($P$), the number of modes ($Q$) used to model it, and the number of iterations of the IPR method. Indeed, $P$ and $Q$ should be chosen sufficiently high to capture the variation of the SV-PSF, in our case $P = 28, Q = 15$ and $N_{iter} = 200$. RMSE and SSIM values depending on $Q$ are provided in Fig.S5. Full FOV reconstructions are provided in Fig. S4.

## 5 CONCLUSION

In this paper, we propose a simple way to consider aberration and space-variant PSF in the commonly used IPR algorithm reconstruction. Using a calibrated PSF on a regular grid, the propagation model of the IPR algorithm is modified to account for SI-PSF or SV-PSF. When PSFs slowly vary in the FOV, modifying the propagation to account for SI-PSF is sufficient to obtain qualitatively better results than obtained with the classical IPR algorithm with a similar computational cost. However, the SV-PSF should be considered to perform more quantitative reconstructions and reduce the twin image artifacts. This approach relies on a modal decomposition of the SV-PSF. As the image formation model is more accurate, both RMSE and SSIM criteria show better values on simulated and experimental holograms. The computational time of the SV-PSF IPR algorithm iterations depends linearly on the number of modes used to capture the variation of the PSF. We have shown on simulated data, that this method leads to more reliable results than the classical IPR algorithm, even with fewer iterations. Both methods have been implemented with MATLAB and the source code is provided (see Sec. 6 and Supplementary File). The proposed aberration correction can be adapted for any other iterative reconstruction algorithms such as the Gerchberg-Saxton algorithm.

## 6 BACKMATTER

**Funding** The authors acknowledge for the grant ANR-23-CE51-0023 (project ATICS), and Swiss National Science Foundation research grant 200021_197107.

**Disclosures** The authors declare no conflicts of interest.

**Data Availability Statement** The source code and calibration data are available at https://github.com/braultd/Iterative-phase-retrieval-algorithm-for-shift-variant-PSF-in-optical-system-with-aberrations.

**Supplemental document** See Supplementary File for supporting content.

## REFERENCES

[1] James R Fienup. Reconstruction of an object from the modulus of its fourier transform. *Optics letters*, 3(1):27–29, 1978.

[2] Ralph W Gerchberg. A practical algorithm for the determination of plane from image and diffraction pictures. *Optik*, 35(2):237–246, 1972.
[3] Jianwei Miao, Pambos Charalambous, Janos Kirz, and David Sayre. Extending the methodology of x-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens. *Nature*, 400(6742):342–344, 1999.

[4] Tatiana Latychevskaia and Hans-Werner Fink. Solution to the twin image problem in holography. *Physical review letters*, 98(23):233901, 2007.

[5] Tatiana Latychevskaia, Conrad Escher, William Andregg, Michael Andregg, and Hans-Werner Fink. Direct visualization of charge transport in suspended (or free-standing) dna strands by low-energy electron microscopy. *Scientific reports*, 9(1):8889, 2019.

[6] Jean-Nicolas Longchamp, Stephan Rauschenbach, Sabine Abb, Conrad Escher, Tatiana Latychevskaia, Klaus Kern, and Hans-Werner Fink. Imaging proteins at the single-molecule level. *Proceedings of the National Academy of Sciences*, 114(7):1474–1479, 2017.

[7] Myung K Kim. Principles and techniques of digital holographic microscopy. *SPIE reviews*, 1(1):018005, 2010.

[8] Joseph W Goodman. *Introduction to Fourier optics*. Roberts and Company publishers, 2005.

[9] Dylan Brault, Thomas Olivier, Ferréol Soulez, Sachin Joshi, Nicolas Faure, and Corinne Fournier. Accurate unsupervised estimation of aberrations in digital holographic microscopy for improved quantitative reconstruction. *Optics Express*, 30(21):38383–38404, 2022.

[10] James G Nagy and Dianne P O'Leary. Restoring images degraded by spatially variant blur. *SIAM Journal on Scientific Computing*, 19(4):1063–1082, 1998.

[11] Ralf C Flicker and François J Rigaut. Anisoplanatic deconvolution of adaptive optics images. *JOSA A*, 22(3):504–513, 2005.

[12] Michael Hirsch, Suvrit Sra, Bernhard Schölkopf, and Stefan Harmeling. Efficient filter flow for space-variant multi-frame blind deconvolution. *2010 IEEE Computer Society Conference on Computer Vision and Pattern Recognition*, pages 607–614, 2010.

[13] Loïc Denis, Eric Thiébaut, Ferréol Soulez, Jean-Marie Becker, and Rahul Mourya. Fast approximations of shift-variant blur. *International Journal of Computer Vision*, 115:253–278, 2015.

[14] William S Cleveland. Robust locally weighted regression and smoothing scatterplots. *Journal of the American statistical association*, 74(368):829–836, 1979.

[15] Yibo Zhang, Hongda Wang, Yichen Wu, Miu Tamamitsu, and Aydogan Ozcan. Edge sparsity criterion for robust holographic autofocusing. *Optics letters*, 42(19):3824–3827, 2017.
# Iterative phase retrieval algorithm for space-variant PSF in optical systems with aberrations - Supplementary file

DYLAN BRAULT<fnref n="1" />,*, CORINNE FOURNIER<fnref n="1" />, AND TATIANA LATYCHEVSKAIA<fnref n="2" />,<fnref n="3" />

Compiled November 25, 2024

## SUPPLEMENTARY INFORMATION

### A. Supplementary tables

<table><caption>**Table S1.** Experimental parameters (*from the manufacturer, ThermoFisher Scientific, Inc.)</caption><tr><td>Wavelength</td><td>637.6 nm</td></tr><tr><td>Magnification</td><td>66.5</td></tr><tr><td>Numerical Aperture (NA)</td><td>1.2</td></tr><tr><td>Pixel pitch</td><td>82.5 nm</td></tr><tr><td>Total field of view</td><td>204×273 µm<fnref n="2" /></td></tr><tr><td>Typical defocus</td><td>17 µm</td></tr><tr><td>Beads diameter*</td><td>(1.0±0.06) µm</td></tr><tr><td>Beads refractive index*</td><td>1.587 (polystyrene)</td></tr><tr><td>Refractive index of immersion medium</td><td>1.47 (glycerol)</td></tr><tr><td>Coverslip thickness</td><td>(0.170±0.001) mm</td></tr></table>

<figure><img src="image_4.png" /><figcaption>**Fig. S1.** Principle of the IPR phase retrieval algorithm</figcaption></figure>

<table><caption>**Table S2.** Zernike coefficients used for the simulations of SI-PSF (all other coefficients are inferior to 0.1)</caption><tr><td>Zernike coefficients</td><td>α<sub>2</sub><sup>2</sup></td><td>α<sub>3</sub><sup>-3</sup></td><td>α<sub>3</sub><sup>-1</sup></td><td>α<sub>3</sub><sup>1</sup></td><td>α<sub>3</sub><sup>3</sup></td><td>α<sub>4</sub><sup>0</sup></td><td>α<sub>4</sub><sup>4</sup></td><td>α<sub>5</sub><sup>-1</sup></td><td>α<sub>6</sub><sup>-2</sup></td><td>α<sub>6</sub><sup>0</sup></td><td>α<sub>6</sub><sup>4</sup></td></tr><tr><td>Value for SI-PSF model</td><td>0.62</td><td>-0.25</td><td>4.91</td><td>-1.28</td><td>0.19</td><td>-4.71</td><td>0.45</td><td>-0.71</td><td>0.15</td><td>-0.8</td><td>0.20</td></tr></table>

### B. Supplementary figures

<figure><img src="image_5.png" /><figcaption>**Fig. S2.** Principle of the PSF interpolation model proposed in [11]</figcaption></figure>

<footnote n="1">¹Université Jean Monnet Saint-Etienne, CNRS, Institut d'Optique Graduate School, Laboratoire Hubert Curien UMR 5516, 42023 Saint-Etienne, France.</footnote>

<footnote n="2">²Center for Life Sciences, Paul Scherrer Institute, Forschungsstrasse 111, 5232 Villigen, Switzerland.</footnote>

<footnote n="3">³Department of Physics, University of Zurich, Winterthurerstrasse 190, 8057 Zurich, Switzerland.</footnote>

<footnote>*dylan.brault@telecom-st-etienne.fr</footnote>
<figure><img src="image_6.png" /><figcaption>**Fig. S3.** Left: Estimation of the phase of the Optical Transfer Function (OTF) for z = 0 depending on the position in the field of view on experimental data, $v_x$ and $v_y$ are the spatial frequencies in Fourier space. Right: Estimation of the defocus depending on the position in the field of view on experimental data. Black dots represent the regular grid on which the OTF and defocus have been estimated. The SV-PSF is computed as a shift-variant convolution between aberrated PSFs in focus and propagation kernels that correspond to the distances z (illustrated on the right side).</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>**Fig. S4.** Full field of view reconstructions of an experimental hologram with classical IPR algorithm, IPR algorithm with SI-PSF and IPR algorithm with SV-PSF. Phase of the reconstructions. Zooms on the reconstructions are provided in Fig.2.</figcaption></figure>
<figure><img src="image_8.png" /><figcaption>Fig. S5. Evolution of the error depending of the Q on simulated data. (P = 49)</figcaption></figure>

### C. Implementation remarks

Some implementation remarks are provided in this section to facilitate the implementation of the SV-PSF IPR algorithm.

**Remark #1:** To ensure that SVD does not introduce frequencies beyond $\frac{NA}{\lambda}$ (where NA is the numerical aperture of the microscope objective), each mode $\underline{m}_p$ is filtered in the Fourier domain such that frequencies beyond $\frac{NA}{\lambda}$ are set to 0.

**Remark #2:** As in holography PSFs are spatially extended, padding of the field of view is required to perform convolution. Since holography is not sensitive to phase offset (as only intensity is recorded on the sensor), it is common to choose the phase origin in the object plane such that one-padding can be performed. Thus, the phase of the (0,0) frequency of the propagator is set to 0. In this case:

$$ \underline{U}_{\text{diff}} = \underline{t} * \underline{h} = (1 + \underline{o}) * \underline{h} = 1 + \underline{o} * \underline{h} \quad (1) $$

where $\underline{o}$ is a complex-valued function that describes the perturbation caused by the sample. For SV-PSF propagation, the weights $\underline{w}_p$ have to be defined over the padded field of view. As it may be challenging to extrapolate them realistically, we suggest setting them to 0 out of the field of view and propagating only the perturbation $\underline{o}$ such that the pixel-wise multiplication between the weights and $\underline{t} - 1$ or $\underline{U}_{\text{diff}} - 1$ is expected to be 0 outside of the field of view. Finally, as the phase of the (0,0) frequency of the propagator should be set to 0, we set the phase of the (0,0) frequency of each mode $\underline{m}_p$ to 0 for the SV-PSF IPR algorithm.

**Remark #3:**

The modes $\underline{m}_p$ are not the PSF at p-th spatial location but an orthogonal set of space-invariant PSF modes, and the $\underline{w}_p$ are the associated weights maps that specify how the modes $\underline{m}_p$ should be

weighted together over the FOV. Indeed, the matrix $\underline{K}$ is composed of all the locally calibrated P PSF. As the modes are sorted according to their associated singular value, the low-rank approximation reduces the number of modes required. Indeed, the truncated SVD provides the minimal mean square error between $\underline{K}$ and the approximated matrix of rank $Q$ (Eckart–Young–Mirsky theorem).

### D. Supplementary algorithm

**Algorithm 1. Classical IPR algorithm**

1: **Input:** Data $\underline{d}$ (L × C image)

2: **Input:** Propagation kernel $\underline{h}$ (L × C image)

3: **Input:** Number of iterations $N_{\text{iter}}$ (integer)

4: **Output:** Transmission function $\underline{t}$ (L × C image)

5:

6: $\underline{U}_{\text{diff}} \leftarrow \sqrt{d} \quad \Leftrightarrow \text{Initialization of the wavefront on the sensor plane}$

7:

8: **for** $n_{\text{iter}} = 1$ **to** $N_{\text{iter}}$ **do**

9: $\underline{t} \leftarrow \underline{U}_{\text{diff}} * \underline{h}^{\dagger} \quad \Leftrightarrow \text{Back-propagation to the sample plane}$

10:

11: **for** $k = 1$ **to** L × C **do**

12: **if** $\underline{t}_k \notin \gamma$ **then** $\quad \Leftrightarrow \text{Applying constraints}$

13: $\underline{t}_k = 1$

14:

15: $\underline{U}_{\text{diff}} \leftarrow \underline{t} * \underline{h}$ $\quad \Leftrightarrow \text{Propagation to the sensor plane}$

16:

17: $\underline{U}_{\text{diff}} \leftarrow \underline{d} \odot \exp(i \times \arg(\underline{U}_{\text{diff}}))$ $\quad \Leftrightarrow \text{Updating modulus}$

18:

19: **return** $\underline{t}$
# MATLAB CODE

## A. SI-PSF Reconstruction

% a_SIPSF_Reconstruction.m

% File to reconstruct a hologram with a SIPSF hypothesis

% The PSF has been calibrated at the center of the field of view

clear all; close all;clc;

% Adding auxiliary functions to path
addpath(genpath('./Utils'))

% Load the hologram load('./Data/data.mat');

% Experimental parameters
pix nb x=1000; pix nb y=1000;
pix dx=82.5e-9;
pix dy=82.5e-9;
lambda=637.6e-9;
n 0=1.47;
ON=1.2;

% Pixel size of the padded images
pix_pad=pix;
pix_pad nb_x=2*pix nb_x;
pix_pad nb_y=2*pix nb_y;

% Defining functions of padding and cropping
cropimg=@(x) x(pix nb y/2+1:pix nb y/2+pix nb_x/2+1:pix nb_x/2+pix nb_x);
padimg=@(x) padarray(x,[pix nb_y/2 pix nb_x/2],1);

% Load the shift-invariant PSF h_SI estimated at the center of the field of view
disp('Loading the shift-invariant PSF...')
load('./Data/PSF_FOV_center.mat')

% Number of iterations of IPR algorithm
n_iter=200;

% Initialization of the wavefront in the hologram plane
disp('Initializing the reconstructed diffracted wavefront...')
U_diff=padimg sqrt(data);

for n=1:n_iter
disp strcat('Processing: ',' ',num2str(floor(100*n/n_iter)),' % ...'))

% Backpropagation of the wavefront into the sample plane
t=backpropagation(h SI,U_diff);

% Imposing constraints in the sample plane
mod t=abs(t);
arg_t=angle(t);
arg_t(arg_t<0)=0; % The phase of the sample is positive
mod t(:)=1; % The objects are transparent
t=mod_t.*exp(li*arg_t);

% Propagation of the wavefront into the hologram plane
U_diff=propagation(h SI,t);
% Imposing data values in the hologram plane
mod U_diff=padimg sqrt(data));
arg U_diff=angle(U_diff);
U_diff=mod U_diff.*exp(li*arg U_diff);

end
disp('Reconstruction process completed')

% Cropping the extended reconstruction
t=cropimg(t);

% Display the reconstructed transmission function
figure(1);
subplot(1,2,1);imshow(abs(t),[]);colormap(pink);colorbar;title('Modulus of the reconstructed transmittance')subplot(1,2,2);imshow(angle(t),[]);colormap(pink);colorbar;title('Phase of the reconstructed transmittance')

## B. SV-PSF Reconstruction

% b_SVPSF_Reconstruction.m

% File to reconstruct a hologram with a SVPSF hypothesis

% The SVPSF has been calibrated on regular grid

clear all; close all;clc;

% Adding auxiliary functions to path
addpath(genpath('./Utils'))

% Load the hologram

load('./Data/data.mat');

% Experimental parameters

pix nb_x=1000;
pix nb_y=1000;
pix dx=82.5e-9;
pix dy=82.5e-9;
lambda=637.6e-9;
n_0=1.47;
ON=1.2;

% Pixel size of the padded images

pix_pad=pix;
pix_pad nb_x=2*pix nb_x;
pix_pad nb_y=2*pix nb_y;

% Defining functions of padding and cropping

cropimg=@(x) x(pix nb_y/2+1:pix nb_y,pix nb_x/2+1:pix nb_x/2+pix nb_x);
padimg=@(x) padarray(x,[pix nb_y/2 pix nb_x/2],1);

% Load PSFs on regular grid of coordinates X_sub, Y_sub (Each column is a PSF)

disp('Loading the calibrated SV-PSF...');
load('./Data/SV Calibrated PSFs.mat');

% Coordinated of each pixels of the padded image

[X_field,Y_field]=meshgrid(-pix nb_x:pix nb_x-1,-pix nb_y:pix nb_y-1);

% Low rank approximation of the PSF

disp('Computing the low-rank approximation of the SV-PSF...');

Rank=15; % Rank of the approximation

[U,S,V]=SVPSF_LowRank(K,Rank);
[m,w]=SVPSF_ConvolutionKernel_WeightsComputation(U,S,V,X_sub,Y_sub,X_field,Y_field);
clear U;clear V;clear S; % Clearing the SVD variable to save memory

% Number of iterations of the IPR algorithm

n_iter=200;

% Initialization of the wavefront in the hologram plane

disp('Initializing the reconstructed diffracted wavefront...');
U_diff=padimg(Jsqrt(data));
% Shift-variant propagation and backpropagation functions
propagation_ =@(x) SVPSF_Propagation(m,w,mtfON,x);
backpropagation_ =@(x) SVPSF_BackPropagation(m,w,mtfON,x);

for n=1:n_iter
disp(strcat('Processing: ',' ','num2str тур(100*n/n_iter)),'% ...'))

% Backpropagation of the wavefront into the sample plane
t=backpropagation_(U_diff);

% Imposing constraints in the sample plane
mod_t=abs(t);
arg_t=angle(t);
arg_t(arg_t<0)=0; % The phase of the sample is positive
mod_t(:)=1; % The objects are transparent
t=mod_t.*exp(li*arg_t);

% Propagation of the wavefront into the hologram plane
U_diff=propagation_(t);

% Imposing data values in the hologram plane
mod_U_diff=padimg(sqrt(data));
arg_U_diff=angle(U_diff);
U_diff=mod_U_diff.*exp(li*arg_U_diff);

end

disp('Reconstruction process completed')

% Cropping the extended reconstruction
t=cropimg(t);

% Display the reconstructed transmission function
figure(1);
subplot(1,2,1);imshow(abs(t),[]);colormap(pink);colorbar;title('Modulus of the reconstructed
transmittance')subplot(1,2,2);imshow(angle(t),[]);colormap(pink);colorbar;title('Phase of the
reconstructed transmittance')

C. Auxiliary functions

function [Uprop]=propagation(propagation_kernel,x)

% This function compute the backpropagation of x using the propagation kernel propagation_kernel

% Inputs:

% propagation_kernel: Propagation kernel

% x: Wavefront to be propagated

% Outputs:

% Uprop: Backpropagated wavefront

Uprop=ifft2(fft2(propagation_kernel).*fft2(x));

end

function [Ubprop]=backpropagation(propagation_kernel,x)

% This function compute the backpropagation of x using the propagation kernel propagation_kernel

% Inputs:

% propagation_kernel: Propagation kernel

% x: Wavefront to be backpropagated

% Outputs:

% Ubprop: Backpropagated wavefront

Ubprop=ifft2(conj(fft2(propagation_kernel)).*fft2(x));

end
function [ConvolutionKernel,Weights] =
SVPSF_ConvolutionKernel_WeightsComputation(U,S,V,X_sub,Y_sub,X_field,Y_field)

% This function compute the mode and weights for SV-PSF

% Inputs :

% USV: Singular Value Decomposition of psfs

% X_sub, Y_sub: Coordinates of position of calibrated PSFs (coarse regular grid)

% X_field,Y_field: Coordinates of the pixel grid

% Outputs :

% ConvolutionKernel: Convolution modes for SV convolution

% Weights: Corresponding weights for SV convolution

Weights = [];

ConvolutionKernel = [];

Weights_vec=S*V; % Each column contains weights for one mode on a coarse regular grid
for i=1=size(Weights_vec,1)

% Reshaping the weights on the regular grid

Weights_sub=shape(Weights_vec(iallo,[size(X_sub,1),size(X_sub,2)]));

% Interpolation of the weights on the pixel grid

Weights (:,:i)=interp2(X_sub,Y_sub,Weights_sub,X_field,Y_field,'cubic',0);

% Reshaping the mode on the pixel grid

ConvolutionKernel (:,:i)=reshape(U(:,i),[size(X_field,1),size(X_field,2)]);

% Each column of U contains a mode

end

end

function [U S V] = SVPSF_LowRank(psfs,rank)

% This function compute the low-rank approximation of the SV-PSF

% Inputs :

% psfs : Calibrated PSFs (each column is a PSF - corresponds to matrix K)

% rank : Rank of the low rank approximation

% Outputs :

% USV: Singular Value Decomposition of psfs (Low Rank Approximation of rank rank)

[U S V] = svd(psfs,'econ');

V_ = V';

U = U(:,1:rank);

S = S(1:rank,1:rank);

V = V_(1:rank, :);

end

function [Uprop] = SVPSF_Propagation(ConvolutionKernel,Weights,mtfON,x)

% This function compute the shift-varying propagation

% Inputs :

% ConvolutionKernel : Convolution kernel obtained by SVD decomposition

% Weights : Weights map associated to the convolution kernel (obtained by SVD decomposition)

% mtfON : MTF correspond to filtering by numerical aperture

% x: Wavefront to be propagated (centered on 1)

% Outputs :

% Uprop: Propagated wavefront

% Suppression of the DC component of the signal

x = x-1;
% Initialization of the wavefront
Uprop=zeros(size(x));

% Low rank approximation of the convolution
for k=1=size(ConvolutionKernel,3)

% Computation of the Fourier transform of k-th mode
FTConvKern=fft2(ConvolutionKernel(:,:,k));
% Setting the phase at the origin of FTConvKern (thus, the phase origin is on the hologram plane)
FTConvKern(1,1)=abs(FTConvKern(1,1));

% Computation of the shift-variant convolution
Uprop=Uprop+ifft2(fft2(Weights(:,:,k).*x).*FTConvKern.*mtfON); end

% Adding DC component of the signal
Uprop=Uprop+1;
end

function [Uprop] = SVPSF_BackPropagation(ConvolutionKernel,Weights,mtfON,x)

% This function compute the shift-varying propagation

% Inputs:

% ConvolutionKernel: Convolution kernel obtained by SVD decomposition

% Weights: Weights map associated to the convolution kernel (obtained by SVD decomposition)

% mtfON: MTF corresponds to filtering by numerical aperture

% x: Wavefront to be propagated (centered on 1)

% Outputs:

% Uprop: Propagated wavefront

% Suppression of the DC component of the signal
x=x-1;

% Initialization of the wavefront
Uprop=zeros(size(x));

% Low rank approximation of the convolution
for k=1=size(ConvolutionKernel,3)

% Computation of the Fourier transform of k-th mode (conjugation for back-propagation)
FTConvKern=conj(fft2(ConvolutionKernel(:,:,k)));

% Setting the phase at the origin of FTConvKern (thus, the phase origin is on the sample plane)
FTConvKern(1,1)=abs(FTConvKern(1,1));

% Computation of the shift-variant convolution (the weights are conjugated for the backpropagation)
Uprop=Uprop+ifft2(fft2(conj(Weights(:,:,k)).*x).*FTConvKern.*mtfON);

end

% Adding DC component of the signal
Uprop=Uprop+1;
end