Linear mathematical model for the unique solution of 3D ptychographic iterative engine

Liqing Wu; Naijie Qi; Chengcheng Chang; Hua Tao; Xiaoliang He; Cheng Liu; Jianqiang Zhu

doi:10.3788/COL202422.050501

Abstract

Diffraction intensities of the 3D ptychographic iterative engine (3PIE) were written as a set of linear equations of the self-correlations of Fourier components of all sample slices, and an effective computing method was developed to solve these linear equations for the transmission functions of all sample slices analytically. With both theoretical analysis and numerical simulations, this study revealed the underlying physics and mathematics of 3PIE and demonstrated for the first time, to our knowledge, that 3PIE can generate mathematically unique reconstruction even with noisy data.

Keywords

3D imaging 3D ptychographic iterative engine linear equation group unique solution

1. Introduction

Recently, computational imaging has become a research hotspot in optical field, especially phase retrieval^[1–4]. Coherent diffraction imaging (CDI)^[5,6] is a kind of phase retrieval technique using various iterative algorithms. The G–S algorithm^[7,8] is the earliest CDI algorithm that records diffraction intensity at two separated planes. The error reduction (ER) algorithm^[9] and Fienup’s hybrid-input–output (HIO) algorithm^[10], which recorded only one frame of diffraction intensity, have much faster convergence and much better reconstruction quality than the G–S algorithm. Ptychography^[11,12] was invented by Walter Hoppe to reconstruct the objects with periodic structure and has been successfully used in material inspection with X-ray and high-energy electrons^[13,14]. By combining the CDI algorithm and ptychography technique together, Rodenburg proposed the ptychography iterative engine (PIE)^[15] to solve problems of twin image and low convergence in the classical CDI. PIE scans samples through a localized light beam to many positions and reconstructs the complex transmission function of the sample from diffractive intensities recorded at all scanning positions. The overlap between adjacent illuminating regions in PIE greatly improves its convergence speed and reconstruction quality, and PIE has been realized with visible light^[16], X-rays^[17,18], high-energy electron beams^[19], and terahertz waves^[20,21]. While the original PIE required exactly known illumination and sample positions, good reconstruction can be achieved by the extended ptychography iterative engine (ePIE) algorithm^[22] to reconstruct sample and illumination wavefront simultaneously and by annealing or cross-correlation algorithms to correct the scanning positions of the sample^[23,24], greatly improving the performance of PIE and extending its applications^[25–27]. Applying the multislice theory of electron microscopy^[28], 3D imaging can also be realized with PIE by regarding a 3D object as a series of 2D infinitely thin layers. Comparing to traditional 3D imaging methods such as optical coherent tomography^[29] and magnetic resonant tomography^[30], which generated intensity images, 3D ptychographic iterative engine (3PIE)^[31,32] can provide a high-quality 3D phase image for a transparent volume object rapidly. While 3PIE was first demonstrated experimentally with X-rays under geometric projection approximation^[33,34], it was also realized using visible light with diffraction taken into consideration^[32]. Single-shot 3PIE^[35,36] was also realized by recording subdiffraction patterns array with one detector exposure, making 3D phase imaging for dynamic imaging possible^[37]. 3PIE has shown good performance in 3D phase imaging; however, there is still no analytical theory to explain why it can work and to illustrate whether its reconstruction has mathematic uniqueness. In experiments, we were always not sure how reconstruction accuracy was affected by the optical alignment used, hindering the further development of 3PIE. Furthermore, since the analytical relationship between recorded diffraction intensities and reconstructed images has not been found by now, we cannot do quantitative and analytical error analysis on the reconstructed phase in experiments, making it impossible for 3PIE to be applied in fields of optical measurement and optical metrology, where the mathematical uniqueness of reconstruction and analytical error analysis are very critical^[38].

To investigate the underlying physics and mathematics of 3PIE algorithm, in this study diffraction intensities were written as a linear equation set of the self-correlations of Fourier components of sample slices, and the spatial components of all sample slices can be analytically determined by solving this linear equation set. Furthermore, an effective computing method that requires only small computer memory and can solve this linear equation set speedily was proposed. The influence of noise on this proposed linear model and computing method was also considered, demonstrating that both the linear model and a speedy computing method have strong noise immunization capability, and the influence of detector noise can be effectively reduced by simply dividing all recorded intensities into groups and adding each group together. In this paper, while theoretical analysis was illustrated, numerical simulations were also carried out to verify the feasibility of the proposed model and computing method. This study proves the mathematical uniqueness of the 3PIE algorithm for the first time and puts forward a speedy computing method to get analytical reconstruction from diffraction intensities, promoting the development of 3PIE in fields of optical measurement or metrology, where a strict unique mathematic solution and quantitative error analysis are required.

2. Theory and Method

2.1 Theoretical analysis

The optical alignment of 3PIE is schematically shown in Fig. 1, where the volume object is assumed to be composed of two layers, and the laser beam $P (x_{1}, y_{1})$ incident on the first layer is generated by a parallel beam passing through a tiny aperture. The interval between two object layers was assumed as $d_{1}$ , and their transmission functions were assumed as $S_{1} (x_{1}, y_{1})$ and $S_{2} (x_{2}, y_{2})$ , respectively. The distance of a CCD to the second object layer is $z$ . The light field leaving the first object layer $U_{1} (x_{1}, y_{1}) = P (x_{1}, y_{1}) S_{1} (x_{1}, y_{1})$ was regarded as the illumination of the second object layer after it propagated the distance of $d_{1}$ . The illumination arriving at the second object layer can be written as $U_{2} (x_{2}, y_{2}) = F^{- 1} {F (P (x_{1}, y_{1}) S_{1} (x_{1}, y_{1})) \cdot H (u, v)} = F^{- 1} {\tilde{P} (u, v) \otimes {\tilde{S}}_{1} (u, v) \cdot H (u, v)},$ (1)where $\tilde{P} (u, v)$ and ${\tilde{S}}_{1} (u, v)$ are Fourier transforms of $P (x_{1}, y_{1})$ and $S_{1} (x_{1}, y_{1})$ , respectively. $H (u, v) = \exp [i k d_{1} \cdot \sqrt{1 - {(λ u)}^{2} - {(λ v)}^{2}}]$ is the transfer function, where $λ$ means wavelength and $k$ is a wave vector. The symbol $F$ and $F^{- 1}$ are Fourier transform and inverse Fourier transform, respectively. $\otimes$ means a 2D convolution operator. The light field arriving at the detector was the Fresnel diffraction of the transmitted light of the second layer $U_{2} (x_{2}, y_{2}) S_{2} (x_{2}, y_{2})$ , and it can be written as $U_{3} (x, y) = \frac{1}{i λ z} \exp (i k z) \exp [i \frac{k}{2 z} (x^{2} + y^{2})] \iint_{- \infty}^{\infty} U_{2} (x_{2}, y_{2}) S_{2} (x_{2}, y_{2}) \cdot \exp [i \frac{k}{2 z} (x_{2}^{2} + y_{2}^{2})] \exp [- i \frac{2 π}{λ z} (x_{2} x + y_{2} y)] d x_{2} d y_{2} .$ (2)

Sign up for Chinese Optics Letters TOC. Get the latest issue of Chinese Optics Letters delivered right to you！Sign up now

Figure 1.Schematic diagram of 3PIE with two slices object.

By defining $S_{2}^{'} (x_{2}, y_{2}) = S_{2} (x_{2}, y_{2}) \exp [i \frac{k}{2 z} (x_{2}^{2} + y_{2}^{2})]$ , $f_{x} = \frac{x}{λ z}$ and $f_{y} = \frac{y}{λ z}$ , Eq. (2) can be rewritten as $U_{3} (λ z f_{x}, λ z f_{y}) = \frac{1}{i λ z} \exp (i k z) \exp [i λ z π (f_{x}^{2} + f_{y}^{2})] \cdot [{\tilde{U}}_{2} (f_{x}, f_{y}) \otimes {\tilde{S}}_{2}^{'} (f_{x}, f_{y})] .$ (3)

The intensity of the $(k, l)$ th pixel received by the detector can be written into discrete form as $I (k, l) = {\sum_{u_{1}, v_{1}} [\sum_{m_{1}, n_{1}} \tilde{P} (f_{m_{1}}, f_{n_{1}}) {\tilde{S}}_{1} (- f_{m_{1}} + f_{m_{2}}, - f_{n_{1}} + f_{n_{2}})] \cdot H (f_{u_{1}}, f_{v_{1}}) {\tilde{S}}_{2}^{'} (- f_{m_{2}} + f_{k}, - f_{n_{2}} + f_{l})} \cdot {\sum_{m_{2}, n_{2}} [\sum_{m_{2}, n_{2}} \tilde{P} (f_{m_{2}}, f_{n_{2}}) {\tilde{S}}_{1} (- f_{m_{1}^{'}} + f_{m_{2}^{'}}, - f_{n_{1}^{'}} + f_{n_{2}^{'}})] \cdot H (f_{u_{2}}, f_{v_{2}}) {\tilde{S}}_{2}^{'} (- f_{m_{2}^{'}} + f_{k}, - f_{n_{2}^{'}} + f_{l})} *,$ (4)where $f_{x_{1}} = \frac{x_{2}}{λ d_{1}}$ and $f_{y_{1}} = \frac{y_{2}}{λ d_{1}}$ , ${\tilde{S}}_{2}^{'} (f_{x}, f_{y})$ represents the Fourier transform of $S_{2}^{'} (x_{2}, y_{2})$ , and $*$ indicates conjugation. For simple discussion, Eq. (4) can be rewritten into a compact form as $I (k, l) = \sum_{\binom{m_{2}, n_{2}}{m_{2}^{'}, n_{2}^{'}}} \sum_{\binom{m_{1}, n_{1}}{m_{1}^{'}, n_{1}^{'}}} {\tilde{S}}_{1} (- f_{m_{1}} + f_{m_{2}}, - f_{n_{1}} + f_{n_{2}}) \cdot {\tilde{S}}_{1}^{*} (- f_{m_{1}^{'}} + f_{m_{2}^{'}}, - f_{n_{1}^{'}} + f_{n_{2}^{'}}) \cdot {\tilde{S}}_{2}^{'} (- f_{m_{2}} + f_{k}, - f_{n_{2}} + f_{l}) {\tilde{S}}_{2}^{' *} (- f_{m_{2}^{'}} + f_{k}, - f_{n_{2}^{'}} + f_{l}) \cdot H (f_{m_{2}}, f_{n_{2}}) H^{*} (f_{m_{2}^{'}}, f_{n_{2}^{'}}) \tilde{P} (f_{m 1}, f_{n 1}) {\tilde{P}}^{*} (f_{m_{1}^{'}}, f_{n_{1}^{'}}) .$ (5)

In 3PIE, the illumination $\tilde{P} (f_{m_{1}}, f_{n_{1}}) {\tilde{P}}^{*} (f_{m_{1}^{'}}, f_{n_{1}^{'}})$ and the transmission function $H (f_{m_{2}}, f_{n_{2}}) H^{*} (f_{m_{2}^{'}}, f_{n_{2}^{'}})$ are known, defined as a coefficient matrix $A = H (f_{m_{2}}, f_{n_{2}}) H^{*} (f_{m_{2}^{'}}, f_{n_{2}^{'}}) \cdot \tilde{P} (f_{m_{1}}, f_{n_{1}}) {\tilde{P}}^{*} (f_{m_{1}^{'}}, f_{n_{1}^{'}})$ . The unknowns are each layer of 3D object $x_{\binom{m_{1}, n_{1}, m_{2}, n_{2}}{m_{1}^{'}, n_{1}^{'}, m_{2}^{'}, n_{2}^{'}}} = \sum_{\binom{m_{2}, n_{2}}{m_{2}^{'}, n_{2}^{'}}} \sum_{\binom{m_{1}, n_{1}}{m_{1}^{'}, n_{1}^{'}}} {\tilde{S}}_{1} (- f_{m_{1}} + f_{m_{2}}, - f_{n_{1}} + f_{n_{2}}) \cdot {\tilde{S}}_{1}^{*} (- f_{m_{1}^{'}} + f_{m_{2}^{'}}, - f_{n_{1}^{'}} + f_{n_{2}^{'}}) {\tilde{S}}_{2}^{'} (- f_{m_{2}} + f_{k}, - f_{n_{2}} + f_{l}) {\tilde{S}}_{2}^{' *} (- f_{m_{2}^{'}} + f_{k}, - n_{2}^{'} + f_{l})$ , and the intensity matrix of detector is defined as $B = I (k, l) |_{(k = 1, \dots, K; l = 1, \dots, L)}$ . All linear equations of Eq. (5) can be written in the matrix form: $A X = B$ , and the solution of Eq. (5) is written as $X = A^{- 1} B .$ (6)

Assuming that the detector records the effective intensity information of $M \times N$ pixels, then there are $M \times N$ linear equations in Eq. (6). Undoubtedly, as long as the number of equations is greater than the number of unknowns, all the $x_{\binom{m_{1}, n_{1}, m_{2}, n_{2}}{m_{1}^{'}, n_{1}^{'}, m_{2}^{'}, n_{2}^{'}}}$ can be calculated. The object information cannot be obtained directly from these computed $x_{\binom{m_{1}, n_{1}, m_{2}, n_{2}}{m_{1}^{'}, n_{1}^{'}, m_{2}^{'}, n_{2}^{'}}}$ , which always include all spectral components of each layer. Then, from all $x_{\binom{m_{1}, n_{1}, m_{2}, n_{2}}{m_{1}^{'}, n_{1}^{'}, m_{2}^{'}, n_{2}^{'}}}$ we can choose specific pixel of $m_{2} = m_{20}$ , $n_{2} = n_{20}$ , $m_{1}^{'} = m_{10}^{'}$ , $n_{1}^{'} = n_{10}^{'}$ , $m_{2}^{'} = m_{20}^{'}$ and $n_{2}^{'} = n_{20}^{'}$ and pick a new vector $x_{1}^{'}$ about $m_{1}, n_{1}$ as variables, $x_{1}^{'} = {\tilde{S}}_{1} (- f_{m_{1}} + f_{m_{2}}, - f_{n_{1}} + f_{n_{2}}) [{\tilde{S}}_{1}^{*} (- f_{m_{1}^{'}} + f_{m_{2}^{'}}, - f_{n_{1}^{'}} + f_{n_{2}^{'}}) \cdot {\tilde{S}}_{2}^{'} (- f_{m_{2}} + f_{k}, - f_{n_{2}} + f_{l}) {\tilde{S}}_{2}^{' *} (- f_{m_{2}^{'}} + f_{k}, - n_{2}^{'} + f_{l})] = C_{\binom{m_{20}, n_{20}, m_{10}^{'}}{n_{10}^{'}, m_{20}^{'}, n_{20}^{'}}} {\tilde{S}}_{1} (- f_{m_{1}}, - f_{n_{1}}) |_{\binom{m_{1} = 1, \dots, M^{'}}{n_{1} = 1, \dots, N^{'}}},$ (7)where $C_{\binom{m_{20}, n_{20}, m_{10}^{'}}{n_{10}^{'}, m_{20}^{'}, n_{20}^{'}}}$ is a constant with value determined by $m_{20}, n_{20}, m_{10}^{'}, n_{10}^{'}, m_{20}^{'}, n_{20}^{'}$ . Physically, the light field multiplied by a constant is essentially the same as the original light field. Therefore, $x_{1}^{'}$ is equal to ${\tilde{S}}_{1} (- f_{m_{1}}, - f_{n_{1}}) |_{\binom{m_{1} = 1, \dots, M^{'}}{n_{1} = 1, \dots, N^{'}}}$ and the first layer $S_{1} (x_{1}, y_{1})$ can be obtained by doing inverse Fourier transform on $x_{1}^{'}$ . Similarly, the second layer $S_{2} (x_{2}, y_{2})$ can be obtained by another vector as $x_{2}^{'} = C_{\binom{m_{20}^{'}, n_{20}^{'}, m_{10}^{'}}{n_{10}^{'}, m_{20}^{'}, n_{20}^{'}}} {\tilde{S}}_{2} (- f_{m_{2}}, - f_{n_{2}}) |_{\binom{m_{2} = 1, \dots, M^{'}}{n_{2} = 1, \dots, N^{'}}} .$ (8)

Since the number of unknown elements $X_{\binom{m_{1}, n_{1}, m_{2}, n_{2}}{m_{1}^{'}, n_{1}^{'}, m_{2}^{'}, n_{2}^{'}}}$ is very huge and is much larger than the pixel number of $M \times N$ of detector in most cases, Eq. (6) cannot be solved with one frame of recorded diffraction intensity. The condition for getting a unique reconstruction is that the coefficient matrix A is of full rank. When the sample is shifted by a distance $m_{x} δ_{x}$ and $n_{y} δ_{y}$ along the $x$ and $y$ directions, respectively, the phase-shifting factor $e^{i [(- f_{m_{1}} + f_{m_{2}}) m_{x} δ_{x} + (- f_{n_{1}} + f_{n_{2}}) n_{y} δ_{y}]}$ as a known term should be multiplied to A, and more linear equations are obtained. It is easy to get A of full rank when the sample is shifted to positions with random intervals. By scanning the sample to many positions, we can get a huge linear equation group in Eq. (5), and we can compute all $x_{\binom{m_{1}, n_{1}, u_{1}, v_{1}}{m_{2}, n_{2}, u_{2}, v_{2}}}$ and corresponding $S_{1} (x_{1}, y_{1})$ and $S_{2} (x_{2}, y_{2})$ analytically.

The above analysis is on a volume object composed of only two layers, but similar analysis can also be carried out on volume object composed of $L$ layers, and the only difference lies in that the above unknown element $X_{\binom{m_{1}, n_{1}, m_{2}, n_{2}}{m_{1}^{'}, n_{1}^{'}, m_{2}^{'}, n_{2}^{'}}}$ will be replaced by $\sum_{\binom{m_{M, \dots}, n_{N, \dots}, m_{2}, n_{2} m_{1}, n_{1}}{m_{M}^{'}, \dots, n_{N}^{'}, \dots, m_{2}^{'}, n_{2}^{'} m_{1}^{'}, n_{1}^{'}}} {\tilde{S}}_{1} (- f_{m_{1}} + m_{2}, - f_{n_{1}} + n_{2}) \cdot {\tilde{S}}_{1}^{*} (- f_{m_{1}^{'}} + m_{2}^{'}, - f_{n_{1}^{'}} + n_{2}^{'}) {\tilde{S}}_{2}^{'} (- m_{2} + m, - n_{2} + l) \cdot {\tilde{S}}_{2}^{' *} (- m_{2}^{'} + k, - n_{2}^{'} + l) \dots .$ The mathematical analysis is the same as that shown in Eq. (1) to Eq. (8). If the volume object is sliced into $L$ layers with size of $M^{'} \times N^{'}$ , there will be ${(M^{'} N^{'})}^{2 L}$ unknowns $x_{\binom{m_{M, \dots}, n_{N, \dots}, m_{2}, n_{2} m_{1}, n_{1}}{m_{M}^{'}, \dots, n_{N}^{'}, \dots, m_{2}^{'}, n_{2}^{'} m_{1}^{'}, n_{1}^{'}}}$ to be solved. The largest number of uncorrelated linear equations available from one frame diffractive patterns is $M N$ ; then, to get ${(M^{'} N^{'})}^{2 L}$ uncorrelated linear equations, the sample should be scanned to at least $\frac{{(M^{'} N^{'})}^{2 L}}{M N}$ positions. It is obvious that with the increasing object layer number $L$ , the required scanning positions exponentially increase. For experiments, where the number of scanning sample positions cannot be very huge because positioning error always accelerates with scanning range, a good reconstruction can be achieved with 3PIE when the sample is always sliced into a very limited number of layers or when the sample has a very limited number of spatial components; that is, $L$ or $M^{'} N^{'}$ always takes small values.

2.2. Efficient computing method

In the above mathematical analysis, to compute all ${(M^{'} N^{'})}^{2 L}$ unknown terms we need ${(M^{'} N^{'})}^{2 L}$ uncorrelated linear equations; then the size of A is ${(M^{'} N^{'})}^{2 L} \times {(M^{'} N^{'})}^{2 L}$ , which is an unreasonably huge number for most computer stations. Thus, it is impossible to compute all unknown terms by directly using Eq. (6).

We can find from the above analysis that only a very small number of computed unknown terms were finally applied to reconstruct the sample slices, and it is not essential to compute all of them. Furthermore, many $x_{\binom{m_{M, \dots}, n_{N, \dots}, m_{2}, n_{2} m_{1}, n_{1}}{m_{M}^{'}, \dots, n_{N}^{'}, \dots, m_{2}^{'}, n_{2}^{'} m_{1}^{'}, n_{1}^{'}}}$ have zero values, and these terms need not be computed. Figures 2(a) and 2(b) show the amplitude transmissions of two sample layers, and Figs. 2(c) and 2(d) show the modulus of their Fourier spectrum ${\tilde{S}}_{1} (f_{m_{1}}, f_{n_{1}})$ and ${\tilde{S}}_{2} (f_{m_{2}}, f_{n_{2}})$ in log scale, respectively. Figure 2(e) shows ${\tilde{S}}_{1} (- f_{m_{1}} + f_{m_{2}}, - f_{n_{1}} + f_{n_{2}}) {\tilde{S}}_{1}^{*} (- f_{m_{1}^{'}} + f_{m_{2}^{'}}, - f_{n_{1}^{'}} + f_{n_{2}^{'}}) {\tilde{S}}_{2}^{'} (- f_{m_{2}} + f_{k}, - f_{n_{2}}) {\tilde{S}}_{2}^{' *} (- f_{m_{2}^{'}} + f_{k}, - f_{n_{2}^{'}})$ , where we can clearly find that most of values are very close to zeros, except pixels around the center. Thus, it is possible to find an efficient computing method without huge computer memory to calculate the sample’s transmission function.

Figure 2.Amplitude transmissions and spectra of two layers.

To illustrate the generation of diffraction intensity with Eq. (5) intuitively, a $3 \times 3$ $\tilde{P} (f_{m}, f_{n})$ , a $3 \times 3$ ${\tilde{S}}_{1} (f_{m}, f_{n})$ , and a $3 \times 3$ ${\tilde{S}}_{2} (f_{m}, f_{n})$ are used in Fig. 3 to show the computation of ${\tilde{U}}_{2} (- 2, - 2)$ , ${\tilde{U}}_{2} (- 2, - 1)$ , ${\tilde{U}}_{2} (- 1, - 2)$ , and ${\tilde{U}}_{2} (- 1, - 1)$ with Figs. 3(a)–3(d), respectively, where orange grids $\bar{{\tilde{S}}_{1}} (f_{m}, f_{n})$ indicating the reversed spatial component matrix of the first layer were shifted by varying unities in the $x$ and $y$ directions with respect to the green grids indicating the illumination spatial components $\tilde{P} (f_{m}, f_{n})$ . ${\tilde{U}}_{2} (- 2, - 2)$ , ${\tilde{U}}_{2} (- 2, - 1)$ , ${\tilde{U}}_{2} (- 1, - 2)$ , and ${\tilde{U}}_{2} (- 1, - 1)$ can be written as Eq. (9), ${\begin{matrix} {\tilde{U}}_{2} (- 2, - 2) = \bar{{\tilde{S}}_{1}} (- 1, - 1) \tilde{P} (- 1, - 1) H (- 2, - 2) \\ {\tilde{U}}_{2} (- 2, - 1) = [\bar{{\tilde{S}}_{1}} (- 1, 0) \tilde{P} (- 1, - 1) + \bar{{\tilde{S}}_{1}} (- 1, - 1) \tilde{P} (- 1, 0)] H (- 2, - 1) \\ {\tilde{U}}_{2} (- 1, - 2) = [\bar{{\tilde{S}}_{1}} (0, - 1) \tilde{P} (- 1, - 1) + \bar{{\tilde{S}}_{1}} (- 1, - 1) \tilde{P} (0, - 1)] H (- 1, - 2) \\ {\tilde{U}}_{2} (- 1, - 1) = [\bar{{\tilde{S}}_{1}} (0, 0) \tilde{P} (- 1, - 1) + \bar{{\tilde{S}}_{1}} (0, - 1) \tilde{P} (- 1, 0) + \bar{{\tilde{S}}_{1}} (- 1, 0) \tilde{P} (0, - 1) + \bar{{\tilde{S}}_{1}} (- 1, - 1) \tilde{P} (0, 0)] H (- 1, - 1) \end{matrix} .$ (9)

Figure 3.Formation of ${\tilde{U}}_{2} (f_{m}, f_{n})$ on varying pixels.

Similarly, Figs. 4(a)–4(d) illustrate the formation of ${\tilde{U}}_{3} (- 3, - 3)$ , ${\tilde{U}}_{3} (- 3, - 2)$ , ${\tilde{U}}_{3} (- 2, - 3)$ , and ${\tilde{U}}_{3} (- 2, - 2)$ , respectively, where $\bar{{\tilde{S}}_{2}} (f_{m}, f_{n})$ in light pink indicating the reversed spatial component matrix of the second layer was shifted by varying unities in the $x$ and $y$ directions with respect to the blue grids indicating the illumination spatial components ${\tilde{U}}_{2} (f_{m}, f_{n})$ ; they can be written as ${\begin{matrix} {\tilde{U}}_{3} (- 3, - 3) = \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 2, - 2) \\ {\tilde{U}}_{3} (- 3, - 2) = \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 2, - 1) + \bar{{\tilde{S}}_{2}} (- 1, 0) {\tilde{U}}_{2} (- 2, - 2) \\ {\tilde{U}}_{3} (- 2, - 3) = \bar{{\tilde{S}}_{2}} (0, - 1) {\tilde{U}}_{2} (- 2, - 2) + \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 1, - 2) \\ {\tilde{U}}_{3} (- 2, - 2) = \bar{{\tilde{S}}_{2}} (0, 0) {\tilde{U}}_{2} (- 2, - 2) + \bar{{\tilde{S}}_{2}} (0, - 1) {\tilde{U}}_{2} (- 2, - 1) + \bar{{\tilde{S}}_{2}} (- 1, 0) {\tilde{U}}_{2} (- 1, - 2) + \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 1, - 1) \end{matrix},$ (10) ${\begin{matrix} I (- 3, - 3) = {| \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 2, - 2) |}^{2} \\ I (- 3, - 2) = {| \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 2, - 1) |}^{2} + {| \bar{{\tilde{S}}_{2}} (- 1, 0) {\tilde{U}}_{2} (- 2, - 2) |}^{2} + {\tilde{U}}_{2} (- 2, - 1) \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2}^{*} (- 2, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, 0) \\ + {\tilde{U}}_{2} (- 2, - 2) \bar{{\tilde{S}}_{2}} (- 1, 0) {\tilde{U}}_{2}^{*} (- 2, - 1) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, - 1) \\ I (- 2, - 3) = {| \bar{{\tilde{S}}_{2}} (0, - 1) {\tilde{U}}_{2} (- 2, - 2) |}^{2} + {| \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 1, - 2) |}^{2} + {\tilde{U}}_{2} (- 2, - 2) \bar{{\tilde{S}}_{2}} (0, - 1) {\tilde{U}}_{2}^{*} (- 1, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, - 1) \\ + {\tilde{U}}_{2} (- 1, - 2) \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2}^{*} (- 2, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (0, - 1) \\ I (- 2, - 2) = {| \bar{{\tilde{S}}_{2}} (0, 0) {\tilde{U}}_{2} (- 2, - 2) |}^{2} + {| \bar{{\tilde{S}}_{2}} (0, - 1) {\tilde{U}}_{2} (- 2, - 1) |}^{2} + {| \bar{{\tilde{S}}_{2}} (- 1, 0) {\tilde{U}}_{2} (- 1, - 2) |}^{2} + {| \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 1, - 1) |}^{2} \\ + \bar{{\tilde{S}}_{2}} (0, 0) {\tilde{U}}_{2} (- 2, - 2) {\tilde{U}}_{2}^{*} (- 2, - 1) {\bar{{\tilde{S}}_{2}}}^{*} (0, - 1) + \bar{{\tilde{S}}_{2}} (0, 0) {\tilde{U}}_{2} (- 2, - 2) {\tilde{U}}_{2}^{*} (- 1, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, 0) \\ + \bar{{\tilde{S}}_{2}} (0, 0) {\tilde{U}}_{2} (- 2, - 2) {\tilde{U}}_{2}^{*} (- 1, - 1) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, - 1) + \bar{{\tilde{S}}_{2}} (0, - 1) {\tilde{U}}_{2} (- 2, - 1) {\tilde{U}}_{2}^{*} (- 2, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (0, 0) \\ + \bar{{\tilde{S}}_{2}} (0, - 1) {\tilde{U}}_{2} (- 2, - 1) {\tilde{U}}_{2}^{*} (- 1, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, 0) + \bar{{\tilde{S}}_{2}} (0, - 1) {\tilde{U}}_{2} (- 2, - 1) {\tilde{U}}_{2}^{*} (- 1, - 1) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, - 1) \\ + \bar{{\tilde{S}}_{2}} (- 1, 0) {\tilde{U}}_{2} (- 1, - 2) {\tilde{U}}_{2}^{*} (- 2, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (0, 0) + \bar{{\tilde{S}}_{2}} (- 1, 0) {\tilde{U}}_{2} (- 1, - 2) {\tilde{U}}_{2}^{*} (- 2, - 1) {\bar{{\tilde{S}}_{2}}}^{*} (0, - 1) \\ + \bar{{\tilde{S}}_{2}} (- 1, 0) {\tilde{U}}_{2} (- 1, - 2) {\tilde{U}}_{2}^{*} (- 1, - 1) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, - 1) + \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 1, - 1) {\tilde{U}}_{2}^{*} (- 2, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (0, 0) \\ + \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 1, - 1) {\tilde{U}}_{2}^{*} (- 2, - 1) {\bar{{\tilde{S}}_{2}}}^{*} (0, - 1) + \bar{{\tilde{S}}_{2}} (- 1, - 1) {\tilde{U}}_{2} (- 1, - 1) {\tilde{U}}_{2}^{*} (- 1, - 2) {\bar{{\tilde{S}}_{2}}}^{*} (- 1, 0) \end{matrix} .$ (11)

Figure 4.Formation of ${\tilde{U}}_{3} (f_{m}, f_{n})$ on varying pixels.

When multiplied by a constant in spatial domain or spatial frequency domain, the sample’s transmission function does not change essentially, and for simplicity we can assume $\bar{{\tilde{S}}_{2}} (- 1, - 1)$ has a value of 1.0 or other given values without losing generality. Then ${\tilde{U}}_{2} (- 2, - 2)$ can be computed as ${\tilde{U}}_{2} (- 2, - 2) = \sqrt{I_{3} (- 3, - 3)}$ with the first equation of Eq. (12), and then $\bar{{\tilde{S}}_{1}} (- 1, - 1)$ can be determined as $\bar{{\tilde{S}}_{1}} (- 1, - 1) = \frac{\sqrt{I_{3} (- 3, - 3)}}{\tilde{P} (- 1, - 1) H (- 2, - 2)}$ with the first equation of Eq. (9). According to Figs. 3(b) and 4(b), $\bar{{\tilde{S}}_{2}} (- 1, 0)$ and $\bar{{\tilde{S}}_{1}} (- 1, 0)$ can be computed using the computed $\bar{{\tilde{S}}_{2}} (- 1, - 1)$ and $\bar{{\tilde{S}}_{1}} (- 1, - 1)$ . For clarity, we defined $\bar{{\tilde{S}}_{1}} (- 1, 0) = x$ , $\bar{{\tilde{S}}_{2}} (- 1, 0) = y$ , and the intensity $I^{t} (- 3, - 2)$ at the $(- 3, - 2)$ th pixel when the sample was scanned to the $t$ th position can be written as Eq. (12), ${\begin{matrix} I^{t} (- 3, - 2) = {| A_{1} y |}^{2} + {| A_{2} x |}^{2} + {| A_{3} |}^{2} + A_{1} A_{2}^{*} y x^{*} + A_{1} A_{3}^{*} y + A_{2} A_{1}^{*} x y^{*} + A_{2} A_{3}^{*} x + A_{3} A_{1}^{*} y^{*} + A_{3} A_{2}^{*} x^{*} \\ A_{1} = \bar{{\tilde{S}}_{1}} (- 1, - 1) e^{i Δ_{t}} \tilde{P} (- 1, - 1) H (- 2, - 2) e^{i ε_{t}} \\ A_{2} = \bar{{\tilde{S}}_{2}} (- 1, - 1) e^{i θ_{t}} \tilde{P} (- 1, - 1) H (- 2, - 1) e^{i δ_{t}} \\ A_{3} = \bar{{\tilde{S}}_{2}} (- 1, - 1) e^{i θ_{t}} \bar{{\tilde{S}}_{1}} (- 1, - 1) e^{i Δ_{t}} \tilde{P} (- 1, 0) H (- 2, - 1) \end{matrix},$ (12)where $e^{i Δ_{t}}$ , $e^{i ε_{t}}$ , $e^{i θ_{t}}$ and $e^{i δ_{t}}$ are additional phase factors of $\bar{{\tilde{S}}_{1}} (- 1, - 1)$ , $\bar{{\tilde{S}}_{2}} (- 1, 0)$ , $\bar{{\tilde{S}}_{2}} (- 1, - 1)$ , and $\bar{{\tilde{S}}_{1}} (- 1, 0)$ caused by the $t$ th shifting of the sample, respectively. Since ${| A_{1} y |}^{2}, {| A_{2} x |}^{2}, {| A_{3} |}^{2}$ do not change with the sample’s positions, these three terms can be eliminated by subtracting $I^{1} (- 3, - 2)$ from $I^{t_{1}} (- 3, - 2) |_{t_{1} = 2, 3, 4, 5, 6, 7}$ to yield six linear equations. Then, with these six linear equations, we can compute six unknowns: ${y x^{*}, y, x y^{*}, x, y^{*}, x^{*}}$ of the Eq. (12); then the values of $\bar{{\tilde{S}}_{2}} (- 1, 0)$ and $\bar{{\tilde{S}}_{1}} (- 1, 0)$ are computed. With the same strategy, the values of $\bar{{\tilde{S}}_{1}} (0, 0)$ and $\bar{{\tilde{S}}_{2}} (0, 0)$ can be computed in the next step, and all other spatial components of two sample layers can be computed in the same way. The transmission functions of the two sample layers can be computed by doing inverse Fourier transform on all computed ${\tilde{S}}_{1} (f_{m_{1}}, f_{n_{1}})$ and ${\tilde{S}}_{2} (f_{m_{2}}, f_{n_{2}})$ . To be suitable for large matrix objects, the point-by-point calculation takes a certain amount of time. However, when calculating objects with large layers, compared with the traditional 3PIE algorithm, which needs to wait for convergence, this method can compute each layer at the same time without extra time cost.

With the above computing method, two spatial components of sample slices can be computed in each step using seven linear equations, and the computer memory required was very small. Then solving Eq. (7) becomes quite easier than directly using Eq. (6). A two-layer sample was used as an example in the above analysis, and the transmission function of the volume object composed of many layers can also be computed in a similar way. When the sample was composed of three layers, four layers, and $L$ layers, the number of diffraction patterns required will be 13, 21, and $L^{2} + L + 1$ , respectively.

3. Numerical Simulations

To check the feasibility of the above theoretical analysis and proposed computing method, a series of numerical simulations were carried out. Two biological images of $512 \times 512$ pixels shown in Figs. 5(a) and 5(b) were used as amplitude transmissions of two layers of a volume sample. Two images shown in Figs. 5(c) and 5(d) were used as phase retardations of two layers, respectively. The interval between two layers was assumed as 1 mm. The probe light illuminating on the sample was a parallel laser beam of 632.8 nm passing through an aperture 0.7 mm in radius, and the distance of this aperture from the sample was 30 mm, equal to the distance from the sample to the detector. The amplitude and phase of illumination are shown in Figs. 5(e) and 5(f), respectively. The strength of the Fourier components of two sample layers and the illumination are shown in Figs. 5(g), 5(h) and 5(i) in log scale, respectively. When the sample was shifted by distances of (450 µm, 450 µm), (446 µm, 900 µm), (450 µm, 1320 µm), (890 µm, 450 µm), (905 µm, 920 µm), (900 m, 1361 µm), and (1330 µm, 450 µm), seven frames of diffraction patterns shown in Fig. 6 can be obtained. The pixel size of the detector was assumed to be 9 µm.

Figure 5.Object and illumination. (a) and (b) are the amplitude of two layers of object; (c) and (d) are the corresponding phase of two layers of object; (e) and (f) are the amplitude and phase of illumination; the spectra in log scale of two layers object and illumination are shown in (g)–(i), respectively. The scale bar of (a) is suitable for (b)–(d) and (g)–(i).

$Seven diffraction patterns used in computation.$

Figure 6.Seven diffraction patterns used in computation.

With the computing method discussed above, the strengths of the Fourier components of two sample layers computed with Eqs. (9 )–(12) are shown in Figs. 7(a) and 7(b) in log scale, respectively. For quantitative comparison, the differences between the reconstructed images and the original images were calculated based on the formula: $error = | | O_{re} | - | O_{origin} | |$ , shown in Figs. 7(c) and 7(d). We can find that the difference is around 0.05%, which is the computing accuracy of a common desktop of 32 bits.

Figure 7.Reconstructed spectra of two layers. (a) and (b) represent the recovered spectra of two layers, and (c) and (d) depict the difference between reconstructed spectrum and original spectrum.

By doing inverse Fourier transform on computed Fourier components in Fig. 7, we get the modulus and phase of two sample layers, shown in Figs. 8(a)–8(d). The differences of computed modulus and phase to their original values are shown in Figs. 8(e)–8(h), which are all at the scale of about $10^{- 2}$ , and are the computing accuracy of a common workstation of 32 bits. Results in Figs. 7 and 8 perfectly match our theoretical expectations of Eqs. (10 )–(12), proving the correctness of the above theoretical analysis and suggested computing methods.

Figure 8.Reconstruction of two layers object. (a)–(d) are the amplitudes and phases of two layers object; (e)–(h) are the differences of modulus and phases to their original values.

In the above studies, we did not touch experimental factors, including detector noises, which are a kind of inevitable error source of PIE experiments. If there is random noise $Δ B$ in the diffraction intensity, the linear equation set will become $A X = B + Δ B$ . Then, spatial components of sample slices cannot be accurately computed by directly using noisy diffraction intensities. If the sample was scanned at many positions to record a large enough number of diffraction intensities, we can add a large number of linear equations corresponding to the same detector pixel together as $A_{1} X + A_{2} X + \dots + A_{N} X = B_{1} + B_{2} + \dots + B_{N} + Δ B_{1} + Δ B_{2} + \dots + Δ B_{N} .$ (13)

When $N$ is large enough, $\sum_{1}^{N} Δ B_{m}$ will become close to zero, and Eq. (13) will become $\sum_{1}^{N} A_{m} X = \sum_{1}^{N} B_{m}$ . Then, $X$ , without the influence of detector noise, can also be computed as $X = {(\sum_{1}^{N} A_{m})}^{- 1} \sum_{1}^{N} B_{m}$ . That means that, by shifting sample to more positions and recording more diffraction patterns, we can remarkably suppress the influence of detector noise and get accurate reconstruction for 3PIE with the above illustrated analytical method. As the approximation holds only when $Δ B$ takes a small value, this method cannot be available when external noises are too large.

To verify the robustness of this anti-noise computing method, another set of simulations was carried out by adding Poisson noise to diffractive intensities, shown in Fig. 6. The sample was shifted by $7 \times 7$ positions, yielding 49 frames of diffraction patterns, and Poisson noise [Fig. 9(a)] with strength between $- 1$ and 2 was added to each pattern, resulting in the 20 dB signal-to-noise ratio (SNR). These noisy intensities are divided into seven groups; after diffraction patterns in each group are summed up, seven frames of new hybrid diffractive intensities, shown in Figs. 9(b)–9(h), are obtained.

$Noise and seven new diffraction patterns. (a) is the Poisson noise. (b)–(h) are seven new hybrid diffraction patterns.$

Figure 9.Noise and seven new diffraction patterns. (a) is the Poisson noise. (b)–(h) are seven new hybrid diffraction patterns.

With our suggested computing method, the spatial components of two sample layers were computed from hybrid diffraction patterns in Fig. (10), where Figs. 10(a) and 10(b) are the modulus of computed spectral components of the two sample layers in log scale. By doing inverse Fourier transform, the complex transmission function of the two sample layers can be obtained. Figures 10(c) and 10(d) are the modulus of the two layers, and Figs. 10(e) and 10(f) are corresponding phases. Figures 10(g) and 10(h) show the differences between the recovered modulus and their corresponding original values, respectively. The maximum difference is about 0.5%, which matches our expectation well that the influence of noise can be effectively eliminated by shifting the sample to more positions and recording more diffraction intensities.

Figure 10.Results reconstructed from noisy data. (a) and (b) are the recovered spectra in log scale; (c) and (d) are the modulus of the reconstructed two-layer object; the corresponding phase is shown in (e) and (f); amplitude differences to original image are shown in (g) and (h).

4. Conclusions

The underlying physical mechanism and mathematics of the 3PIE imaging method was revealed by writing its diffraction intensities as a linear equation set. The spatial components of all sample slices can be analytically determined using an efficient computing method to solve this linear equation set. The robustness of this suggested computing method in dealing with noisy data was also studied, and it was demonstrated that the influence of detector noise can be effectively eliminated by simply dividing many recorded intensities into groups and summing each group up. While theoretical analysis was illustrated, numerical simulations were also carried out to verify the feasibility of the proposed model and computing method by taking a two-slice thick sample as an example. This study clarified the mathematical uniqueness of the 3PIE algorithm for the first time and suggested a speedy computing method to get analytical reconstruction from diffraction intensity, breaking the theoretical bottleneck that hinders the application of 3PIE in fields of optical measurement or metrology, where mathematical uniqueness and error analysis are very crucial.

References

[1] J. X. Li, T. Y. Gan, B. J. Bai et al. Massively parallel universal linear transformations using a wavelength-multiplexed diffractive optical network. Adv. Photonics, 5, 016003(2023).

[2] L. P. Lu, J. J. Li, Y. F. Shu et al. Hybrid brightfield and darkfield transport of intensity approach for high-throughput quantitative phase microscopy. Adv. Photonics, 4, 056002(2022).

[3] Y. M. Xu, X. C. Pan, M. Sun et al. Single-shot ultrafast multiplexed coherent diffraction imaging. Photonics Res., 10, 1937(2022).

[4] A. Saba, C. Gigli, A. B. Ayoub et al. Physics-informed neural networks for diffraction tomography. Adv. Photonics, 4, 066001(2022).

[5] J. W. Miao, P. Charalambous, J. Kirz et al. Extending the methodology of X-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens. Nature, 400, 342(1999).

[6] J. M. Miao, T. Ishikawa, I. K. Robinson et al. Beyond crystallography: diffractive imaging using coherent X-ray light sources. Science, 348, 530(2015).

[7] R. W. Gerchberg. A practical algorithm for the determination of phase from image and diffraction plane pictures. Optik, 35, 237(1972).

[8] W. O. Saxton, J. M. Cowley. Computer techniques for image processing in electron microscopy. Phys. Today, 32, 74(1979).

[9] J. R. Fienup. Reconstruction of an object from the modulus of its Fourier transform. Opt. Lett., 3, 27(1978).

[10] J. R. Fienup. Phase retrieval algorithms: a comparison. Appl. Opt., 21, 2758(1982).

[11] R. Hegerl, W. Hoppe. Phase evaluation in generalized diffraction (ptychography). Proceedings of the Fifth European Congress Electron Microscopy(1972).

[12] J. M. Rodenburg. Ptychography and related diffractive imaging methods. Adv. Imag. Elect. Phys., 150, 87(2008).

[13] P. D. Nellist, B. C. McCallum, J. M. Rodenburg. Resolution beyond the ‘information limit’ in transmission electron microscopy. Nature, 374, 630(1995).

[14] H. N. Chapman. Phase-retrieval X-ray microscopy by Wigner-distribution deconvolution. Ultramicroscopy, 66, 153(1996).

[15] J. M. Rodenburg, H. M. L. Faulkner. A phase retrieval algorithm for shifting illumination. Appl. Phys. Lett., 85, 4795(2004).

[16] J. M. Rodenburg, A. C. Hurst, A. G. Cullis. Transmission microscopy without lenses for objects of unlimited size. Ultramicroscopy, 107, 227(2007).

[17] P. Thibault, M. Dierolf, A. Menzel et al. High-resolution scanning X-ray diffraction microscopy. Science, 321, 379(2008).

[18] J. M. Rodenburg, A. C. Hurst, A. G. Cullis et al. Hard X-ray lensless imaging of extended objects. Phys. Rev. Lett., 98, 034801(2007).

[19] S. Cao, P. Kok, P. Li et al. Modal decomposition of a propagating matter wave via electron ptychography. Phys. Rev. A, 94, 063621(2016).

[20] L. Valzania, T. Feurer, P. Zolliker et al. Terahertz ptychography. Opt. Lett., 43, 543(2018).

[21] L. Rong, F. R. Tan, D. Y. Wang et al. High-resolution terahertz ptychography using divergent illumination and extrapolation algorithm. Opt. Laser Eng., 147, 106729(2021).

[22] A. M. Maiden, J. M. Rodenburg. An improved ptychographical phase retrieval algorithm for diffractive imaging. Ultramicroscopy, 109, 1256(2009).

[23] A. M. Maiden, M. J. Humphry, M. C. Sarahan et al. An annealing algorithm to correct positioning errors in ptychography. Ultramicroscopy, 120, 64(2012).