Open Access
Add to BibSonomy
Abstract
Coherent modulation imaging is a lensless imaging technique, where a complex-valued image can be recovered from a single diffraction pattern using the iterative algorithm. Although mostly applied in two dimensions, it can be tomographically combined to produce three-dimensional (3D) images. Here we present a 3D reconstruction procedure for the sample’s phase and intensity from coherent modulation imaging measurements. Pre-processing methods to remove illumination probe, inherent ambiguities in phase reconstruction results, and intensity fluctuation are given. With the projections extracted by our method, standard tomographic reconstruction frameworks can be used to recover accurate quantitative 3D phase and intensity images. Numerical simulations and optical experiments validate our method.
© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement
1. Introduction
Three-dimensional (3D) imaging techniques are highly desirable in both biology and materials science applications. Typical 3D imaging techniques provide 3D information about the absorption or reflection property of the sample, such as computerized tomography (CT) [1], 3D laser scanning [2], optical coherent tomography (OCT) [3], and Magnetic Resonance Imaging (MRI) [4]. In addition, the information about the sample is also available as a phase shift of the incident wave, which provides higher contrast images for phase objects like transparent biological samples, making 3D phase imaging an attractive technique. 3D phase imaging generally involves tomography, which requires two steps: first, the 2D phase is computed at each angle, and second, the data are used as input for tomographic reconstruction to obtain the final 3D image [5]. 2D phase imaging methods, like holography [6–8], transport of intensity equation (TIE) [9,10], and coherent diffractive imaging (CDI) [11–19], are easily combined with tomography. However, these methods usually require multiple measurements at each angle, or a reference beam.
A good case in point is CDI, which uses diffraction patterns to reconstruct the wavefield distribution through iterative algorithms [20]. Nevertheless, the conventional iterative methods are usually time-consuming, and plaguing with the uniqueness of the solutions and the stagnation of the algorithm [21–23]. Ptychographic CDI [24–28] can faithfully reconstruct a complex wavefield from multiple partially overlapped diffraction patterns, but the data collection is time-consuming, as usually tens to hundreds of diffraction patterns need to be recorded.
Another feasible route to recover phase quantitatively is coherent modulation imaging (CMI), which is a novel variant of single-shot CDI [29–31]. Compared to the conventional CDI methods, CMI introduces an additional phase modulation constraint to eliminate ambiguous solutions during the iterative process, leading to faster convergence. CMI can recover the complex wavefield from one single measurement, thus having great potential for 3D phase imaging. It has been successfully demonstrated in 2D phase reconstruction with visible light [30] and X-rays [31], and current researches mainly focus on the optimization of imaging systems and algorithms [32–42], as well as the search for applications [43–45]. Recently, we showed that spatiotemporal constraints can be applied to CMI reconstruction, which can make 2D dynamic measurements more feasible [40].
Here we report a 3D phase imaging method based on CMI, termed 3D CMI, in which only one diffraction pattern is recorded at each angle. Because the reconstructions are complex-valued, 3D CMI recovers both 3D phase and intensity images, providing the sample’s 3D absorption and phase-shift properties simultaneously. However, to obtain accurate 3D results, several processing steps need to be carried out before tomographic reconstruction. First, what CMI recovers is the exit wavefield of the sample, therefore the illumination probe should be removed to obtain the transmissivity of the sample. Second, CMI reconstructions also suffer from additive constant and linear phase terms, as well as phase wrap. And last, fluctuation in the recovered intensity can be present if the light source is unstable. In this work, we present a projection extraction procedure to deal with the problems mentioned above. In our method, the removal of the illumination probe relies on measuring a diffraction pattern without the sample, from which the illumination function can be directly reconstructed. For the removal of constant and linear phase terms we rely on measuring a region of free space outside the sample, thus knowing that for this region the phase should be zero, providing no phase shift to the incident wave. The removal of intensity fluctuation (power correction) also benefits from the region of free space, because the intensity of this region should be consistent. After the removal of all these terms, accurate phase and intensity projections are obtained for further tomographic reconstruction. We carried out numerical simulations and proof-of-principle optical experiments to verify the feasibility of our proposed method.
The paper is organized as follows: Section 2 describes the forward model of the 3D reconstruction problem with CMI measurements. Section 3 introduces the phase retrieval algorithm of CMI. In Section 4, the projection extraction method is described with probe removal (Section 4.1), phase ramp removal (Section 4.2), and power correction (Section 4.3). Section 5 and Section 6 use simulation and experimental results to validate the method. Section 7 concludes.
2. Forward model
Figure 1 shows the schematic diagram of our 3D CMI method. In a single-shot CMI experiment, the sample is illuminated with a coherent light beam that is restricted by a pinhole. The sample exit wave then propagates to the modulator, interacts with the modulator, and propagates to the detector where the diffraction pattern is recorded. The whole forward model of a coherent modulation imaging measurement can be expressed as:
where $I_{\theta }$ denotes the measured diffraction pattern at projection angle $\theta$ , $\psi _{\theta }$ denotes the sample exit wave, and $M$ is the modulator function; $\mathcal {F}$ denotes the wave propagation operator, $z1$ is the sample-to-modulator distance, $z2$ is the modulator-to-detector distance, $\cdot$ denotes the element-wise multiplication. In CMI experiments, the sample is illuminated by a finite probe, therefore the support region constraints can be applied in reconstruction. Under thin sample conditions, the sample exit wave can be modeled as a multiplication of the probe and the sample:Here we use $O_{\theta }$ to denote the sample’s complex transmission function at projection angle $\theta$. The probe function $P$ is independent of $\theta$ as we assume the illumination is constant during data acquisition. The sample’s transmission function $O$ can be further defined as a function of the sample’s complex refractive index $n$:
where $i$ is the imaginary unit, $\lambda$ is the wavelength, $(x,y,z)$ is the spatial coordinate of the sample; $\delta$ and $\beta$ are real and imaginary parts of the complex refractive index, representing the phase-shift and absorption properties of the sample, respectively. By combining Eq. (1)–(4), the whole forward model of the recorded diffraction pattern $I_{\theta }$ can be summarized as:3. Reconstruction algorithm of CMI
Three-dimensional phase imaging usually contains two steps: phase retrieval and tomographic reconstruction. Before applying tomographic reconstruction, iterative phase retrieval algorithm is performed to obtain the sample exitwaves from measured diffraction patterns. Figure 2(a) illustrates the phase retrieval algorithm of CMI, where the running estimate of the wave at $j^{th}$ iteration for projection angle $\theta$ is denoted as $\psi _{\theta,j}$, $\psi _{\theta,j}^{m}$, $\psi _{\theta,j}^{M}$, and $\psi _{\theta,j}^{D}$ for exitwave at support plane, the incident and exit wave of the modulator, and the wave at the detector plane, respectively. The algorithm starts with an initial guess of the wavefield at the support plane, where support constraint is applied. Reconstruction is achieved by iteratively propagating the wavefield estimate between the support plane and the detector plane, via the modulator, and applying constraints.
Fig. 2. Flowchart of the reconstruction algorithm. (a) Phase retrieval algorithm of CMI; (b) projection extraction procedure for further tomographic reconstrucion.
In the support plane, the wavefield is known to have non-zero values within a finite region. For the $j^{th}$ iteration, the support region constraint is on the estimated wavefield according to:
In the detector plane, modulus constraint is applied to the estimated wavefield $\psi _{\theta,j}^{D}$. The modulus of the estimated detector wave is replaced by the square root of the measured pattern $I_{\theta }$, while the phase is kept unchanged, yielding the revised detector wave $\hat {\psi }_{\theta,j}^{D}$:
In the modulator plane, phase modulation is applied to the estimated incident wave $\psi _{\theta,j}^{m}$ during forward propagation and removed from the revised modulator exit wave $\hat {\psi }_{\theta,j}^{M}$ during backward propagation, according to the following equations:
It should be noted that the modulator function $M$ is prior knowledge, which is usually calibrated with ptychography [26]. The algorithm ends until either a predetermined number of iterations is completed or a preset value of some error metric has been achieved. The wave propagation operator $\mathcal {F}$ can be angular spectrum propagation or Fourier transform, depending on the experimental setup. For optical near-field configuration, angular spectrum propagation is used in the simulations and experiments below.
4. Projection extraction method
To get the sample’s 3D refractive indices, the illumination probe should be removed from the sample’s exitwaves. In addition, inherent ambiguities in phase reconstruction, including phase ramps and phase wrapping, constitute obstacles for accurate quantitative tomographic reconstruction. Another undesirable factor is intensity fluctuation, which may be caused by the instability of the light source and the acquisition setting of the detector. Therefore, we propose a projection extraction method to solve all these problems, detailed descriptions of our method are as follows.
4.1 Probe removal
The object’s transmission function $O_{\theta }$ is obtained by dividing the exitwave $\psi _{\theta }$ with the probe function $P$:
In most CDI experiments, the probe function is not available. But in CMI, the probe function $P$ can be recovered from a diffraction measurement with the sample removed. When recovering the probe function, the pinhole plane is selected as the support plane, therefore the complex wavefront at the pinhole plane is recovered. The recovered wavefront is then propagated to the sample plane, giving the illumination probe function $P$.
4.2 Phase ramp removal
Although the intensity and phase can be directly extracted from the reconstructed object function, there are phase ramps that must be removed before tomographic reconstruction to get accurate quantitative results [14,17]. The reconstructed object function $O_{\theta }$ can additionally have constant and linear phase terms:
Here $t_{\theta }(x,y)$ denotes the unambiguous sample transmissivity. The linear phase terms $a_{\theta }$ and $b_{\theta }$ are present in reconstructions if the diffraction patterns are not centered. The constant phase term $c_{\theta }$ arises from an inherent ambiguity in phase retrieval reconstructions and can take any value. The removal of these terms is detailed as follows.
4.2.1 Removal of linear phase terms
To avoid the linear phase terms $a_{\theta }$ and $b_{\theta }$, the measured diffraction patterns $\{I_{\theta }(x',y')\}$ should be aligned before phase retrieval reconstruction, here $(x',y')$ is the diffraction plane coordinate. During the experiment, a diffraction pattern without the sample should be recorded. Regions of free space adjacent to the sample should be measured, for this region the phase should be zero. Because the diffraction effect of the modulator is much stronger than that of the sample, and the light in free space is only diffracted by the modulator, the pixel values near the edge are nearly consistent in all diffraction patterns. To align the diffraction patterns, we apply a subpixel image registration algorithm [46] to the consistent region, a step-by-step description is as follows:
- 1) Find a windowed region $I_{\theta }(x_{c}',y_{c}')$ that is nearly consistent in all diffraction patterns $\{I_{\theta }(x',y')\}$, here $(x_{c}',y_{c}')$ denotes a subset of $(x',y')$.
- 2) For projection at angle $\theta$, extract the windowed region $I_{p}(x_{c}',y_{c}')$ and $I_{\theta }(x_{c}',y_{c}')$ from $I_{p}(x',y')$ and $I_{\theta }(x',y')$, respectively. Here $I_{p}(x',y')$ is the diffraction pattern without the sample, and the extracted regions will be used to derive global drifts.
- 3) Estimate the global drifts $d$ by the registration algorithm of the extracted regions with subpixel accuracy.
- 4) Shift the diffraction pattern $I_{\theta }(x',y')$ by the negative of the estimated value of $d$ from step 3).
- 5) Continue step 2)-4) until all projections have been processed.
4.2.2 Removal of constant phase term
Let us define $B_{\theta }(x,y)$ as the phase of $O_{\theta }(x,y)$, it is commonly computed as $B_{\theta }(x,y)=arg(O_{\theta }(x,y))$. $B_{\theta }$ then has a shift term:
Here $c_{\theta }$ is a real-value term that varies among all projections. We refer to $B_{\theta }(x,y)$ as non-normalized phase projection data. The influence of the term $c_{\theta }$ must be subtracted from $B_{\theta }(x,y)$ to obtain $B_{\theta }''(x,y)$. To do this, we also rely on having measured a region of free space $B_{\theta }(x_{c},y_{c})$, which is supposed to have two useful properties:
- 1) As a subset of the phase $B_{\theta }(x,y)$, the pixel values within $B_{\theta }(x_{c},y_{c})$ satisfies the previous formula, therefore: $B_{\theta }(x_{c},y_{c})=B_{\theta }'(x_{c},y_{c})+c_{\theta }$.
- 2) As a region of free space, the corrected pixel values $B_{\theta }'(x_{c},y_{c})$ is supposed to be equal to zeros for all phase projections: $B_{\theta }'(x_{c},y_{c})=0$.
Considering the first property, we select a projection $B_{0}(x,y)$ as the reference projection, and pixel values of other projections within $B_{\theta }(x_{c},y_{c})$ should be corrected to be equal to $B_{0}(x_{c},y_{c})$ according to:
Here $\gamma =\frac {\sum B_{0}(x_c,y_c)B_{\theta }(x_c,y_c)}{\sum |B_{\theta }(x_c,y_c)|^{2}}=\frac {\sum B_{0}(x_c,y_c)B_{0}(x_c,y_c)+c_{\theta }}{\sum |B_{0}(x_c,y_c)+c_{\theta }|^{2}}=\frac {B_{0}(x_c,y_c)}{B_{0}(x_c,y_c)+c_{\theta }}$. After this operation, $B_{\theta }(x_c,y_c)$ is corrected because $\gamma B_{\theta }(x_c,y_c)=\frac {B_{0}(x_c,y_c)}{B_{0}(x_c,y_c)+c_{\theta }}\cdot (B_{0}(x_c,y_c)+c_{\theta }) = B_{0}(x_c,y_c)$.
Considering the second property, the averaged pixel value of $B_0(x_c,y_c)$ should be subtracted from each projection, we can get the normalized phase projection according to:
Here $\overline {B}_{0c}$ is the averaged pixel value of $B_{0}(x_c,y_c)$, the formula is divided by $\gamma$ because the pixel values of sample have been magnified $\gamma$ times in Eq. (14).
It should be noted that the removal of constant phase is operated on the sample plane after phase retrieval, while the removal of linear phase terms is operated on the diffraction plane before phase retrieval. In addition to the phase ramps, phase wrap occurs when phase shift exceeds $2\pi$, which is usually the case when measuring thick samples. To solve this, we apply phase unwrapping method [47,48] to the normalized phase projection $B_{\theta }''(x,y)$ to obtain $\tilde {B}_{\theta }(x,y)$ that will be further used for tomographic reconstruction.
4.3 Power correction
During data acquisition, the instability of the light source can cause fluctuation in the power of diffraction patterns. Therefore, a power correction step is performed to the reconstructed amplitude $A_{\theta }(x,y)$. To do this, we select an intensity projection as the reference projection $A_0(x,y)$, and extract a region of free space $A_0(x_c,y_c)$. The power of other intensity projections $A_{\theta }(x,y)$ is corrected according to:
Here $\overline {A}_{0c}$, $\overline {A}_{\theta c}$ are the averaged pixel values of $A_0(x_c,y_c)$ and $A_{\theta }(x_c,y_c)$.
After the projection extraction process, we can obtain the corrected phase and intensity projections $\tilde {B}_{\theta }$ and $\tilde {A}_{\theta }$, and then the $\delta$ and $\beta$ projections are computed according to:
The extracted projections $\delta _{\theta }$ and $\beta _{\theta }$ are related to the 3D volume of either quantity by a Radon transform. Once projections at each angle have been reconstructed, standard tomographic reconstruction methods, such as the filtered back-projection algorithm (FBP), the algebraic reconstruction technique (ART) [5], and the simultaneous iterative reconstruction technique (SIRT) [49,50], can be used to obtain the 3D volumes that provide quantitative information about the phase-shift and absorption properties of the sample. In this work, the 3D reconstruction is solved with the SIRT algorithm.5. Simulation
In this section, simulations are carried out to verify and analyze the performance of the proposed 3D CMI method. The tomographic reconstruction was performed using the ASTRA toolbox [51,52], which is an open-source, GPU-accelerated software platform for tomographic algorithms.
A $200^3$-voxel discretized complex-valued 3D Shepp-Logan phantom was used to produce 180 transmission projections over 180 degrees. In the generated complex-valued phantom, different structures are used to simulate the real and imaginary parts, central slices of the simulated phantom are shown in Figs. 3(a) and (b). To evaluate the performance of our method, we simulated a phase-wrapped sample, where the $\beta$ values range from 0 to $10^{-5}$ and the $\delta$ values range from 0 to $10^{-4}$. The diffraction patterns are then generated according to the forward model of CMI, where the support at the projection plane was $240$ pixels in diameter, the modulator that randomly produces $0-\pi$ phase delay was placed at $20 mm$ from the support, and the detector with $512\times 512$ pixels of a size $5.5 \mu m$ was placed at $10 mm$ from the modulator. The simulated wavelength was $637 nm$.
Fig. 3. Complex-valued phantom used in simulations: (a) central slice for the simulated $\beta$ volume, (b) central slice for the simulated $\delta$ volume; simulation reconstruction results. Central slice of the reconstruced $\delta$ volume (c) without constant phase removal, (d) with constant phase removal.
5.1 Processing of projections
The generated diffraction patterns are first processed with 100 iterations of CMI algorithm, and the illumination probe is removed. For the recovered CMI results, the tomographic reconstruction is mainly affected by the constant phase term and phase wrapping, because intensity fluctuation and drifts in diffraction pattern are not simulated. To validate the removal of constant phase term (RCP), tomographic reconstructions with and without this step are compared. Figure 3(c) shows a central slice of the reconstructed $\delta$ volume without RCP, from which we can see obvious circular background. Voxel values within this region are supposed to be zero, while it is not the case in this reconstruction. After RCP, constant shifts in phase projections are corrected, and the unwanted background in tomographic reconstruction is successfully removed, as shown in Fig. 3(d).
5.2 Convergence analysis
We then evaluate the convergence of the 3D CMI method. It is obvious that the quality of tomographic reconstruction will be affected by the quality of CMI reconstruction. However, CMI reconstruction often contains noise, which not only originates from the Poisson noise in measured signal, but also comes from noise due to inconvergence. The convergence of CMI is affected by various experimental parameters (such as the size of support region, the distances from support to modulator and modulator to detector, and the phase depth of modulator) and reconstruction parameters (such as the iteration number), which has been analyzed in detail by He et al. [35]. In our work, we choose to control the significance of noise by running the CMI algorithm with different iterations, as the noise is supposed to be reduced with the increasing of iterations.
We run the CMI algorithm with iteration number N = 60, 70, 80, 90, and 100, to obtain the projection images. For the tomographic reconstruction, we run 100 iterations of the SIRT algorithm. Figure 4 compares the reconstructed phase and intensity tomograms when N = 60, 80, and 100, from which we can clearly see improvements in reconstruction accuracy with the increase of CMI iteration. The results when N = 60 is severely degraded by artifacts and blurring. The results after 100 CMI iterations are significantly better and there are no artifacts visible. For quantitative analysis, we calculate relative root-mean-square error (RRMSE) between the reconstructed tomograms and the original phantom according to:
Fig. 4. Simulation reconstruction results. Central slice of the reconstruced $\beta$ volume with 100 SIRT iterations and (a) 60 CMI iterations, (b) 80 CMI iterations, and (c) 100 CMI iterations. Central slice of the reconstruced $\delta$ volume with 100 SIRT iterations and (d) 60 CMI iterations, (e) 80 CMI iterations, and (f) 100 CMI iterations. (g) RRMSE of the tomographic reconstructions with the iteration number of CMI. (h) PSNR and NRMSE of the CMI reconstruction results.
Figure 4(h) plots the evolution of peak signal to noise ratio (PSNR) and normalized root mean square error (NRMSE) over 100 CMI iteration. The PSNR is calculated between the reconstructed intensity and the original one, it is usually used to measure image quality [40]. The NRMSE is used to monitor the convergence of phase retrieval, which is calculated according to:
5.3 Reconstruction with limited angles
In this section, we compare the tomographic reconstruction performance with different number of angular projections $N_{\theta }$. The tomographic dataset was collected with 180 projection angles over 180 degrees with a step size of 1 degree. To further analyze 3D CMI, we tried less angular projections: 90, 60, 30, and 15 (also over 180 degrees). We run 100 iterations of CMI algorithm for each projection and 100 iterations of SIRT algorithm for tomographic reconstruction. Figure 5 shows the reconstructed phase slice with fewer angular projections. For quantitative analysis, we calculated the PSNR and mean square error (MSE) between the reconstructed central phase slice and the original one. The PSNR of the reconstructed central phase slice with 90, 60, 30, and 15 projections are 23.97, 23.92, 22.97, and 20.78 (dB), respectively; and the MSEs are 0.0040, 0.0041, 0.0051, and 0.0084, respectively. Clearly, they are not as good as the reconstructions with complete dataset (see Fig. 4(f)). However, the results show that acceptable tomographic reconstruction can be obtained using only 30 angular projections. The errors become visually obvious with only 15 angular projections.
Fig. 5. Simulation reconstruction results. Central slice of the reconstruced $\delta$ volume with (a) 90 projections, (b) 60 projections, (c) 30 projections, and (d) 15 projections.
6. Experiment
As a proof-of-principle experiment, we present optical results of an insect wing sample. The experimental setup is shown in Fig. 6(a). The collimated beam of 637nm (Coherent, OBIS 637LX) is restricted by a 1.5mm diameter pinhole near the sample. The sample is adhered to the top of a glass capillary and mounted on a rotation stage (DAHENG, GCD-012100M). The sample exit wave passes through the beam splitter, the modulator, and reaches the detector (IMPERX 6620B, $5.5\mu m\times 5.5\mu m$ pixel size) where the diffraction patterns are recorded. The sample-to-modulator distance z1 is $100mm$ and the modulator-to-detector distance z2 is $28.6mm$. The sample is rotated a total of 180 degrees and one diffraction pattern of $512\times 512$ pixels is recorded every one degree. A diffraction pattern with the sample removed is also recorded.
Fig. 6. (a) Optical experimental setup for 3D CMI: a collimated beam passes through a pinhole and illuminated on the sample which is mounted on a rotation stage, the modulator is placed just before the camera where diffraction patterns are recorded. The (b) amplitude and (c) phase of the modulator’s transmission function calibrated using ptychography.
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$ $\mu m$ $\times$ $16$ $\mu m$. It is worth mentioning that CMI has no strict requirement on the distribution of the modulator, as continuous phase modulation [33] and amplitude modulation [34] are also reported. The modulator’s transmission function is calibrated using ptychography [53,54] with the ePIE algorithm [25], and the reconstructed amplitude and phase are shown in Figs. 6(b) and (c), respectively.
6.1 Processing of projections
Before phase retrieval reconstruction, the diffraction patterns need to be aligned to avoid the linear phase terms according to Section 4.2.1. Figures 7(a) and (b) show the diffraction patterns with sample removed and with the sample, respectively. Removal of linear phase terms (RLP) involves selecting a windowed region to derive the global drifts, as shown in Figs. 7(c) and (d).
Fig. 7. Recorded diffraction patterns with (a) the sample removed and (b) with the sample; (c) and (d) are the selected windowed region from (a) and (b), respectively.
After RLP, the aligned diffraction patterns are processed with 100 iterations of the CMI algorithm. One of the reconstructed results are shown in Figs. 8(a) and (b), from which we can clearly see the phase wraps in the recovered phase image (Fig. 8(b)) and the illumination probe in both. The reconstructed amplitude and phase are then processed according to the projection extraction procedure. The final extracted projections are shown in Figs. 8(c) and (d).
Fig. 8. Central $400\times 400$ pixels of the reconstructed (a) amplitude and (b) phase from a diffraction pattern after 100 CMI iterations. Projection extraction results for the (c) $\beta$ projection and (d) $\delta$ projection. Tomographic slice of the reconstructed $\beta$ volume (e) with linear phase removal, (f) without linear phase removal; tomographic slice of the reconstructed $\delta$ volume (g) with linear phase removal, (h) without linear phase removal.
Because the sample is not located perfectly at the rotation center, there are displacements in the horizontal dimension that must be corrected. For the $\beta$ projections, we apply the center of mass (COM) method [55] in a region of interest (ROI) to align the projections in horizontal direction. The COM method returns the displacement in pixels for each $\beta$ projection. These displacements are then used to align the corresponding $\delta$ projections.
6.2 Tomographic reconstruction
For the tomographic reconstruction, we run 100 iterations of the SIRT algorithm with GPU acceleration. Minimum constraints (the voxel values must be larger than zero) are enforced to help the convergence. The projections are cropped to their central $200\times 200$ pixels. Figures 8(e-h) compare the tomographic reconstruction slice with and without RLP, from which we can see the artifacts are reduced obviously after RLP. Figure 8(f) indicates that the intensity reconstruction is also affected, this is because RLP aligns the diffraction patterns. Without RLP, the probe may not be removed correctly from the sample’s exit wave, resulting in ring artifacts in the reconstructed tomogram.
The final tomographic results of $200^{3}$ voxels are shown in Fig. 9, including the 3D volume rendering of the reconstructed (a) $\beta$ and (b) $\delta$ volumes. Figures 9(c), (d), and (e) are the y-z, x-y, and x-z projections of the 3D $\beta$ reconstruction, while Figs. 9(f), (g), and (h) are the y-z, x-y, and x-z projections of the 3D $\delta$ reconstruction, respectively. Compared with the 3D $\beta$ reconstruction, the 3D $\delta$ reconstruction shows higher contrast because of the higher phase sensitivity of the transparent tissue. The real and imaginary parts of the reconstructed refractive index at different axial depths are shown in Fig. 10. Figures 10(a-c) show axial slices from the reconstructed $\beta$ volume, and Figs. 10(d-f) the corresponding slices from the reconstructed $\delta$ volume. The separation between the slices is 110um. An animation of the full stack of 200 axial slices can be found in Visualization 1, which demonstrates the tomographic ability of 3D CMI. Here the 200 slices are plotted sequentially and are labelled by their slice number.
Fig. 9. Experimental reconstruction results. (a) 3D rendering of the reconstructed $\beta$ volume, (b) 3D rendering of the reconstructed $\delta$ volume; (c), (d), and (e) show the projections of the reconstructed $\beta$ volume; (f), (g), and (h) show the projections of the reconstructed $\delta$ volume.
Fig. 10. Selected axial slice from the reconstructed $\beta$ volume, (a)-(c), and the $\delta$ volume, (d)-(f). The slices correspond to planes 80, 100, and 120 out of a total of 200. The separation of the selected slices is 110$\mu m$ (see Visualization 1).
It should be noted that the voxel size of our reconstruction is 5.5$\mu m$, which is the same as the pixel size of the employed detector because we adopted the near-field configuration and plane-wave illumination. The spatial resolution of 3D reconstruction is decided by the spatial resolution of 2D projections and tomographic reconstruction accuracy. For optical near-field CMI configuration, the spatial resolution of 2D reconstruction is mainly limited by the detector pixel size [27]. Modifications to the experimental setup, such as physically reducing the pixel size of a sensor, and using a divergent illumination [33], can improve the spatial resolution. For far-field geometry, the spatial resolution can be much smaller than the detector pixel itself, and high-resolution images can be obtained using radiations with shorter wavelengths such as x-rays. In the existing literatures, Zhang et al. [31] have achieved a 2D spatial resolution of 37nm with CMI using x-rays of energy 6.2 keV.
In addition, we should note that there are still noticeable artifacts in the tomographic reconstruction, which also affects the reconstruction resolution. These artifacts may arise from linear and angular motion errors of the experimental tomography setup. We will investigate on optimizing these errors in our future work. Nevertheless, the 3D CMI method is successfully experimentally demonstrated.
7. Conclusion
In conclusion, we present a 3D reconstruction procedure for sample’s complex refractive indices from 2D coherent modulation measurements. We propose a projection extraction method and validate it with numerical simulations and optical experiments. In the simulations, tomographic reconstruction results with and without constant phase removal are compared; the convergence of 3D CMI is analyzed and the results are quantified with RRMSE and PSNR; the 3D reconstruction performance with limited projections is also explored. In the experiments, we present 3D reconstruction results of an insect wing sample, in which the removal of linear phase term is validated; 3D volume rendering and axial slices of the reconstructed complex refractive indices are given to demonstrate the tomographic ability of 3D CMI.
The simulation and experimental results presented suggest that the proposed projection extraction method is able to obtain the sample transmissivity, remove the constant and linear phase terms, and calibrate the light source power, thus paving the way from 2D CMI reconstruction to 3D tomographic reconstruction. Although these methods are demonstrated with optical near-field configuration and plane-wave illumination, they are in principle applicable to other tomographic configurations, including for example far-field x-ray CMI [31], and divergent illumination CMI [33].
Further improvements to our work will focus on the tomographic alignment and reconstruction, including the optimization of tomography setup to remove possible linear and angular motion errors, and the use of compressive sensing [56] to reduce the number of projections required for reconstruction.
Funding
National Natural Science Foundation of China (61975205, 6207522, 62131011); Youth Innovation Promotion Association of the Chinese Academy of Sciences (2017489); Natural Science Foundation of Hebei Province (F2018402285); Funded Project of Hebei Province Innovation Capability Improvement Plan (No. 20540302D); Fundamental Research Funds for the Central Universities; Fusion Foundation of Research and Education of CAS.
Disclosures
The authors declare no conflicts of interest.
Data availability
Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.
References
1. Z. Yousaf, P. Withers, and P. Potluri, “Compaction, nesting and image based permeability analysis of multi-layer dry preforms by computed tomography (CT),” Compos. Struct. 263, 113676 (2021). [CrossRef]
2. S. W. Paddock and K. W. Eliceiri, “Laser Scanning Confocal Microscopy: History, Applications, and Related Optical Sectioning Techniques,” in Confocal Microscopy: Methods in Molecular Biology, vol. 1075, S. W. Paddock, ed. (Springer, 2014), pp. 9–47.
3. J. Kweon, S.-J. Kang, Y.-H. Kim, J.-G. Lee, S. Han, H. Ha, D. H. Yang, J.-W. Kang, T.-H. Lim, O. Kwon, J.-M. Ahn, P. H. Lee, D.-W. Park, S.-W. Lee, C. W. Lee, S.-W. Park, and S.-J. Park, “Impact of coronary lumen reconstruction on the estimation of endothelial shear stress: in vivo comparison of three-dimensional quantitative coronary angiography and three-dimensional fusion combining optical coherent tomography,” Eur. Hear. J. - Cardiovasc. Imaging 19(10), 1134–1141 (2018). [CrossRef]
4. M. Engel, L. Kasper, B. Wilm, B. Dietrich, L. Vionnet, F. Hennel, J. Reber, and K. P. Pruessmann, “T-Hex: Tilted hexagonal grids for rapid 3D imaging,” Magn. Reson. Med. 85(5), 2507–2523 (2021). [CrossRef]
5. A. C. Kak and M. Slaney, Principles of Computerized Tomographic Imaging (Society for Industrial and Applied Mathematics, 2001).
6. F. Charrière, A. Marian, F. Montfort, J. Kuehn, T. Colomb, E. Cuche, P. Marquet, and C. Depeursinge, “Cell refractive index tomography by digital holographic microscopy,” Opt. Lett. 31(2), 178 (2006). [CrossRef]
7. V. Balasubramani, A. Kus, H.-Y. Tu, C.-J. Cheng, M. Baczewska, W. Krauze, and M. Kujawinska, “Holographic tomography: techniques and biomedical applications [Invited],” Appl. Opt. 60(10), B65 (2021). [CrossRef]
8. D. Kim, S. Lee, M. Lee, J. Oh, S.-A. Yang, and Y. Park, “Holotomography: Refractive Index as an Intrinsic Imaging Contrast for 3-D Label-Free Live Cell Imaging,” in Advanced Imaging and Bio Techniques for Convergence Science: Advances in Experimental Medicine and Biology, vol. 1310, J. K. Kim, J. K. Kim, and C.-G. Pack, eds. (Springer, 2021), pp. 211–238.
9. L. Tian, J. C. Petruccelli, Q. Miao, H. Kudrolli, V. Nagarkar, and G. Barbastathis, “Compressive x-ray phase tomography based on the transport of intensity equation,” Opt. Lett. 38(17), 3418 (2013). [CrossRef]
10. Y. Wen and A. Asundi, “3D profile measurement for stepped microstructures using region-based transport of intensity equation,” Meas. Sci. Technol. 30(2), 025202 (2019). [CrossRef]
11. 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]
12. Y. Nishino, Y. Takahashi, N. Imamoto, T. Ishikawa, and K. Maeshima, “Three-Dimensional Visualization of a Human Chromosome Using Coherent X-Ray Diffraction,” Phys. Rev. Lett. 102(1), 018101 (2009). [CrossRef]
13. C. T. Putkunz, M. A. Pfeifer, A. G. Peele, G. J. Williams, H. M. Quiney, B. Abbey, K. A. Nugent, and I. McNulty, “Fresnel coherent diffraction tomography,” Opt. Express 18(11), 11746 (2010). [CrossRef]
14. T. Ramos, J. S. Jørgensen, and J. W. Andreasen, “Automated angular and translational tomographic alignment and application to phase-contrast imaging,” J. Opt. Soc. Am. A 34(10), 1830 (2017). [CrossRef]
15. T. Ramos, B. E. Grønager, M. S. Andersen, and J. W. Andreasen, “Direct three-dimensional tomographic reconstruction and phase retrieval of far-field coherent diffraction patterns,” Phys. Rev. A 99(2), 023801 (2019). [CrossRef]
16. M. Dierolf, A. Menzel, P. Thibault, P. Schneider, C. M. Kewish, R. Wepf, O. Bunk, and F. Pfeiffer, “Ptychographic X-ray computed tomography at the nanoscale,” Nature 467(7314), 436–439 (2010). [CrossRef]
17. 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 (2011). [CrossRef]
18. M. Holler, A. Diaz, M. Guizar-Sicairos, P. Karvinen, E. Färm, E. Härkönen, M. Ritala, A. Menzel, J. Raabe, and O. Bunk, “X-ray ptychographic computed tomography at 16 nm isotropic 3D resolution,” Sci. Rep. 4(1), 3857 (2015). [CrossRef]
19. P. Li and A. Maiden, “Multi-slice ptychographic tomography,” Sci. Rep. 8(1), 2049 (2018). [CrossRef]
20. J. R. Fienup, “Phase retrieval algorithms: a comparison,” Appl. Opt. 21(15), 2758 (1982). [CrossRef]
21. J. R. Fienup and C. C. Wackerman, “Phase-retrieval stagnation problems and solutions,” J. Opt. Soc. Am. A 3(11), 1897 (1986). [CrossRef]
22. M. Guizar-Sicairos and J. R. Fienup, “Understanding the twin-image problem in phase retrieval,” J. Opt. Soc. Am. A 29(11), 2367 (2012). [CrossRef]
23. D. Yang, J. Zhang, Y. Tao, W. Lv, S. Lu, H. Chen, W. Xu, and Y. Shi, “Dynamic coherent diffractive imaging with a physics-driven untrained learning method,” Opt. Express 29(20), 31426 (2021). [CrossRef]
24. J. M. Rodenburg and H. M. L. Faulkner, “A phase retrieval algorithm for shifting illumination,” Appl. Phys. Lett. 85(20), 4795–4797 (2004). [CrossRef]
25. A. M. Maiden and J. M. Rodenburg, “An improved ptychographical phase retrieval algorithm for diffractive imaging,” Ultramicroscopy 109(10), 1256–1262 (2009). [CrossRef]
26. 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 (2013). [CrossRef]
27. W. Xu, H. Lin, H. Wang, and F. Zhang, “Super-resolution near-field ptychography,” Opt. Express 28(4), 5164 (2020). [CrossRef]
28. W. Lv, J. Zhang, H. Chen, D. Yang, T. Ruan, Y. Zhu, Y. Tao, and Y. Shi, “Resolution-enhanced ptychography framework with equivalent upsampling and precise position,” Appl. Opt. 61(10), 2903–2909 (2022). [CrossRef]
29. 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]
30. F. Zhang and J. M. Rodenburg, “Phase retrieval based on wave-front relay and modulation,” Phys. Rev. B 82(12), 121104 (2010). [CrossRef]
31. 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]
32. 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 (2018). [CrossRef]
33. 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]
34. X. Pan, C. Liu, and J. Zhu, “Coherent amplitude modulation imaging based on partially saturated diffraction pattern,” Opt. Express 26(17), 21929 (2018). [CrossRef]
35. 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]
36. B. Wang, Q. Wang, W. Lyu, and F. Zhang, “Modulator refinement algorithm for coherent modulation imaging,” Ultramicroscopy 216, 113034 (2020). [CrossRef]
37. 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 (2020). [CrossRef]
38. B. Wang, Z. He, and F. zhang, “Coherent modulation imaging using unknown modulators,” Opt. Express 29(19), 30035 (2021). [CrossRef]
39. Y. Geng, X. Wen, X. Zhou, Y. Li, J. Tan, W. Ding, S. Liu, and Z. Liu, “Multi-rotation coherent imaging by a phase mask,” Opt. Lasers Eng. 139, 106511 (2021). [CrossRef]
40. J. Zhang, D. Yang, Y. Tao, Y. Zhu, W. Lv, D. Miao, C. Ke, H. Wang, and Y. Shi, “Spatiotemporal coherent modulation imaging for dynamic quantitative phase and amplitude microscopy,” Opt. Express 29(23), 38451 (2021). [CrossRef]
41. C. Xu, H. Pang, A. Cao, and Q. Deng, “Coherent modulation imaging realized by coaxial shift of an amplitude mask,” Opt. Commun. 513, 128118 (2022). [CrossRef]
42. Y. Xu, X. Pan, C. Liu, and J. Zhu, “Single-shot phase reconstruction based on beam splitting encoding and averaging,” Opt. Express 29(26), 43985 (2021). [CrossRef]
43. 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]
44. 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]
45. 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 (2015). [CrossRef]
46. M. Guizar-Sicairos, S. T. Thurman, and J. R. Fienup, “Efficient subpixel image registration algorithms,” Opt. Lett. 33(2), 156 (2008). [CrossRef]
47. Z. Zhao, H. Zhang, Z. Xiao, H. Du, Y. Zhuang, C. Fan, and H. Zhao, “Robust 2D phase unwrapping algorithm based on the transport of intensity equation,” Meas. Sci. Technol. 30(1), 015201 (2019). [CrossRef]
48. Z. Zhao, H. Zhang, C. Ma, C. Fan, and H. Zhao, “Comparative study of phase unwrapping algorithms based on solving the Poisson equation,” Meas. Sci. Technol. 31(6), 065004 (2020). [CrossRef]
49. A. van der Sluis and H. van der Vorst, “SIRT- and CG-type methods for the iterative solution of sparse linear least-squares problems,” Linear Algebr. its Appl. 130, 257–303 (1990). [CrossRef]
50. J. Gregor and T. Benson, “Computational Analysis and Improvement of SIRT,” IEEE Trans. Med. Imaging 27(7), 918–924 (2008). [CrossRef]
51. W. van Aarle, W. J. Palenstijn, J. De Beenhouwer, T. Altantzis, S. Bals, K. J. Batenburg, and J. Sijbers, “The ASTRA Toolbox: A platform for advanced algorithm development in electron tomography,” Ultramicroscopy 157, 35–47 (2015). [CrossRef]
52. W. van Aarle, W. J. Palenstijn, J. Cant, E. Janssens, F. Bleichrodt, A. Dabravolski, J. De Beenhouwer, K. Joost Batenburg, and J. Sijbers, “Fast and flexible X-ray tomography using the ASTRA toolbox,” Opt. Express 24(22), 25129 (2016). [CrossRef]
53. 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]
54. R. Ma, S.-Y. Zhang, T.-H. Ruan, Y. Tao, H.-Y. Wang, and Y.-S. Shi, “Scanning-Position Error-Correction algorithm in Dual-Wavelength Ptychographic Microscopy,” Chin. Phys. Lett. 37(4), 044201 (2020). [CrossRef]
55. S. Azevedo, D. Schneberk, J. Fitch, and H. Martz, “Calculation of the rotational centers in computed tomography sinograms,” IEEE Trans. Nucl. Sci. 37(4), 1525–1540 (1990). [CrossRef]
56. E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Phys. Med. Biol. 53(17), 4777–4807 (2008). [CrossRef]
