Compensation for the setup instability in ptychographic imaging

The high-frequency vibration of the imaging system degrades the quality of the reconstruction of ptychography by acting as a low-pass filter on ideal diffraction patterns. In this study, we demonstrate that by subtracting the deliberately blurred diffraction patterns from the recorded patterns and adding the properly amplified subtraction to the original data, the high-frequency components lost by the vibration of the setup can be recovered, and thus the image quality can be distinctly improved. Because no prior knowledge regarding the vibrating properties of the imaging system is needed, the proposed method is general and simple and has applications in several research fields.


Introduction
Coherent diffraction imaging (CDI) is a promising technology for obtaining the complex transmission function of a specimen from the recorded diffraction intensity. As a lens-free technique, CDI can bypass the resolution limits imposed by the poor focusing optics available at short wavelengths [1,2] and can theoretically reach the diffraction-limited resolution. With X-ray and high-energy electrons, a resolution of nanometers or angstroms can be achieved; thus, CDI is becoming an important tool in material and biological sciences [3][4][5][6][7]. Because the performance of traditional CDI algorithms is not very satisfying with regard to convergence, accuracy, and reliability, several improved CDI methods have been proposed [8][9][10][11]. The ptychographic iterative engine (PIE) [12] is a scanning version of the CDI technique where the specimen is scanned through a localized illumination beam to a grid of positions and the resulting diffraction patterns are recorded. Using an iterative scheme with a proper overlapping ratio between two adjacent scanning positions, the modulus and phase of the transmission functions of the specimen and the illumination beam can be reconstructed accurately and rapidly [13,14]. In theory, PIE can easily generate images with a resolution only limited by the numerical aperture of the detector; however, in practice, the resolution is affected by the flaws of the imaging system, especially for imaging with X-ray and electron beams. The coherence of synchrotron radiation sources is not as good as that of a laser beam and has been proven to be the largest barrier for PIE to reach the theoretical diffraction-limited resolution [15,16]. Several algorithms have been investigated to improve the image quality of PIE with partially coherent illumination [17][18][19][20][21][22]. Although these methods require complete or partial prior knowledge of the properties of illumination or timeconsuming computation and may slightly compromise the spatial resolution, the quality of the reconstruction can be distinctively improved in many cases. Furthermore, the recently developed technique of Fourier ptychography has potential to avoid the influence of the incoherence for achieving high-quality reconstruction [23,24]. Another factor that limits the practical resolution of PIE is the instability of the imaging system, including the vibration of the mechanical scanning system, the tiny pointing-direction change, and the transverse shifting of the radiation beam [25,26]. Because the wavelengths of the X-ray and the electron beams are much smaller than 1 nm, the tiny departure of the illumination beam from the right position and direction can generate obvious errors in the final reconstruction. Numerous methods have been proposed to efficiently correct the low-frequency vibration with a period far longer than the exposure time of the detector [27-30]. However, there is no ideal method for effectively handling the highfrequency vibration of the experimental system, which makes the recorded diffraction intensity a summation of many diffraction intensities formed by the changing illumination during the exposure of the detector [31]. Although the influence of the high-frequency vibration can be treated as a type of incoherence of the illumination, the aforementioned methods [17][18][19][20][21][22]31] require complete or partial prior knowledge of the properties of vibration-including the frequencies and amplitude-or massive data processing. However, in practice, the parameters of the vibration of the whole imaging system are difficult to measure in real time. To deal with the image-quality degradation induced by the setup instability, we must examine its influence on the recorded data and the final reconstruction and then develop a method to circumvent this problem.
In this study, we mathematically analyze how the recorded diffraction patterns and the final reconstruction are blurred by the instability of the imaging system and then propose a simple numerical method for enhancing the lost high-frequency components. The proposed method does not require any prior knowledge regarding the characteristics of the high-frequency vibration of the imaging system and can be extended to other CDI methods using X-ray and electron beams to solve the problems related to imaging-system instability. In the PIE method, the specimen O(r) is fixed on a two-dimensional (2D) translation stage and illuminated by a localized probe P(r). Assuming that the specimen is sufficiently thin, the exiting wave from the specimen is ϕ(r, R) = O(r − R)P(r), and the recorded diffraction intensity is I(k, R) ∝ |ℑ[ϕ(r, R)]| 2 in most experiments with short-wavelength sources, where k is the reciprocal coordinate with respect to the real space coordinate r in the specimen plane, and R denotes the position of the specimen during the raster scanning. With the recorded diffraction patterns, the complex amplitudes of the specimen and the probe can be reconstructed. The flowchart of the reconstruction process of PIE is shown in Fig. 1, and a detailed description is found in the literature [12,13].

Influence of the setup instability
In the PIE algorithm, the illumination probe is assumed to be absolutely static during the data acquisition; however, its direction and position change continually within a small range with respect to the specimen and detector owing to the instability of the imaging system. In most cases, the direction and position of the illumination beam change simultaneously during the exposure of the detector, but for simplicity, we analyze them separately to determine how the recorded data are influenced mathematically and then consider their combined effects.
According to the principles of Fourier optics, any illumination beam can be decomposed into a series of spherical waves of different strengths; thus, we can assume a spherical illuminating probe without a loss of generality. The spherical illumination is expressed in Fourier form as P(k) = W(k)exp(− jπλzk 2 ), where λ is the wavelength of the probe, and z is the distance between the focal spot of the probe and the specimen, and W(k) is determined by the numerical aperture of the illuminating optics. For a short wavelength or in the far field, the diffraction distribution in the detector plane is the Fourier transform of the wave exiting the specimen. The complex amplitude of the diffraction pattern in the detector plane can be expressed as where O(k) is the Fourier transform of the transmission function of the specimen. The recorded intensity with static illumination is The first item I 0 (k) is the intensity summation of all diffracted beams, and the second item m n A m (k)A n (k)cos[θ(k)] represents the interference between different spatial-frequency components.
Here, k m − k n is the frequency of the interference fringe formed by the m th and n th diffraction orders, and φ m,n is the additional phase introduced by the specimen. Considering the coordinate transformation r c = λLk on the CCD plane, where L is the distance between the specimen and the CCD, the diffraction pattern intensity can be described in the frequency domain as where u is the reciprocal coordinate with respect to the real space coordinate r c in the CCD plane.
When considering the pointing instability of the imaging system, the illumination beam tilted an angle of α can be expressed as, So the Fourier transform of the probe is P ′ (k) = P(k − ∆k), where ∆k = sinα/λ. The recorded intensity with tilted illumination becomes Assuming that the vibration in the pointing direction of the illumination beam follows the normal distribution H 1 (∆k) = exp(−∆k 2 /K 2 ), where K is a constant related to the vibrating properties of the imaging system. The recorded intensity with high-frequency vibration in the pointing direction of the illuminating beam can be interpreted as a summation of the diffraction patterns of all possible illuminations with different pointing directions. Thus, it can be expressed as Considering the coordinate transform r c = λLk in the detector plane, the recorded intensity can be expressed in the frequency domain as Thus, the influence of the vibration in the pointing direction of the illuminating beam acts as a low-pass filter on the ideal diffraction patterns. On the other hand, when the illumination probe suffers from transverse positioning vibration, the recorded diffraction pattern can be expressed as a summation of the diffraction intensities formed by all possible illumination beams with diverse transverse shifts. The probe with a transverse shift of δ is expressed as P ′′ (r) = P(r + δ). And its Fourier transform is P ′′ (k) = P(k)exp( j2πkδ). The corresponding diffraction-pattern intensity is Assume that the transverse shifting of the illuminating beam has a normal distribution , where D is determined by the standard deviation of ∆ = δ/λz. The recorded intensity with high-frequency positioning vibration is a summation of different intensities for different illumination shifts: The Fourier transform of the recorded diffraction intensity becomes Compared with Eq. (4), the position vibration degrades the contrast of the interference fringe in the recorded intensity, acting as a low-pass filter.
In practice, the pointing and transverse vibration of the illumination beam can occur simultaneously; thus, the Fourier transform of the recorded diffraction intensity is Because I 0 (u) has a narrow frequency band for spherical illumination, Eq. (12) can be approximated as where C = πλL √ K 2 + D 2 . Clearly, the instability of the illumination beam during the exposure of the detector leads to the loss of high-frequency components of the diffraction intensity, and the contrast of the interference fringes is differently suppressed depending on their frequency. As reported in the literature [22], change in the contrast of the diffraction patterns leads to mathematical ambiguity and generates errors in the final reconstruction.

Compensation method
At first glance, it appears that the Fourier component of the ideal diffraction pattern can be obtained by dividing exp(−C 2 u 2 ) to the I v (u) by each possible C until a satisfactory reconstruction is achieved. However, because this operation can seriously amplify the errors or noise, and the expired resolution improvement can be submerged by them, it cannot be adopted in practice.
In the proposed method, the recorded intensities I v (k) are modified before the conventional PIE process, as shown in Fig.1. The recorded diffraction intensities are deliberately blurred via convolution with the Gaussian function exp(−k 2 /K 2 v ), where K v is a constant chosen to slightly blur the recorded diffraction pattern. The deliberately blurred intensity pattern is expressed in the frequency domain as where C v = πλLK v . We subtract the blurred pattern I ′ v (k) from the original recorded diffraction I v (k) and then add the subtracted pattern to I v (k) after multiplying it by a constant β: The Fourier transform of the modified intensity pattern is Clearly, when (16) is wider than exp(−C 2 u 2 ), I c (u) is close to the Fourier transform of the vibration-free diffraction pattern, and then the influence of the imaging system instability on the reconstructed image can be remarkably suppressed. A numerical simulation is performed using the proposed algorithm to check its feasibility. A divergent spherical wave with a wavelength of 632.8 nm is used for the illumination. Two pictures are used as the modulus and phase, respectively, of the specimen. The diameter of the illuminated area on the specimen is 0.74 mm, and 10 × 10 diffraction intensities are recorded with a step size of 0.185 mm. The charge-coupled device (CCD) camera has a resolution of 256 × 256 pixels and a pixel size of 7.4 µm. For comparison, the diffractive intensities with stable illumination are calculated, one of which is shown in Fig.2(a). When the unstable imaging system results in diffraction pattern vibration with variance of 46µm in the detector plane, the intensity distributions are calculated by adding diffraction intensities under many slightly different illuminations with varying incident angles and transverse shifts. The diffraction patterns corresponding to the same illuminating position of Fig.2(a) is shown in Fig.2(b), where the intensity is obviously blurred. This coincides with Eq.(13). The reconstructed modulus and phase of the specimen with stable illumination are very clear as shown in Figs.2(d) and (g), respectively. Figs.2(e) and (h) show the reconstructed complex transmission of the specimen with the unstable illumination, where the resolution is apparently decreased compared with Figs.2(d) and (g). Corrected diffraction patterns are obtained from the raw recorded data using the proposed method. A Gaussian function is used to slightly blur the recorded intensities with K v = 1.24mm −1 . And the obtained blurred diffraction patterns are subtracted from the raw recorded diffraction patterns. Then the modified intensity patterns are obtained by adding the amplified subtractions to the raw data with β = 2.5. As shown in Fig.2(c), the contrast of the corrected diffraction pattern is obviously improved compared with the raw data and the distribution is similar to that of Fig.2(a). The reconstructed modulus and phase of the specimen using the corrected intensities are shown in Figs.2(f) and (i), respectively, where the resolution is remarkably improved compared with Figs.2(e) and (h). The insets in Figs.2(d)-(f) and (g)-(i) show the reconstructed modulus and phase, respectively, of the illumination field, indicating that the quality of the retrieved illumination is also improved using the proposed method.To quantify the performance of the proposed method, the normalized root-mean-square error metric [13] is calculated, as shown in Fig. 3. The accuracy of the reconstruction are obviously improved when the proposed method is used to correct the recorded diffraction patterns before the iterative computation. Fig. 3. Progress of the error of the conventional PIE method for a static imaging system (blue dashed line),an an unstable system (green dashed line); the proposed method for an unstable system(red line). The proof-of-principle experiment is conducted with visible light as shown in Fig. 4, where the divergent laser beam from a He-Ne laser is slightly diffused by a rotating plastic diffuser before irradiating on the sample to simulate the instability of the imaging system. A cross section of a monocotyledon placed on a 2D translation stage is used as the sample. The illumination area on the specimen is ∽ 2.5 mm, and 10 × 10 diffraction patterns are recorded by the CCD camera while the sample is scanned relative to the instable beam with a step size of 0.37 mm. Recorded intensities obtained with stable and unstable laser beams are shown in Figs.5(a) and (b), respectively. It is clear that the recorded intensities using unstable imaging system are totally blurred. In order to overcome this limitation, the newly proposed method is applied to modify the recorded patterns. The Gaussian function with K v = 0.9355mm −1 is adopted to slightly blur the recorded diffraction intensities. According to the proposed method, the corrected diffraction patterns are obtained with β = 5. As shown in Fig.5(c), the contrast of the corrected diffraction pattern is remarkably improved and similar to that shown in Fig.4(a) Here, the quality of the reconstructed images is distinctly improved, and the images obtained are roughly identical to those for the stable illumination.

Experimental result
To quantify the resolvability of the proposed method, a USAF1951 target is measured. As shown in Fig.6, group 5, which is resolvable for the stable illumination, is obviously blurred for the unstable system. After the proposed method is applied, the elements in Group 5 become distinguished, as shown in Fig.6(c) and (f).

The influence of β
To show the dependency of the resolution-improvement effect on the value of β, a numerical simulation is performed with the same parameters that were used for the simulation in Fig.2. Fig.7 shows the effect of the proposed method with various curves, where the blue dashed line indicates the low-pass filter related to the vibration of the imaging system, and the other curves show the increased width of the low-pass filter with varying β. The width of the low-pass filter increases remarkably with increasing β, and for β = 2.5, the width of is about twice that of the exp(−C 2 u 2 ), and accordingly the resolution of the final reconstruction with these generated diffraction patterns can be remarkably improved. However, when β = 5, the low-pass filter has the shape of the black dashed line in Fig.7; that is, the strength of the high-frequency components is over-amplified. Thus, the constant β should be carefully selected to assure a satisfactory final reconstruction. Fig.8 shows the reconstruction results with different β values according to the aforementioned analysis, where Fig.8 (a) and (b) are the reconstructed modulus and phase images, respectively, obtained using the raw data acquired with the unstable system. Fig.8(c) to (h) show the reconstructed modulus and phase images with β = 2.5, 1, and 5. These results indicate that the modulus and phase reconstructed using the proposed method produce a better reconstruction quality than those obtained directly with unstable data. On the other hand, the effect of the proposed method depends on the properly chosen value of β, which is coincident with the curves shown in Fig.7.

Robustness of the proposed method
In the aforementioned theoretical analysis, the vibration of the illumination is assumed to have a normal distribution, and another Gaussian function is used to convolve with the recorded intensities to slightly blur the patterns. But the characteristics of a real imaging system are unknown and may be more complex, so it appears that the proposed method is difficult to implement for real experiments. However, the fundamental principle of the proposed method is to recover the high-frequency components lost because of the instability of the radiation source, and other functions that can realize this purpose can replace the Gaussian function in the analysis. For example, circular functions, triangular functions, and para-curves can be adopted to slightly blur the diffraction patterns for improving the resolution without considering the exact properties of the vibration of the imaging system. This feature makes the proposed method more applicable to real experiments and is demonstrated by the following simulations and experiments. When the vibration of the imaging system follows a normal distribution, the parameters are the same as that were used for the simulations in Fig.2. As previously discussed, when the diffraction patterns are deliberately blurred via convolution with another Gaussian function exp(−k 2 /K 2 v ), the width of the low-pass filter shown by the red line in Fig.9(a) becomes much wider than that for the dashed line related to the instability. The same results can be obtained by Fig. 9. When the vibration follows normal distribution, filter window (a) modified with Gaussian function (K v = 1.24mm −1 , β = 2.5); (b) modified with circular function (K v = 1.1694mm −1 , β = 3.5); (c) modified with triangular function (K v = 0.8771mm −1 , β = 4); when the vibration follows random distribution, filter window (d) modified with Gaussian function (K v = 0.8268mm −1 , β = 2.6); (e) modified with circular function (K v = 0.8771mm −1 , β = 3.5); (f) modified with triangular function (K v = 0.7309mm −1 , β = 3.5); when the distribution of the vibration follows triangular function , filter window (g) modified with Gaussian function (K v = 0.8268mm −1 , β = 2); (h) modified with circular function (K v = 0.8771mm −1 , β = 3.2); (i) modified with triangular function (K v = 0.7309mm −1 , β = 2.4).
using the circular function circ(k/K v ) ( Fig.9(b)) and the triangular function Λ(k/K v ) ( Fig.9(c)). For the case where the vibration of the imaging system follows a random distribution with the maximum extent of 52.5µm in the detector plane, the low-pass filter window of the recorded intensity pattern is shown by the black dashed line in Fig.9(d). The modified results obtained by the deliberately blurring diffraction patterns via convolution with a Gaussian function, circular function, and triangular function are shown as the red lines in Fig.9(d) to (f), respectively. When the vibration of the imaging system has a triangular distribution and results in an extent of 36.2µm in the detector plane, the low-pass filter window of the recorded intensity pattern has the profile of the black dashed line shown in Fig.9(g). The modified results obtained by the deliberately blurring diffraction patterns via convolution with a Gaussian function, circular function, and triangular function are shown as the red lines in Fig.9(g) to (i), respectively.
And an experiment is carried out to verify the practicality of the proposed method, and the parameters are the same with Fig.5. The reconstructed images shown in Fig.10(a) and (b) are the reconstructions using the raw data recorded with instable illumination, which are seriously blurred. The reconstructed images shown in Fig.10 (c) and (d) are obtained by using the Gaussian function with K v = 0.9355mm −1 and β = 5 to do the resolution improvement with the proposed method. Fig.10 (e) and (f) are obtained by using intensities modified with circular function with a diameter K v = 1.4033mm −1 and β = 8 to do the reconstruction. Fig.10 (g) and (h) are obtained by using intensities modified with triangular function with a diameter K v = 1.1694mm −1 and β = 6 to do the reconstruction. We can find that their effects on the resolution improvement are almost the same. These results shows that the effect of the proposed method is independent of the vibrating model of the imaging system; thus, the method does not require prior knowledge of the vibrating properties.

Conclusion
In conclusion, a simple method is proposed to relax the requirement of ptychography for the stability of the imaging system. The influence of the instability of the imaging system can be described as a low-pass filter acting on the ideal diffraction intensities; the high-frequency components are lost in the recorded data, generating a blurred image in the final reconstruction. Using the proposed method, the recorded intensity patterns are corrected by subtracting the deliberately blurred diffraction patterns from the original ones and then adding the amplified subtracted patterns to the recorded data. The resolution of the final reconstruction can be significantly improved by using the corrected diffraction patterns. The feasibility of the proposed method is demonstrated via both numerical simulations and experiments on an optical bench. Because the proposed method does not require any prior knowledge of the characteristics of the illumination beam or massive calculation during the iterative process, it is an easy approach for acquiring a high-quality reconstruction with an unstable imaging system. The method can be extended to other CDI techniques for imaging with short-wavelength irradiation, such as free electrons or soft X-ray lasers.

Funding
National Natural Science Foundation of China (No. 61675215).