Open Access
Add to BibSonomy
Abstract
The single-shot capability of coherent modulation imaging (CMI) makes it have great potential in the investigation of dynamic processes. Its main disadvantage is the relatively low signal-to-noise ratio (SNR) which affects the spatial resolution and reconstruction accuracy. Here, we propose the improvement of a general spatiotemporal CMI method for imaging of dynamic processes. By making use of the redundant information in time-series reconstructions, the spatiotemporal CMI can achieve robust and fast reconstruction with higher SNR and spatial resolution. The method is validated by numerical simulations and optical experiments. We combine the CMI module with an optical microscope to achieve quantitative phase and amplitude reconstruction of dynamic biological processes. With the reconstructed complex field, we also demonstrate the 3D digital refocusing ability of the CMI microscope. With further development, we expect the spatiotemporal CMI method can be applied to study a range of dynamic phenomena.
© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement
1. Introduction
Coherent diffractive imaging (CDI) is a powerful lensless imaging technique with simple concepts and easy implementation. It collects diffraction intensities and then uses iterative feedback algorithms to reconstruct the sample exiting wavefields by solving the phase problem [1]. Since its first demonstration [2], CDI has gained rapid development and has been widely used in materials and biological sciences [3–8]. However, current CDI techniques still suffer from several shortcomings when investigating in a dynamic process. Conventional CDI methods are accompanied by complementary problems, plaguing with the uniqueness of the solutions and the stagnation of the algorithm [9–11]. Its scanning variant, ptychography [12–14], can faithfully reconstruct a complex object, but it is not suitable for dynamic samples due to the time-consuming data acquisition.
Coherent modulation imaging is a novel CDI method that can reconstruct a complex wavefield from a single measurement, thus has great potential for imaging dynamic processes [15–17]. It has been shown to address the drawbacks of CDIs: free from the inherent ambiguities, reducing the dynamic range, and enhancing the convergence. CMI has been demonstrated with visible light [16] and X-ray [17], and current researches mainly focus on the optimization of imaging systems and algorithms [18–24], as well as search for applications [25–27]. To reduce the speckle noise in CMI reconstructions, He et al [18] introduced a second detector to apply additional intensity constraints. Wang et al [23] proposed a modulator refinement algorithm to relax the requirement on exactly knowing the modulator. Various variations of conventional CMI, such as multi-wavelength CMI [19], continuous phase modulation CMI [20], amplitude modulation CMI [21], and cascade modulation CMI [22] have been reported. In addition, CMI has been applied to the diagnosis of laser beams [25,26] and the measurement of large optical elements [27].
Speckle noise is the main problem in CMI reconstruction. Although it can be improved by introducing a second detector, it undoubtedly increases the complexity of the CMI experiments. Recently, some methods in CDI have made use of the redundant information in a dynamic process itself to improve performance [28–30]. In time-series frames of a dynamic process, there usually is a time-variant region that changes over time and a time-invariant region that remains stationary. Therefore, these methods apply the time-invariant region as a powerful spatiotemporal constraint to reconstruct dynamic processes, achieving fast and robust convergence. The identification of the time-invariant region is achieved either by introducing a constant ’reference’ object [28] or restricting the object to a sub-region of the probe [29]. However, limitations are that the location of the reference object must be known and fixed, and the object must be confined to a region far smaller than the probe. Fortunately, a recently developed method has overcome the above limitations [30]. It uses an unsupervised approach to identify the time-invariant region without any changes to the original experimental arrangement, which permits the use of spatiotemporal constraints in other imaging methods for improving dynamic measurements.
Here, we report a method of using the spatiotemporal constraint in CMI, termed as spatiotemporal CMI, or ST-CMI. The spatiotemporal constraint is applied in the support plane. The time-invariant region is identified using the unsupervised method in [30] and constrained to equal values, while the time-variant region is updated normally during the iterative reconstruction process. We validate the spatiotemporal CMI method with numerical simulations and optical experiments. In the experiment, the CMI module is combined with an optical microscope to achieve quantitative phase and amplitude reconstruction of dynamic biological processes. With the reconstructed complex wavefield, we further demonstrate the 3D digital refocusing ability of our CMI microscope, which is very useful when imaging dynamic 3D processes.
The article is organized as follows: Section 2 describes the details of the spatiotemporal CMI method, Section 3 shows the simulation results to demonstrate the feasibility of the method, Section 4 describes the experimental setup and presents the experimental results, Section 5 concludes.
2. Method
2.1 Forward imaging model
Figure 1(a) shows a typical forward imaging process of CMI in optical near-field configuration. In a single-shot experiment, a complex-valued object is illuminated with coherent light, and then the exit wavefield $U_{0}(x,y)$ restricted by a pinhole $S$ propagates a distance $z_{1}$ to a modulator $M(x,y)$. The modulated wave $U_{2}(x,y)$ further propagates a distance $z_{2}$ to the detector $U_{3}(x,y)$, resulting in the diffraction pattern $I(x,y)$ that is recorded. The whole forward imaging process of CMI can be summarized as follow:
where $(x,y)$ denotes the spatial coordinate, $\mathcal {P}$ denotes the wave propagation operator, $z1$ is the support-to-modulator distance and $z2$ is the modulator-to-detector distance, and $\cdot$ denotes the element-wise multiplication. For simplification, the spatial coordinate $(x,y)$ is omitted in the rest of this article.
Fig. 1. (a) Schematic of the CMI set-up. (b) Flowchart of the CMI phase retrieval algorithm; the arrows indicate wave-propagation steps.
In the CMI reconstruction algorithm, one key prior information is the transmission function of the modulator $M$, which can be calibrated using ptychography. With the calibrated modulator function and recorded diffraction pattern, the CMI algorithm can be used to reconstruct the exit wavefield of the complex-valued object. Flowchart of the CMI reconstruction algorithm is shown in Fig. 1(b). It starts with an initial guess of the wavefield at the support plane $U_{0}$, which is known to have non-zero values within a finite region. Reconstruction is achieved by iteratively propagating the wavefield estimate $U_{0}^{k}$ between the support plane and the detector plane, via the modulator. The use of a modulator in CMI can reduce the possible solutions that backpropagate to the support plane, yielding rapid convergence without twin-image and spatial shift ambiguities. For a full explanation of the conventional CMI algorithm, the reader can refer to [17].
2.2 Spatiotemporal CMI algorithm
In dynamic CMI experiments, time series frames of diffraction patterns $I_{n}$ are recorded with a fast-framing detector. Conventional CMI algorithm treats each diffraction pattern $I_{n}$ independently to reconstruct $U_{0}^{n}$. One can see from the reconstructions that consecutive frames have a high similarity: the time-variant region changes over time while the time-invariant region remains unchanged. Recently developed methods have made use of this redundant information in time-series frames, termed as spatiotemporal constraints, to improve the performance of traditional CDI [28–30]. Spatiotemporal CMI follows the same idea. In the spatiotemporal CMI reconstruction algorithm, the spatiotemporal constraint is applied in the support plane, the time-invariant region is identified and constrained to equal values, while the time-variant region is updated normally.
The pseudocode of the reconstruction algorithm is depicted in Algorithm 1, where $I_{n}$ denotes the $n^{th}$ diffraction pattern of the $N$ recorded frames; the running estimate at $k^{th}$ iteration for the $n^{th}$ object is denoted as $U_{0}^{n,k}$, $U_{1}^{n,k}$, $U_{2}^{n,k}$ and $U_{3}^{n,k}$ for exit wave at support plane, the incident and exit wave of the modulator, and the wave at the detector plane; $\mathcal {P}_{z}$ and $\mathcal {P}_{z}^{-1}$ are the forward and backward propagation operators for distance $z$, respectively; corresponding revised wave at each plane is denoted as $u_{0}^{n,k}$, $u_{1}^{n,k}$, $u_{2}^{n,k}$ and $u_{3}^{n,k}$; $G_{k}$ denotes the time-average object estimate at $k^{th}$ iteration; $\sigma (x,y)$ denotes the standard deviation of the phase or amplitude of the revised object exit wave and $D$ denotes the time-invariant region updated according to $\sigma (x,y)$. The superscript $*$ denotes the complex conjugation.
The reconstruction algorithm starts with random initial guesses of object exit waves $\{U_{0}^{n,0}\}$ that are constrained by support $S$. At the beginning of each iteration loop, the time-averaged object estimate $G_{k}$ is obtained by calculating the average of all the updated exiting wave $u_{0}^{n,k-1}$ at the previous iteration (line 2). Then calculate the standard deviation $\sigma (x,y)$ from the phase or amplitude of $u_{0}^{n,k-1}$ over all values of $n$ (line 3). The value of $\sigma (x,y)$ within the time-invariant region is relatively small and close to zero, which permits us to identify the time-invariant region $(D=1)$ using some binarization method (line 4). The binarization threshold $\beta$ can be determined manually or by using automatic segmentation approaches. For each measurement, the spatiotemporal constraint is applied by updating $U_{0}^{n,k}$ within the areas where $D=0$, while the values within the time-invariant region are averaged (line 6). Here the term $(1-\gamma )u_{0}^{n,k-1}+\gamma G_{k}$ is a weighted update function within the identified time-invariant region and $\gamma$ is a tuning parameter that weights the amount of information shared by $G_{k}$. When $\gamma = 0$, it becomes the conventional CMI algorithm. The support plane wave $U_{0}^{n,k}$ is then propagated to the modulator plane (line 7) and interacted with the function of the modulator $M$ (line 8). The modulated wave $U_{2}^{n,k}$ is then propagated to the detector plane (line 9) where modulus constraint is applied, yielding the corrected wave $u_{3}^{n,k}$ (line 10); $\epsilon$ is a small number used here to avoid division by zero. The corrected wave is then back-propagated to the modulator plane (line 11) to remove modulation (line 12), the resulted wave $u_{1}^{n,k}$ is further back-propagated to the support plane and constrained by the support $S$ (line 13), obtaining the revised support plane wave $u_{0}^{n,k}$. The algorithm continues to be processed with a predetermined iteration number $K$.
In spatiotemporal CMI reconstruction, an important step is to correctly identify the time-variant and time-invariant regions. We find that smoothing over the standard deviation $\sigma (x,y)$ with average filter helps to reduce the noise, yielding better identification of time-invariant region $D$. The wave propagation operator $\mathcal {P}_{z}$ can be angular spectrum propagation or Fourier transform, depending on the experimental setup. For optical near filed configuration, angular spectrum propagation is used in the simulations and experiments below.
3. Simulations
3.1 Simulated imaging of dynamic object
We first demonstrate the spatiotemporal CMI algorithm with numerical simulations, as we show in Code 1 (Ref. [34]). The diffraction patterns are generated according to the setup in Fig. 1 and the images used in the simulation are $512\times 512$ pixels. The complex object is generated from a ’Hela cells’ image and its inversion, respectively defining the amplitude and phase of the object. The support region is $380$ pixels in diameter and is illuminated with light of wavelength $\lambda = 632.8$ nm. The modulator is placed $z_{1} = 20$ mm downstream of the object. Its amplitude is assumed to be arrays of 1 and its phase distribution is random values of 0 or $\pi$. The inset on the top-right in Fig. 2(a) is a zoom-in of the central $40\times 40$ pixels region. The detector is placed $z_{2}=10$ mm further downstream of the modulator. The pixel size is $2.74$ $\mathrm{\mu}$m. To simulate a dynamic biological process, we rotate the simulated object one degree at a time and a total of 30 frames of diffraction patterns are generated, one of them is shown in Fig. 2(d).
Fig. 2. The images used in the simulation. (a) The random phase modulator, inset is the room-in of the central $40 \times 40$ pixels region. (b) Amplitude and (c) phase of the simulated object wave. (d) One of the generated diffraction patterns.
Then we use the spatiotemporal CMI algorithm to reconstruct the dynamic process from the generated diffraction patterns. The tuning parameter $\gamma$ is set to 1 and the threshold $\beta$ is set to 0.035, manually. For comparison, we also perform the conventional CMI algorithm, which can be achieved by simply setting $\gamma = 0$. We run 100 iterations of the algorithms, as more iterations do not bring significant improvement. The reconstructed results of Frame = 1, 15, and 30 are depicted in Figs. 3(a) and (b). The final standard deviation $\sigma (x,y)$ and the binary map of it are shown in Figs. 3(c) and (d), respectively.
Fig. 3. Reconstruction results of conventional CMI (a1,a2,a3) and spatiotemporal CMI (b1,b2,b3). (c) The standard deviation in reconstructed amplitude and (d) the binary-valued map of it. (e) Convergence curves of the algorithms. (f) PSNR of reconstruction results. (g) Solid line, Fourier ring correlation (FRC) between the reconstructed object waves and the ground truths. Dashed line, the 1/2 bit resolution curve.
As shown in Figs. 3(a1) and (b1), the central $200\times 200$ pixels within the red dashed box of the reconstructed object are extracted for quantitative analysis. To measure the convergence of each algorithm, we use the Normalised Root Mean Square Error (NRMSE) as the error metric, as shown in Fig. 3(e). The NRMSE values are calculated between the reconstructed wavefields $\{U_{0}^{n,k}\}$ and the ground truth $\{U_{0}^{n}\}$:
We observe that the noise in reconstruction with spatiotemporal constraint seems to be reduced, thus we calculate the Peak to Signal Noise Ratio (PSNR) between the reconstructed amplitude $\{|U_{0}^{n,K}|\}$ and the corresponding ground truths $\{|U_{0}^{n}|\}$:
where $MSE$ is the mean square error, $MaxPV$ is the maximum of the pixel value, the size of $U_{0}^{n,k}$ and $U_{0}^{n}$ are $P\cdot Q$. Figure 3(f) shows the results, where the PSNR using ST-CMI (38.3dB on average) is generally higher than that using conventional CMI (26.4dB on average). The spatial resolution is evaluated by Fourier ring correlation (FRC) [31] between the reconstructed object waves and the ground truths, and the 1/2 bit resolution criteria are used. The correlation curve of ST-CMI is above the 1/2 bit threshold across the detected spatial frequency spectrum, while the curve of CMI drops to the threshold at about 0.7 of the Nyquist frequency. Thus, the reconstruction using ST-CMI can achieve higher spatial resolution, because the higher spatial frequencies are affected by the noise in CMI results.It should be noted that the spatiotemporal constraints in fact bring improvement to the reconstruction of the whole field of view, including both the variant and invariant region. We have also conducted quantifications within different regions, for example, using only the pixels in the variant region; and the results are basically the same as the ones described above, as shown in Fig. 4.
Fig. 4. (a) Using the pixels within the variant region for quantitative analysis. (b) Convergence curves of the algorithms. (c) PSNR of reconstruction results.
3.2 Analysis of the tuning parameter $\gamma$
To evaluate the influence of the tuning parameter $\gamma$ on the performance of the spatiotemporal CMI algorithm, we apply the reconstruction in the case of $\gamma =0,0.2,0.6,1$, keeping other parameters unchanged. The convergence curve within 100 iterations and the final PSNR values of the reconstruction results are shown in Figs. 5(a) and (b), respectively. From the convergence curve, we can see that in the first 40 iterations, the error converges faster with the increase of $\gamma$. With the increase of the number of iterations, the errors converge to the same order of magnitude except when $\gamma =0$, which indicates that a smaller $\gamma$ can also help the convergence. Moreover, the PSNR values of the reconstruction results when $\gamma \neq 0$ is generally higher than that when $\gamma =0$, which proves the superiority of the spatiotemporal CMI method.
3.3 Analysis of increased dynamics
As our method is based on spatiotemporal constraints, it is meaningful to study whether the ratio of time-invariant/variant regions influences the reconstruction quality. To do this, we change the step size of each rotation angle $\Delta \phi$ of the object while keeping other simulation parameters unchanged as described in Section 3.1. We simulate 25 groups of dynamic diffraction patterns, and the total rotation angle $\phi$ from the last frame to the first frame increases from 0 degrees to 120 degrees with an interval of 5 degrees. The spatiotemporal CMI method (setting $\gamma =1$) is used to reconstruct each group. The reconstructed object amplitudes of the first frame when $\phi =5^{\circ }, 40^{\circ }, 75^{\circ }, 110^{\circ }$ are given in Fig. 6(a), and the corresponding identification results of time-invariant/variant region are shown in Fig. 6(b). The percentages in the lower right of Fig. 6(b) refer to the proportion of the identified time-invariant region to the field of view (FOV). To quantify the reconstruction quality in these cases, similarly, we select the central 200 $\times$ 200 pixels of the reconstruction amplitude of the first frame in each group to calculate the PSNR, and the results are shown in Fig. 6(c).
Fig. 6. (a) Reconstructed amplitude of the first frame in each group using the spatiotemporal CMI method. (b) Identification results of time-invariant/variant region, which are binary maps of the corresponding standard deviation in reconstructed amplitude. (c) PSNR results of the first frame with respect to the increased total rotation angle $\phi$.
As shown in Fig. 6, the reconstruction quality is improved when $\phi = 5^{\circ }$ compared to the conventional one ($\phi = 0$), although the identification is not completely correct. The spatiotemporal constraints perform best when $\phi =15^{\circ }$, the corresponding ratio of the time-invariant region to the FOV is $61.88\%$. As the dynamics increase, the reconstruction quality decreases slowly, because the time-invariant area where spatiotemporal constraints can be imposed decreases. When the dynamics increase to such an extent that spatiotemporal constraints cannot be imposed on any region, we expect the proposed method degrades to the conventional one. In practical issues, we can control the number of input frames $N$ each time to adjust the proportion of the time-invariant/variant region, so as to achieve high-quality reconstruction as much as possible.
4. Experiments
4.1 Set-up
In order to observe dynamic biological processes on various scales, we combine the CMI module with a lab-assembled optical microscope. As illustrated in Fig. 7(a), a sample is illuminated by a coherent collimated beam with a wavelength of 637 nm (COHERENT, OBIS 637LX); an interchangeable objective forms a magnified image of the sample; a 1mm diameter circular aperture is placed at the image plane to crop the wave; the modulator is mounted on a motorized mechanical x-y stage (DAHENG, GCD0401M) and placed downstream of the aperture; the diffraction detector (BASLER, ace2, a2A4504-18umPRO) with $4512\times 4512$ pixels and a 2.74 $\mathrm{\mu}$m pixel size is placed downstream of the modulator. Note that the pinhole is placed at the image plane instead of the sample plane. Placing at the sample plane can bring extra experimental difficulties to locate the living sample as the pinhole have to be much smaller because of the magnification.
Fig. 7. (a) Experimental setup of the CMI microscope, inside the red box is the CMI module. (b) A typical experimental measured diffraction pattern; the calibrated (c) amplitude and (d) phase of the random phase modulator using ptychography, the inset of (c) is the recovered probe amplitude. (e) The intensity distribution in the support plane, the red dashed circle indicates the support constraint.
In our experiment, the modulator is made of silica glass with a photoetching technique; different etching depths on the substrate brings different phase delay. It roughly has a $0-\pi$ binary structure for the employed wavelength and is designed to be randomly distributed. The physical size of each pixel is $16$ $\mathrm{\mu}$m $\times$ $16$ $\mathrm{\mu}$m. It is worth mentioning that CMI has no strict requirement on the distribution of the modulator, as continuous phase modulation [20] and amplitude modulation [21] are also reported. The exact transmission function of the modulator is calibrated using ptychography [14] with the ePIE algorithm [13], and the reconstructed amplitude and phase are shown in Figs. 7(c) and (d), respectively. The inset in Fig. 7(c) is the recovered probe, backpropagate which can obtain the wavefield at the support plane, as shown in Fig. 7(e). The propagation distances $z1$ and $z2$ can also be determined during calibration. The calibration process is detailed in Appendix A.
Before recording the CMI dataset, the focusing detector placed on a separate path is used to focus the magnified image to the correct plane. In continuous acquisition mode, time-series diffraction patterns of a dynamic process can be recorded, an example from the system is shown in Fig. 7(b).
4.2 Dynamic biological sample imaging
To demonstrate the efficacy of our CMI microscope, investigations of live Paramecium are carried out in this section. Using the tension of water, the living Paramecium sample is sandwiched between two cover slides with a thin steel sheet. For this test, we fitted a $10\times$ objective lens with NA=0.25, and the measured magnification at the plane of the diffraction detector is $13.15\times$. In continuous acquisition mode, we collected 1800 diffraction patterns of $1024\times 1024$ pixels at 30 frames per second in one minute. In the process of collection, the Paramecium often moved out of the field of view, resulting in a lot of empty data (diffraction patterns without the sample). We thus select three groups of the recorded continuous diffraction patterns for reconstruction.
The diffraction patterns are processed using the spatiotemporal CMI algorithm, with the tuning parameter $\gamma =1$ and the threshold $\beta =0.035$. The standard deviation $\sigma$ is calculated from the reconstructed intensities and the iteration number $K$ is set to 50. A loose support constraint of 500 pixels in diameter is used in the reconstruction, as shown in the red dashed circle in Fig. 7(e). The corrected propagation distances $z1$ and $z2$ are $30$ mm and $26$ mm, respectively. Three groups of data are processed separately. A diffraction pattern without sample is added to each group to reconstruct the illumination wavefield. Note that the reconstructed wavefield at the support plane may be defocused. Fortunately, we can propagate the reconstructed wavefield to any axial planes to obtain the in-focus images, which is detailed in the next section. After refocusing, subtraction of the illumination wavefield from the in-focus exit wavefield yields the final sample wavefield shown in Fig. 8. These results are also visualized in Visualization 1 and Visualization 2.
Fig. 8. Experimental results of dynamic biological sample. The recovered (a1,a2,a3) amplitude and (b1,b2,b3) phase of the selected frames. Corresponding (c1,c2,c3) standard deviation in reconstructed amplitude and (d1,d2,d3) the binary-valued map of it. See Visualization 1 and Visualization 2 for the reconstructed video.
Figures 8(a) and (b) show the reconstructed amplitude and phase of the sample wavefield, respectively. It can be seen that the movement of young Paramecium (the smaller one) is better observed in the phase images. Figures 8(c) and (d) show the final standard deviation and the corresponding binary maps. For comparison, we also perform the conventional CMI algorithm (setting $\gamma =0$) on the experimental data. The reconstructed phase of one frame using the two methods is compared as shown in Fig. 9. To quantitatively investigate the improvement of the SNR, we use the relative standard deviation (RSD), which is calculated as follows:
where $x_{i}$ is the pixel value, $n$ is the total pixel number and $\bar {x}$ is the corresponding average value. Here, we calculate the RSD using the pixels within the white solid box in Fig. 9. The RSD with conventional CMI is $3.34\%$, while the RSD with spatiotemporal CMI is $2.24\%$, which is about $67\%$ of the conventional CMI. The improvement of the SNR could be visually observed according to the zoomed-in color-mapped images in Fig. 9.Fig. 9. Reconstructed phase of one frame using the spatiotemporal CMI algorithm and the conventional CMI algorithm. The zoomed-in images are in HSV color scales.
4.3 3D digital refocusing
In Fig. 10, we demonstrate the 3D digital refocusing ability of the CMI microscope. We first test it with a static bee wing sample, using the same experimental setup as before. Figures 10(a1)-(a2) show the recovered complex wavefield at the support plane. Note that for the static sample we can only use the conventional CMI reconstruction algorithm (setting $\gamma = 0$), thus speckle noise is obvious in the images. Similar to many holographic imaging techniques, we can propagate the recovered complex wavefield to any plane along the optical axis. Figures 10(a3)-(a4) show the intensity images at the plane of $z=-1$ mm and $z=1$ mm, respectively. Note that we are using the angular spectrum propagation method at these working distances. We can see that different parts of the sample are brought into focus in different axial positions, as shown in the blue dash box in Figs. 10(a3)-(a4).
Fig. 10. Demonstration of 3D refocusing using the CMI microscope scheme. The recovered (a1) amplitude and (a2) phase of a static bee wing sample at support plane; the recovered section at (a3) z = -1 mm and (a4) z = 1 mm. The recovered (b1) amplitude and (b2) phase of one frame in the dynamic movement of Paramecium at support plane; the recovered section at (b3) z = -25 mm and (b4) z = -11 mm. Inside the blue boxs is the in-focus region.
The digital refocusing ability is very useful for imaging dynamic 3D samples. Take the Paramecium, for example. Paramecium’s motion is three-dimensional, and it can move along the optical axis, so we can’t guarantee that the image plane is exactly on the plane of the circular aperture. Thanks to the 3D refocusing ability, we can digitally refocus the recovered defocused images, as shown in Figs. 10(b1)-(b4). Figures 10(b1)-(b2) show the recovered amplitude and phase at the support plane using the spatiotemporal CMI algorithm, which are obviously defocused. We then propagate the defocused wavefield to different planes along the axis using the angular spectrum propagation method. The intensity images at the plane of $z = -25$ mm and $z = -11$ mm is shown in Figs. 10(b3)-(b4), where the in-focus part is highlighted in the dashed blue circle.
5. Conclusion
In this article, we propose a spatiotemporal CMI method for imaging dynamic processes. The method is validated through numerical simulations and optical experiments. In the simulations, comparisons between the spatiotemporal CMI and the conventional CMI are given and the results are quantified with PSNR and FRC; the tunning parameter $\gamma$ and the effect of increased dynamics are also analyzed. In the experiments, we build a CMI microscope and present experimental reconstruction results of dynamic biological processes; the 3D digital focusing ability of the CMI microscope is also demonstrated. The main contributions are summarized as follows:
Firstly, we demonstrate that the spatiotemporal constraint can be applied to dynamic CMI reconstruction. Through numerical simulations and optical experiments, we have shown that the spatiotemporal CMI method can reconstruct complex wavefields with fast convergence and higher quality.
Secondly, the implementation of our CMI microscope is simple, as long as the modulator is well calibrated. The CMI module can be acted as an ’add-on’ to standard microscopes, which provides an alternative solution to dynamic quantitative phase imaging. Moreover, we have demonstrated the 3D digital refocusing ability of CMI, which is of great benefit for the imaging of 3D dynamic processes.
Although the spatiotemporal CMI method is demonstrated for dynamic biological samples using visible light, the approach is in principle applicable to a broad range of wavelengths and radiations, such as X-rays. Possible applications of our method include the imaging of various dynamic processes, such as live cell monitoring, imaging of crystal formation, multi-particle tracking, and so on. With further development, we expect this general spatiotemporal CMI method could pave the way for its wide use in the study of fast sample dynamics.
Appendix A: Modulator calibration
The modulator transmission function is measured by ptychography. A ptychographic dataset of 49 diffraction patterns over an area of approximately 2 mm $\times$ 2 mm is recorded as the modulator is scanned over a 7 $\times$ 7 regular grid with random offsets added to avoid the raster grid pathology artifact [32]. The sharpness-statistics-based auto-focusing algorithm [33] is used to calibrate the modulator-detector distance $z2$ and then the extended ptychographical iterative engine algorithm [13] is run on the ptychographic dataset for 100 iterations, yielding the results shown in Figs. 7(c) and (d). The recovered probe is then backpropagated to find the distance $z1$ when the pinhole boundary is sharpest, which is the calibrated pinhole-modulator distance.
Funding
National Natural Science Foundation of China (61575197, 61975205, 6207522); Fusion Foundation of Research and Education of CAS; Youth Innovation Promotion Association of the Chinese Academy of Sciences (2017489); University of Chinese Academy of Sciences; Fundamental Research Funds for the Central Universities; Natural Science Foundation of Hebei Province (F2018402285); Funded Project of Hebei Province Innovation Capability Improvement Plan, China (No. 20540302D).
Acknowledgments
Junhao Zhang thanks Dr. Wenhui Xu (Southern University of Science and Technology, China) for the valuable and helpful discussions in the early stage of this work.
Disclosures
The authors declare no conflicts of interest.
Data availability
The Matlab simulation code is made available in the Supplementary Material, Code 1 (Ref. [34]).
References
1. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758–2769 (1982). [CrossRef]
2. J. Miao, P. Charalambous, J. Kirz, and D. Sayre, “Extending the methodology of X-ray crystallography to allow imaging of micrometre-sized non-crystalline specimens,” Nature 400(6742), 342–344 (1999). [CrossRef]
3. J. M. Zuo, “Atomic Resolution Imaging of a Carbon Nanotube from Diffraction Intensities,” Science 300(5624), 1419–1421 (2003). [CrossRef]
4. A. Barty, S. Marchesini, H. N. Chapman, C. Cui, M. R. Howells, D. A. Shapiro, A. M. Minor, J. C. H. Spence, U. Weierstall, J. Ilavsky, A. Noy, S. P. Hau-Riege, A. B. Artyukhin, T. Baumann, T. Willey, J. Stolken, T. van Buuren, and J. H. Kinney, “Three-Dimensional Coherent X-Ray Diffraction Imaging of a Ceramic Nanofoam: Determination of Structural Deformation Mechanisms,” Phys. Rev. Lett. 101(5), 055501 (2008). [CrossRef]
5. D. Shapiro, P. Thibault, T. Beetz, V. Elser, M. Howells, C. Jacobsen, J. Kirz, E. Lima, H. Miao, A. M. Neiman, and D. Sayre, “Biological imaging by soft x-ray diffraction microscopy,” Proc. Natl. Acad. Sci. 102(43), 15343–15346 (2005). [CrossRef]
6. M. A. Pfeifer, G. J. Williams, I. A. Vartanyants, R. Harder, and I. K. Robinson, “Three-dimensional mapping of a deformation field inside a nanocrystal,” Nature 442(7098), 63–66 (2006). [CrossRef]
7. H. Jiang, C. Song, C.-C. Chen, R. Xu, K. S. Raines, B. P. Fahimian, C.-H. Lu, T.-K. Lee, A. Nakashima, J. Urano, T. Ishikawa, F. Tamanoi, and J. Miao, “Quantitative 3D imaging of whole, unstained cells by using X-ray diffraction microscopy,” Proc. Natl. Acad. Sci. 107(25), 11234–11239 (2010). [CrossRef]
8. T. Kimura, Y. Joti, A. Shibuya, C. Song, S. Kim, K. Tono, M. Yabashi, M. Tamakoshi, T. Moriya, T. Oshima, T. Ishikawa, Y. Bessho, and Y. Nishino, “Imaging live cell in micro-liquid enclosure by X-ray laser diffraction,” Nat. Commun. 5(1), 3052 (2014). [CrossRef]
9. J. R. Fienup and C. C. Wackerman, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3(11), 1897–1907 (1986). [CrossRef]
10. M. Guizar-Sicairos, A. Diaz, M. Holler, M. S. Lucas, A. Menzel, R. A. Wepf, and O. Bunk, “Phase tomography from x-ray coherent diffractive imaging projections,” Opt. Express 19(22), 21345–21357 (2011). [CrossRef]
11. M. Guizar-Sicairos and J. R. Fienup, “Understanding the twin-image problem in phase retrieval,” J. Opt. Soc. Am. A 29(11), 2367–2375 (2012). [CrossRef]
12. J. M. Rodenburg and H. M. L. Faulkner, “A phase retrieval algorithm for shifting illumination,” Appl. Phys. Lett. 85(20), 4795–4797 (2004). [CrossRef]
13. A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109(10), 1256–1262 (2009). [CrossRef]
14. F. Zhang, I. Peterson, J. Vila-Comamala, A. Diaz, F. Berenguer, R. Bean, B. Chen, A. Menzel, I. K. Robinson, and J. M. Rodenburg, “Translation position determination in ptychographic coherent diffraction imaging,” Opt. Express 21(11), 13592–13606 (2013). [CrossRef]
15. F. Zhang, G. Pedrini, and W. Osten, “Phase retrieval of arbitrary complex-valued fields through aperture-plane modulation,” Phys. Rev. A 75(4), 043805 (2007). [CrossRef]
16. F. Zhang and J. M. Rodenburg, “Phase retrieval based on wave-front relay and modulation,” Phys. Rev. B 82(12), 121104 (2010). [CrossRef]
17. F. Zhang, B. Chen, G. R. Morrison, J. Vila-Comamala, M. Guizar-Sicairos, and I. K. Robinson, “Phase retrieval by coherent modulation imaging,” Nat. Commun. 7(1), 13367 (2016). [CrossRef]
18. X. He, H. Tao, X. Pan, C. Liu, and J. Zhu, “High-quality laser beam diagnostics using modified coherent phase modulation imaging,” Opt. Express 26(5), 6239–6248 (2018). [CrossRef]
19. X. Dong, X. Pan, C. Liu, and J. Zhu, “Single shot multi-wavelength phase retrieval with coherent modulation imaging,” Opt. Lett. 43(8), 1762–1765 (2018). [CrossRef]
20. X. Pan, C. Liu, and J. Zhu, “Phase retrieval with extended field of view based on continuous phase modulation,” Ultramicroscopy 204, 10–17 (2019). [CrossRef]
21. X. Pan, C. Liu, and J. Zhu, “Coherent amplitude modulation imaging based on partially saturated diffraction pattern,” Opt. Express 26(17), 21929–21938 (2018). [CrossRef]
22. Z. He, B. Wang, J. Bai, G. Barbastathis, and F. Zhang, “High-quality reconstruction of coherent modulation imaging using weak cascade modulators,” Ultramicroscopy 214, 112990 (2020). [CrossRef]
23. B. Wang, Q. Wang, W. Lyu, and F. Zhang, “Modulator refinement algorithm for coherent modulation imaging,” Ultramicroscopy 216, 113034 (2020). [CrossRef]
24. I. Kang, F. Zhang, and G. Barbastathis, “Phase extraction neural network (PhENN) with coherent modulation imaging (CMI) for phase retrieval at low photon counts,” Opt. Express 28(15), 21578–21600 (2020). [CrossRef]
25. X. He, S. P. Veetil, C. Liu, S. Gao, Y. Wang, J. Wang, and J. Zhu, “Accurate focal spot diagnostics based on a single shot coherent modulation imaging,” Laser Phys. Lett. 12(1), 015005 (2015). [CrossRef]
26. X. Pan, S. P. Veetil, C. Liu, H. Tao, Y. Jiang, Q. Lin, X. Li, and J. Zhu, “On-shot laser beam diagnostics for high-power laser facility with phase modulation imaging,” Laser Phys. Lett. 13(5), 055001 (2016). [CrossRef]
27. H. Tao, S. P. Veetil, J. Cheng, X. Pan, H. Wang, C. Liu, and J. Zhu, “Measurement of the complex transmittance of large optical elements with modulation coherent imaging,” Appl. Opt. 54(7), 1776–1781 (2015). [CrossRef]
28. X. Tao, Z. Xu, H. Liu, C. Wang, Z. Xing, Y. Wang, and R. Tai, “Spatially correlated coherent diffractive imaging method,” Appl. Opt. 57(22), 6527–6533 (2018). [CrossRef]
29. Y. H. Lo, L. Zhao, M. Gallagher-Jones, A. Rana, J. J. Lodico, W. Xiao, B. C. Regan, and J. Miao, “In situ coherent diffractive imaging,” Nat. Commun. 9(1), 1826 (2018). [CrossRef]
30. G. N. Hinsley, C. M. Kewish, and G. A. van Riessen, “Dynamic coherent diffractive imaging using unsupervised identification of spatiotemporal constraints,” Opt. Express 28(24), 36862–36872 (2020). [CrossRef]
31. M. van Heel and M. Schatz, “Fourier shell correlation threshold criteria,” J. Struct. Biol. 151(3), 250–262 (2005). [CrossRef]
32. P. Thibault, M. Dierolf, O. Bunk, A. Menzel, and F. Pfeiffer, “Probe retrieval in ptychographic coherent diffractive imaging,” Ultramicroscopy 109(4), 338–343 (2009). [CrossRef]
33. R. Ma, D. Yang, T. Yu, T. Li, X. Sun, Y. Zhu, N. Yang, H. Wang, and Y. Shi, “Sharpness-statistics-based auto-focusing algorithm for optical ptychography,” Opt. Lasers Eng. 128, 106053 (2020). [CrossRef]
34. J. Zhang, D. Yang, Y. Tao, Y. Zhu, W. Lv, D. Miao, C. Ke, H. Wang, and Y. Shi, “Supplemental code for spatiotemporal coherent modulation imaging,” figshare (2021), https://doi.org/10.6084/m9.figshare.14832810.
