Single-shot Phase Retrieval from
a Fractional Fourier Transform Perspective

In the realm of optics, the Fresnel diffraction equation emerges as an approximation of the Kirchhoff–Fresnel diffraction formula, tailored for characterizing optical propagation in the near field. The diffraction pattern denoted as $U_{d}$ can be calculated when light traverses through an object aperture represented as $U_{0}$ by

where $d$ denotes the propagation distance, $\lambda$ is the wavelength of light, and $i$ is the imaginary unit.

Normally, this Fresnel integral (1) can be expressed by the single Fourier transform, called SFT-Fresnel, written as

where $\mathcal{F}$ represents the two-dimensional Fourier transform.

Note that there is a quadratic phase factor in the Fourier transform and that the phase oscillates rapidly over short propagation distances, which poses a great challenge for accurate sampling and computing. Even worse, when employing Fast Fourier Transform (FFT) for discrete calculations, it suffers from serious aliasing artifacts [22].

As a dual version, this Fresnel integral (1) can be also seen as a convolution operation [31] and computed with the help of Fourier transform, stated as

where $H(f_{x},f_{y})=\mathcal{F}[\frac{e^{i\frac{2\pi}{\lambda}d}}{i\lambda d}e^{i\frac{\pi(x^{2}+y^{2})}{\lambda d}}]$ denotes the transfer function of Fresnel diffraction and $(f_{x},f_{y})$ are the Fourier coordinates conjugate to the real space coordinates $(x,y)$.

Moreover, this transfer function has an analytical expression as $H(f_{x},f_{y})=e^{i\pi d(\frac{2}{\lambda}-\lambda(f_{x}^{2}+f_{y}^{2}))}$. However, it is increasingly challenging to appropriately sample the transfer function, characterized by a swift oscillation of the phase component over long propagation distances. Numerical errors still exist when using FFT calculations duet to the fixed sampling pitch [31]. While this problem can be solved by utilizing non-uniform FFT [32], the computational complexity also increases, sacrificing the calculation efficiency.

In summation, it is exceedingly arduous to formulate an efficient model that can compute the whole near-field diffraction without sampling problems. In the following, we will present an innovative FrFT-based model to address the sampling problem in the Fresnel integral without aliasing errors during the treatment of the chirp function, featuring adaptability to both short-distance and long-distance propagation scenarios.

The fractional Fourier transform (FrFT) is the generalized form of the Fourier transform, and its kernel function is a quadratic phase term. The $p$th-order FrFT of a continuous signal $f(x)$ is defined as [33]

where $\mathcal{F}^{p}$ denotes the FrFT operator and $K_{\alpha}(u,x)$ is the transform kernel with $\alpha=\frac{\pi}{2}p$ given as follows

with $A_{\alpha}$ defined as $A_{\alpha}\triangleq\sqrt{1-i\mathrm{cot}\alpha}$ and $\delta(t)$ being the Dirac delta function. Especially, the FrFT reduces to the (inverse) Fourier transform when $p=\pm 1$.

Correspondingly, the inverse FrFT of $X_{\alpha}(u)$ in (4) is

where $\mathcal{F}^{-p}$ and $K_{-\alpha}(u,x)$ denote the inverse of $\mathcal{F}^{p}$ and the kernel $K_{\alpha}(u,x)$ obtained from (4) and (5), respectively.

Now, we consider introducing the two-dimensional FrFT of the input aperture $U_{0}$, which can be defined as

where $K_{\alpha}(x,x^{\prime})$ and $K_{\alpha}(y,y^{\prime})$ are transform kernels on the x-axis and y-axis, respectively.

To connect the Fresnel integral (1) and FrFT (7), we introduce the scaled fields $\hat{U}_{d}(x,y)\equiv U_{d}(s_{2}x,s_{2}y),\hat{U}_{0}(x^{\prime},y^{\prime})\equiv U_{0}(s_{1}x^{\prime},s_{1}y^{\prime})$ with $s_{1}=\sqrt{\frac{\lambda d}{\mathrm{tan}\alpha}},s_{2}=\sqrt{\frac{\lambda d}{\mathrm{sin}\alpha\mathrm{cos}\alpha}}$.

Then we can obtain the scaled diffraction field

Considering the intensity-only measurement, the spherical phase factor of (10) disappears and we can have

Thereby, we conclude that the scaled amplitude distribution of the near-field diffraction can be interpreted as the continuous FrFT magnitude. Given the scale factor $s_{1}$ and the physical propagation distance $d$, we can know exactly the other scale factor $s_{2}=s_{1}\sqrt{1+(\lambda d)^{2}/s_{1}^{4}}$ and the corresponding fractional order $p=2/\pi\times\mathrm{arctan}(\lambda d/s_{1}^{2})$. It can be seen that when the propagation distance $d$ increases, the corresponding fractional order $p$ gradually becomes larger with the range is [0,1]. This is consistent with the Fourier transform corresponding to propagation into the far field.

Discrete Computing. In practical scenarios, the acquisition of the diffraction field is inherently limited to sampling. Therefore, accurate discretization of both the input field and the transform kernel is crucial for numerical calculations to avoid aliasing artifacts. Benefiting from previous studies, various types of discrete fractional Fourier transform (DFrFT) have been developed with distinct strategies and properties [27] that could be applied here. Prominent examples include the eigenvector decomposition-type DFrFT (ED-DFrFT) [34], the improved sampling-type DFrFT (IP-DFrFT) [35], and the closed-form sampling-type DFrFT (CF-DFrFT) [36]. Specifically, ED-DFrFT offers orthogonality, additivity, and reversibility properties at the expense of high computational complexity while IP-DFrFT enjoys high discrete accuracy and can be implemented efficiently with lacking unitarity. CF-DFrFT achieves reversible properties by carefully considering sampling interval limits and involves low-complexity calculations in $O(NlogN)$ time. Nevertheless, it is essentially similar to SFT-Fresnel and also faces the problem of numerical errors. Therefore, we select the IP-DFrFT approach in this work and present some details as follows.

Specifically, the samples of the transformed function in (11) spaced at the interval $\triangle x$ and $\triangle y$ are obtained as

where $m,n$ goes from $-N/2$ to $N/2$, and $N$ is the sampling number.

Following the computation scheme of IP-DFrFT, we can get

where $\triangle x^{\prime}=\triangle x,\triangle y^{\prime}=\triangle y$ is the spacing interval and $m^{\prime},n^{\prime}$ represents the discrete grid in the source field.

Its calculation can be recognized that the discrete source field is first modulated with a chirp function $e^{i\pi(\mathrm{cot}\alpha-\mathrm{csc}\alpha)(m^{\prime 2}\triangle x^{\prime 2}+n^{\prime 2}\triangle y^{\prime 2})}$ and then convoluted by another chirp function $e^{i\pi\mathrm{csc}\alpha[(m^{\prime}\triangle x^{\prime})^{2}+(n^{\prime}\triangle y^{\prime})^{2}]}$ which can be efficiently implemented by FFT. Note that it can avoid the aliasing error due to the rapid oscillations of the kernel by exploiting the periodicity and additivity of the continuous FrFT¹¹1Note that there is an assumption as $0.5\leq|p|\leq 1.5$. Taking advantage of the additivity property of FrFT, we can extend the range of parameter $p$ to cover all its values. For example, for the range $0<p<0.5$, we have that $\mathcal{F}^{p}=\mathcal{F}^{p-1+1}=\mathcal{F}^{p-1}\mathcal{F}^{1}$.. Given this, we can get an efficient, accurate, and unified method to calculate the whole diffraction field. It is worth mentioning that it is necessary to perform the dimensional normalization on the fractional order $p$ and the scale factor $s_{2}$ during the numerical calculations, in order to eliminate the sampling-related influences:

where $L$ is the length of the input aperture.

Scalable Sampling. Considering the scaling operator introduced in (11), it can be addressed by obtaining the rescaling samples in practice. Specifically, the scale factor $s_{1}$ can usually be set to 1 in order to be consistent with the real input field, i.e., $\hat{U}_{0}=U_{0}$. Thereby, the original near-field measurement can be obtained by

where $s_{2}=\sqrt{1+(\lambda dN/L^{2})^{2}}$ can be seen as a sampling rate conversion factor.

Consistent with previous considerations, the sampling interval within the FrFT domain remains $\triangle x^{\prime}$, consequently leading to the pixel pitch of $s_{2}\triangle x^{\prime}$ within the real observation plane. This, in turn, enables us to directly acquire the intensity measurements of the FrFT through near-field diffraction. To clarify, Table I shows a summary of the methods mentioned within this work, encompassing pertinent details such as the computational complexity, pixel pitch, and the respective suitable range. In addition, we can also adopt digital computation in the FrFT domain to eliminate the scaling operator. It can be noticed that the scaling factor $s_{2}$ gradually becomes larger as the propagation distance $d$ increases and is always greater than 1, which corresponds to the upsampling case. According to the digital sampling rate conversion, we can implement it through interpolation, low-pass filtering, and decimation in sequence, achieving flexible sampling intervals. Overall, the proposed FrFT-based measurement model is illustrated in Fig. 2.

Based on the proposed FrFT-based measurement model, the SFrFPR problem can be mathematically stated as

where $I$, $\mathcal{F}^{p}$, and $O$ represent the intensity of the diffraction pattern, the corresponding $p$th-order FrFT measurement model, and the underlying object, respectively.

Typically, the inverse problem of SFrFPR can be solved within a regularized optimization framework by minimizing the following cost functional:

where $\mathcal{R}(O)$ indicates the regularization term associated with the prior knowledge of objects, and $\beta$ is a parameter that controls the weight of the regularizer.

This is a non-convex and non-linear ill-posed problem, mainly caused by the loss of phase. Existing iterative optimization PR methods such as Wirtinger gradient descent [10] and plug-and-play methods [38, 30], are expected to provide transferable ideas to solve this problem. However, these methods all rely on accurate forward and backward projections, resulting in the use of ED-DFrFT. The high computational complexity prevents its application to large-scale and real-time reconstruction. The fast discrete FrFT algorithm represented by IP-DFrFT can solve this dilemma, but it lacks unitarity and suffers from numerical errors during the inverse transformation [27], making it unsuitable for these iterative methods.

To address this issue, we put forward an alternative solution, leveraging the concept of the deep image prior (DIP) [44] and employing an untrained neural network (UNN) for SFrFPR. The key idea is designing a neural network $f_{NN}(I,\Theta)$ to directly perform the inverse mapping from captured intensity measurement $I$ to the desired object $O$ by adjusting the network’s weight $\Theta$ based upon the following empirical risk:

Compared with (17), the main difference is that we optimize the neural network’s parameters instead of the estimated object directly. This direct problem transformation has brought many benefits. First of all, thanks to the development of neural network technology, this optimization process can be well solved by the auto-differentiation technique. In this way, the proposed method solely requires the differentiability of the forward function while avoiding the need for the backward function. In addition, the proposed method operates in a self-supervised manner, utilizing the network’s weight adjustment to reconstruct the desired amplitude or phase object, guided by the captured intensity measurement and incorporating the FrFT-based measurement model. Therefore, the proposed method does not rely on obtaining a large number of pairs of ground truth data and corresponding observations.

Overall Pipeline. The overall pipeline of our method is outlined in Fig. 3. The input to the neural network is a diffraction pattern of an amplitude or phase object, captured in a single snapshot. The neural network processes this input and generates an estimated object. Subsequently, the estimated object is numerically propagated to simulate the resulting diffraction pattern using the proposed FrFT-based measurement model. Engaging a loss function computed between the measurement and the estimated diffraction pattern, the parameters of the neural network are adjusted via the auto-differentiation technique. Notably, this training process only involves the FrFT forward function and unlabeled simulated/measured diffraction patterns. After updating, the trained network can perform the direct inversion from a single intensity pattern to the real space object without requiring the iterative process.

Network Implementation. In UNN, the architecture of the neural network is very important, because it will introduce the neural network structure prior as a regularization term. Witnessing the powerful representation ability of the Transformer network in large models, we explore its potential as an untrained network for SFrFPR. To be specific, we adopt a general U-shaped Transformer architecture, as delineated by [45], and build a small hierarchical encoder-decoder network with a tiny computational burden. Given a snapshot pattern as the network input, a convolutional layer with LeakyReLU is used to extract low-level features. Subsequently, these features traverse through two encoder phases, each of which encompasses a LeWin Transformer block and one down-sampling layer. The LeWin Transformer block captures long-range dependencies via non-overlapping windows instead of global self-attention, resulting in the low computational cost of high-resolution feature maps [45]. Proceeding further, a bottleneck phase characterized by a LeWin Transformer block is appended to culminate the encoding progression. Then, two decoder phases are followed to recover the features, each of which contains an up-sampling layer and a LeWin Transformer block. Finally, a convolutional layer is applied to output the reconstructed object.

While Transformer models have undoubtedly demonstrated remarkable performance in numerous supervised learning tasks, their untapped potential as untrained neural networks remains to be further explored. To the best of our knowledge, this is the first attempt to introduce the Transformer-based architecture into an untrained neural network for solving the PR problem. In Section III, experimental results demonstrate that the Transformer structure priors leverage both local and global dependencies with better reconstruction performance than convolutional neural networks in a self-supervised scheme.

To facilitate a comprehensive assessment of the proposed method, we conducted a computational simulation involving the diffraction propagation of a two-dimensional rectangular aperture under typical physical parameter configurations. Specifically, the spatial length of the sampling windows and the number of samplings are $1000$ um and $512$ in the source and destination plane, respectively. The width of the rectangular aperture and the wavelength are $500$ um and $500$ nm, respectively. The range of propagation distance is $1\sim 50$ mm. The detailed parameter values in numerical calculations are listed in Table II. Moreover, the proposed FrFT-based measurement model was realized on Pytorch, thus embracing the modern GPU acceleration.

To verify the accuracy of the proposed FrFT-based measurement model, we conduct two-dimensional diffraction computations using three distinct methods, i.e., single Fourier-transform-based Fresnel model (SFT-Fresnel), Fresnel transfer function model (Fresnel-TF), and the proposed FrFT-based measurement model (FrFT). For quantitative analysis, the reference field is provided by numerical integration of the Fresnel integral (1) via the trapezium rule. And we achieved the scalable sampling via post-processing to ensure consistent resolution and pixel pitch size for comparison.

The accuracy of each method is assessed by comparing the peak-signal-to-noise ratio (PSNR) of the diffraction pattern against the reference, with the results presented as a function of propagation distance in Fig. 4. It can be observed that the accuracy of the Fresnel-TF method deteriorates with increasing propagation distance, while the SFT-Fresnel method is unsuitable for numerical propagation over short distances. In contrast, the proposed FrFT method demonstrates remarkable accuracy and features suitability for both short-distance and long-distance propagation. An illustrative examination of the diffraction intensity patterns computed through these methods is presented in Figure 5. Notably, the SFT-Fresnel method manifests pronounced numerical errors within regions corresponding to short propagation distances. Conversely, while the Fresnel-TF method shows good performance when close to the source field, numerical errors become increasingly apparent at longer propagation distances. In contrast, the proposed FrFT-based model demonstrates consistency with numerical integration, providing similar results over the entire propagation distance range without apparent aliasing errors in both short- and long-range scenarios.

The proposed reconstruction method was implemented based on the PyTorch version $2.0.0$ platform with Python $3.9$ via one Nvidia GeForce GTX 1080 Ti GPU. During the training process, the model was optimized using the Adam optimizer for a total of $10000$ epochs, with the learning rate empirically set at $2e^{-4}$. Notably, our method does not require any external training data apart from the input measurements themselves. To evaluate the effectiveness of the proposed method, we utilized two testing datasets Set12 [48] and Cell8³³3This dataset comprises a collection of eight distinct cell images, namely Brown fat cell, Peritoneal macrophage, Intestinal epithelial cell, Epithelial cell, Blood cell, Auditory hair cell, Acinar cell, and Chromoffin cell. These original cell images are available for download from the website http://www.cellimagelibrary.org/browse/celltype.. For the sake of unification, all images were resized to a uniform dimension of $256\times 256$.

To verify the performance of the proposed method, we mainly compare it against one classic PR approach, namely Wirtinger Flow (WF) [10], two plug-and-play approaches GAP-tv [30] and prDeep [38], and two state-of-the-art untrained neural network (UNN) approaches PhysenNet [43] and DeepMMSE [41], which are implemented to be extended to solve the SFrFPR problem. Among them, various hyperparameters of these algorithms are set to the optimal. Here we further consider the reconstruction of amplitude-only and phase-only objects from measurements in different fractional Fourier orders⁴⁴4prDeep is limited to real-valued reconstruction, so phase-only objects are not considered here [38].. Table III reports the average recovery accuracy in terms of mean peak-signal-to-noise ratio (PSNR) and structure similarity index measure (SSIM) on two testing datasets Set12 and Cell8. It can be observed that optimization-based iterative algorithms such as WF, GAP-tv, and prDeep, can effectively reconstruct objects from some fractional Fourier measurements ($p=0.5,0.6,0.8$) while all fail to recover images from the Fourier transform measurement ($p=1$). Unfortunately, due to the inverse transform error of the fast fractional Fourier transform discrete algorithm, these iterative algorithms cannot achieve good results at some specific FrFT orders ($p=0.2,0.4$). On the contrary, all UNN methods consistently perform well on FrFT measurements with different fractional orders except $p=1$ which denotes the Fourier transform measurement. Moreover, the reconstructing performance can be dramatically improved by the proposed transformer-based UNN when using the FrFT measurement.

For visual comparison, Fig. 6 and Fig. 7 present the reconstruction results of different PR methods on an amplitude-only and phase-only object from various FrFT measurements, respectively. It can be seen that all UNN methods can reconstruct the satisfactory amplitude and phase of objects from the single FrFT measurement but fail from the Fourier transform measurement. Note that the proposed method can significantly improve the reconstruction quality. In addition, the classic iterative method WF can also reconstruct the object from a single FrFT measurement with a suitable order, which cannot be achieved in the Fourier case.

To further investigate the effectiveness of the proposed method, we show the convergence behaviors of the proposed method from different FrFT measurements. Fig. 8 presents that the measurement loss calculated by mean square error (MSE) and reconstruction quality indicated by PSNR varies with training epochs of the proposed method on testing an amplitude object. It can be found that the measurement loss can gradually decrease and converge to a very small value for each FrFT order. However, the proposed method suffers from serious stagnation and produces a poor reconstruction quality (low PSNR) from the Fourier transform measurement. On the contrary, the FrFT measurement can be well-combined with the UNN priors to effectively avoid the stagnation problem and improve the reconstruction quality with increasing PSNR.

In order to verify that the proposed FrFT measurement can effectively overcome the ambiguity problem, we use two special but representative initializations to illuminate it. Specifically, we shift and flip the original image adding Gaussian random noise as the input of reconstruction methods. Then we utilize the proposed UNN-based method to reconstruct the original image from its Fourier transform ($p=1$) and FrFT ($p=0.5$) measurement, respectively. Fig. 9 presents the reconstructed results of the intermediate reconstruction process. It can be found that the proposed method can gradually reconstruct a correct object by eliminating the effect of translation or inversion from the FrFT measurement. On the contrary, the reconstruction method suffers from serious stagnation and produces an inaccurate solution from the Fourier transform measurement. While such trivial ambiguities are probably acceptable, they will compete with the correct solution and confuse the reconstructing algorithms in practice. Thus, the removal of ambiguities via the FrFT measurement can greatly alleviate the stagnation problem and improve the reconstruction performance.