Far-field diffraction computational imaging based on parameter-robust illumination and direct phase optimization

The Fraunhofer diffraction or far-field diffraction can be easily achieved in most scenarios, including visible light, X -ray and electron[1]. CDI utilizing far-field diffraction is widely applied in observing structure of microorganism[2, 3, 4, 5], molecule[6, 7] and material[8, 9, obinson2001reconstruction]. The mathematical form of the diffraction field is concisely the Fourier transformation of the light field at the object plane. However, because the phase part of the light field is hardly directly measurable[1], it needs to be retrieved for the Fourier inversion. In practice, as a multi -solution problem, phase retrieval is usually conducted in two ways. One is multiple measurements: several masks cast on the target[10, 2, 11] or overlap scanning such as ptychography[12, 13, 9, 14] to provide more information. The other way is to use certain generals of the object field in specific application scenarios, such as real-valued condition[15, 6, 16, 17], sparsity in certain transform[18, 19, 20, 21] and deep learning methods[22, 23]. Multiple measurement needs to ensure that the object will not change between the measurements, which limits dynamic imaging and high energy imaging[24]. Utilization of priori information surely cannot faithfully reconstruct the object that violate these characteristics, for example, a random complex-valued field. If the far-field diffraction measurement and computational reconstruction is taken as a black box, it is hoped to be as simple as a mirror or lens: simply reflects any original field without multiple measurements or priori information. Efforts have been made for the ideal black box and achieve good simulation result[fienup1978reconstruction，fienup1982phase, 25]. But when applied in practice, extra restriction seems indispensable. Intensity object is still more welcome in experiments[8, 15, 26].

The phase retrieval of single-shot far-field diffraction without prior information is inherently a problem with multiple equivalent solutions[27]. Some equivalent solutions cannot be excluded by simple oversampling, including:

• •

  overall phase shift;

• •

  translation paired with phase modulation that matches the moving distance;

• •

  inversion of complex conjugate;

• •

  combinations of the operations above.

In the presence of higher noise, the resulting final phase retrieval solution tends to be degraded by superimposed mixture of the equivalent solutions.

In order to suppress the degradation, there have been some good simulation works using special illumination pattern to destroy certain symmetry[28, 29]. These methods can also be classified into adding general information about the object field. Because the special illumination is usually built-in the set-up, the information remains unchanged, providing more feasibilty for different samples. But it appears difficult to achieve satisfying result in practical applications. Special field restriction is still more viable for intensity sample[15, 30] . The perfect pixel mapping and size parameters of the pattern and object can be easily obtained in the simulation calculation. However, due to the optical path distortion, aberration, fluctuation of refractive index, photosensitive array spacing and other practical problems in experiment, the parameter and pixel mapping used in the reconstruction algorithm are hardly as precise as those in simulation. These errors can greatly degrade the retrieval result under high noise as shown in Fig. 2(b). Here, we propose a class of parametric-robust apertures and a direct phase optimization algorithm for the corresponding special-shaped apertures, together with a measurement correction algorithm.

For phase retrieval without priori information, all we can rely on is that the intensity of the light field outside the passable region is zero. In the optimization procedure, if the imagined aperture is smaller than the true one, the actual light intensity in the assumed zero zone will not be zero, which will mislead the algorithm to solve an impossible problem. Therefore, it is necessary to ensure that the theoretical aperture can cover the actual aperture. Due to parameter errors, the ideal aperture needs to be slightly enlarged. . We choose linear magnification operation to discuss hereinafter. The reason is: 1.Operations such as dilation and convolution will lead to the emergence of translation equivalent solutions; 2.Errors like boundary error on Fourier plane will lead to linear zoom on image plane; 3.The linear magnification operation retains most of the symmetry embodied in the measurement.

The linear magnification operation is: for a enlarge center point $O$ and any point $A$ in the space, the corresponding new point $A^{\prime}$ needs to satisfy $\vec{OA^{\prime}}=a*\vec{OA}\ \ (a>1)$. We require that for any $a>1$, the enlarged space covers the original one, there is $\{A\}\subset\{A^{\prime}\}$. Let the boundaries be $\{L_{i}\}$ for the aperture space. These $\{L_{i}\}$ are Jordan curves. Jordan curve theorem will be used a lot below: every continuous path connecting ‘outside’ and ‘inside’ points intersects with the curve somewhere.

For $O$ outside a certain $L_{i}$, if $O$ is inside $L_{i}\textquoteright$, because $O=O\textquoteright$, $O\textquoteright$ is inside $L_{i}\textquoteright$; $O$ must be inside $L_{i}$; but $O$ is outside $L_{i}$; contradiction occurs. Therefore, If $O$ is outside $L_{i}$, $O$ must be outside $L_{i}\textquoteright$.

If $O$ is outside a certain $L_{i}$, $Q$ is a point inside $L_{i}$, the nearest intersection of line segment $OQ$ and $L_{i}$ to $O$ is $A$(Fig. 1(b)). If $A$ is on or inside $L_{i}\textquoteright$, because $O$ must be outside $L_{i}\textquoteright$, $OA$ and $L_{i}\textquoteright$ should intersect at some point $G\textquoteright$. $OA>OG\textquoteright>OG$ ($OG$ and $OG^{\prime}$ must be in the same direction, outerwise $O$ is inside $L_{i}\textquoteright$(Fig. 1(c))), and $G$ is on $L_{i}$, so $A$ is not the nearest intersection, which is a contradiction. So $A$ is outside $L_{i}\textquoteright$, which means $L_{i}$ cannot be completely surrounded by $L_{i}\textquoteright$. For multi-connected graphs, there may be an a large enough that $L_{j}^{\prime}$ can surround $L_{i}$. But we need $a$ to be any number more than 1. For a multiply connected graph, for any $O$, there is always some $L_{i}$ that $O$ is outside $L_{i}$. Therefore, only a simply connected graph can meet our needs.

And for this simply connected graph, we require that: there is at least one $O$ inside $L$ that any ray from $O$ intersects $L$ at only one point. We call $O$ a simple encircled center. For this kind of graph, suppose that $C$ is in $\{A\}$ but not in $\{A^{\prime}\}$, let $B_{1}^{\prime}$ and $B_{2}^{\prime}$ be the intersection of line $OC$ and $L^{\prime}$ , $OB_{1}^{\prime}$ and $OC$ are in the same direction. If $B_{1}$ does not belong to $\{A\textquoteright\}$, $B_{1}O$ intersects $L\textquoteright$ on $G\textquoteright$, $OG<OG\textquoteright<OB_{1}$, $G\neq B_{1}$ ,contradicting simple encircling. Hence $B_{1}\in\{A^{\prime}\}$. That means, there is at least one intersection $D^{\prime}$ on $B_{1}C$ and $L^{\prime}$ . The point $D$ corresponding to $D^{\prime}$ is the intersection of line $OC$ and $L$ , contradicting simple encircling (Fig. 1(a)). Therefore, for a simply connected graph that meets the above conditions, its linear magnification of $O$ according to any $a>1$ can completely cover the original graph. Incidentally, not every simply connected graph has a simple encircled center, for example, a slim $\Omega$.

In order to satisfy the requirements above, we can generate several sectors centre on point $O$ to generate the desired aperture shape. The radius, angle and number of these sectors are not limited, which means shapes like rectangle and triangle can also be generated.

In addition to the requirement that linear magnification fully covers the original one. The aperture shape had better have the following properties:

• •

  It has a small self-convolution maximum value (corresponding to the maximum overlap with the inversion after translation);

• •

  When enlarged, it limits the possible translation of the original graphics in it;

• •

  Continuous area in it is as large as possible so that the final reconstructed image is more identifiable.

Based on the conditions described above, we choose the aperture of the equal-radius three-sector structure as an example below. The maximum value of self-convolution for three-sectors is 0.484 while that for five-sectors only decreases to 0.477. The maximum self-convolution value for more number of odd-numbered sectors even increases.

Fig. 1 shows the noiseless performance of multiple connected non-centrosymmetry aperture and simple connected three-sectors aperture. When the assumed aperture is exactly the actual one, both of them works well in attaining the original image. However, when the assumed aperture is slightly enlarged, the result is dramatically degraded in the multiple connected case, because the enlarged pattern cannot cover the original one. While three-sectors after enlarged contains the original graph, leading to robustness towards size error.

In this paper, no prior information about the target object is given. The only information is the intensity of the light field outside the aperture is zero or as low as possible. Without the assistance of other information, traditional phase retrieval algorithms such as alternative-projection and gradient algorithms are difficult to stably converge to a better solution for complex-valued object under high noise[27, 31]. Let the light field at the wave vector $k$ in the Fourier plane be ${m_{k}}{e^{i{\theta_{k}}}}$. Many algorithm works take ${e^{i{\theta_{k}}}}$ as a whole to optimize the loss value, so that many linear optimization methods can be utilized[32, 33, 34]. But in the end, The quadratic condition that the modulus of ${e^{i{\theta_{k}}}}$ is 1 has to be used, which is the principal trouble. In order to maximize the use of aperture information, we take $\theta_{k}$ as the optimization object here.

The definitions that will be used later are collected here to be checked.

Let $Z$ denote the set of points unilluminated in object plane.The loss value for the aperture information is expressed as

and its partial derivative with respect to $\theta_{k}$ is:

Since $\gamma_{k}$ does not contain $\theta_{k}$ (Eq. (9)),

The value of $\theta_{k}$ at the extreme point can be directly solved by the condition that the partial derivative is 0. Let

At the extreme point

Because

At the local minimum, we have

An equation for $\theta_{k}$ at the local minimum is already achieved here. However, if $\gamma_{k}$ is calculated separately for each wave vector pixel in the Fourier measurement, the overall phase retrieval may require millions of Fourier transforms, so we need to go further to reduce the calculation amount.

Let

Here $\hat{F}^{\prime}$ and $\hat{f}^{\prime}$ are two demensional FFT (Fast Fourie Transform) and iFFT (inverse Fast Fourie Transform).

Since it is a multivariate problem, the solution obtained after each batch is the extreme point of the previous batch of phases, so it still takes several iterations to converge to the extreme point. For each pixel, if only one $\theta_{k}$ is changed at a time, the new overall loss value obtained after each calculation will only be less than or equal to the value of the previous batch. Therefore we can have confidence in the convergence to this problem. Of course, considering the actual calculation efficiency, it is not necessary to be too cautious here, and the phase value can be updated more comprehensively.

In high noise, the extreme point is not necessarily the global lowest point. Since we are confident in the convergence, we can get away from local minima by resetting the phase values of a large percentage of random locations to random phases when the solution process stalls.

On the other hand, we can also optimize the measurement to some extent in a similar way. (loss to mk)

The measurement value at extreme point under the current phase can be directly obtained utilizing the linear expression of $M$. However, the direct calculation will lead to a wrong measurement when a wrong phase distribution is given. Error in phase cause error in measurement; error in measurement cause error in phase, causing error to be fixed in subsequent iterations. So only limited optimization can be implemented based on the existing gradient. At the same time, because loss value corresponding to all 0 measurement is 0, there is always a tendency to reduce the overall measurement intensity. Hence, change on sum of the measurement values must be restricted. So we adopt the following update strategy here to assign the change value according to the gradient.

Where $q$ is the the variable controlling the overall gradient. In addition, we limit the maximum change during the iteration process to make measurement as faithful as possible to the original one.

In Fourier measurements, there is a significant difference between the intensity in the low -frequency area and the high -frequency area. The traditional Signal-Noise Ratio (SNR) can hardly reflect the impact of noise on the object plane. We define Inverse Fourier Signal-Noise Ratio (IFSNR):

Here ${m^{o}}$ is the noiseless measurement; ${m^{r}}$ is the noisy measurement; $\theta^{s}$ is the supplementary phase. The supplementary phase can be selected as the original Fourier phase, random phase, or the phase of the retrieval. In the simulation, the supplementary is the phase of the noise -free original Fourier transform .

Fig. 3 shows the performance of OPM and OMP under high noise measurement. The procedure in this simulation is OPM $\rightarrow$ OMP $\rightarrow$ OPM $\rightarrow$ OMP. As shown in the Fig. 3(a)(c), even though the high -frequency part has been severely damaged, OPM can still obtain a relatively good result. While OMP also reconstructs many of the destroyed characteristics in high -frequency area.

The simplified optical path of the set-up is shown in Fig. 4.

The light source is a polarized He-Ne laser light source. The aperture is realized by a pattern loaded on the digital micromirror device (DMD). The aperture is imaged on the spatial light modulator (SLM) through a convex lens. Complex-valued object is a pattern on the SLM. The far-field diffraction pattern is measured with a 14-bit cmos at the focal plane of the lens after SLM reflection. The exposure time is 20ms.

The measurement on CMOS contains some structured noise, hence we conduct OMP first with strong restriction to decrease the noise out of aperture Fourier’s transformation. The supplementary phase here is the mean of the Fourier phases of several random patterns with different granularity.

Distortions on the light field that actually projected on the target cannot be perfectly described by linear magnification. Therefore, we first take an assumed aperture with larger angle and radius for reconstruction, and obtain a new aperture based on the reconstruction. Several steps of erosion and dilation[35] are conducted to eliminate scattered pixels to meet the requirements of simple connection.

The coordinate of the Fourier transform center is hardly an just integer on CMOS, result in ripple phase distribution correlated with the offset. Meanwhile, due to surface unevenness, aberration and oblique light path, there is a base phase pattern on the target. These problems are solved by subtracting the reconstructed phase of a full 1 target, in other word, a mirror.

The experiment result is demonstrated in Fig. 5. The OPM faithfully reconstruct the phase morphology of the complex-valued object, showing more details than Hybrid Input Output (HIO). The pillars are also better distinguishable as shown in Fig. 5(e).

In summary, we demonstrate that for the parametrically robustness of special illumination method in CDI, the aperture pattern should be simple connected and have a simple encircled center. Three-sector scheme is raised as an appropriate choice. Meanwhile, we propose an algorithm with Fourier phase as the direct optimization object and a coupling measurement correction algorithm. These methods make it possible to achieve phase retrieval with good convergence for complex-valued objects in high-noise measurements. Finally, the experiment result demonstrates the feasibility and efficacy of our scheme and algorithm.

Our aperture scheme only needs blocking and transmitting, hence it can be easily achieved in various wavelengths. Because no information except for the aperture is required in our algorithm, the target is not necessary to meet certain nature such as real-valued and sparsity. Therefore, the algorithm has good compatibility with various application scenarios. A single-shot imaging without multi -frame modulation enables it to be applied to high-energy destructive imaging and imaging of moving objects. The better imaging performance of complex-valued object is of great use in 3D imaging and non-staining biological imaging.

In the algorithm, we use random resetting method to avoid local optimum. Simulated annealing without cooling down consumes a lot calculation time. Better optimization is anticipated here. Subsequent jobs might take deep learning and other methods to fill up the overexposure. Combination of prior knowledge should bring better image quality in specific practical scenarios.

Funding funding

Acknowledgments We thank Rong-ping Deng, Wen-kai Yu, Zu-xie Hu and Biao-xv Peng for their help in writing and theory assistance.

Disclosures The authors declare no conflicts of interest.

Data Availability The data that support the plots within this paper and other findings of this study are available from the corresponding author upon reasonable request.