Electron tomography can reveal three-dimensional (3D) structure of biological specimens David and Klug (1968); Fernandez et al. (2010) and inorganic materials Yu et al. (2011); Sai et al. (2013); Ercius et al. (2009) across the nano- to meso- scale. However, 3D reconstruction commonly suffers from an insufficient number of projections and a limited angular range inside the microscope, commonly referred to as the missing wedge Jiang et al. (2017). To overcome these experimental limitations, advancements in tomographic reconstruction algorithms now utilize compressed sensing (CS) methods that leverage the notion of sparsity to provide superior 3D resolution with limited angles and at lower doses Donoho (2006); Candés et al. (2006); Saghi et al. (2011). While compressed sensing algorithms provide higher quality reconstructions, they require substantially more computation time to complete and bottleneck electron tomography. Even worse, incomplete convergence and poor selection of tunable parameters are often discovered at the end of an arduous reconstruction—multiplying the total time required to produce a full 3D tomogram. Scientists would benefit from immediate feedback to examine intermediate results and optimize experimental parameters that accelerate the end-to-end scientific process Bicer et al. (2015).
Here, we have developed a dynamic CS framework that offers a 3D specimen reconstruction in real-time as projection data is collected. It enables direct feedback and on the fly optimization of experimental parameters. The reconstruction algorithm begins immediately upon acquiring the first projection and dynamically updates the 3D structure as new projections arrive—unlike traditional schemes which start after the experiment is complete. This means researchers can start analysis and characterization with high-fidelity tomograms before an experiment is complete. Using scanning transmission electron microscope (STEM) tomography Midgley and Weyland (2003), we demonstrate our method accelerates the final convergence by a factor of 2-3 over conventional CS and provides insight into 3D nanostructure within 62% of the total experimental acquisition time. Dynamic reconstruction reduced the reconstruction error for Au/SrTiO nanoparticles by 27% and converged 100% faster than a traditional approach. Moreover, researchers can use dynamic CS to manipulate the data-tolerance throughout the reconstruction and efficiently explore tunable parameters without having completely reset the algorithm. Implementing dynamic CS requires complete parallelization that includes the 3D total variation regularization for the isotropic norm. Tomograms () can reconstructed dynamically on modest multi-core laptops during an electron tomography experiment and larger reconstructions () are achievable with high performance computing.
The notion of sparsity has become widely used in signal processing and image reconstruction as a prior knowledge to regularize solutions in underdetermined problems. It was greatly popularized by the theory of compressed sensing Sidky et al. (2011) that demonstrates the possibility to accurately recover the 3D structure of specimens () from an insufficient number of projections () with -norm optimizations. One of the most representative sparsity-exploiting algorithms is the total variation minimization (TV-min), which was originally proposed for image denoising Rudin et al. (1992) and widely used to reduce tomographic artifacts for reconstructions from a limited number of projections Sidky et al. (2006); Leary et al. (2013); Guay et al. (2016). The technique can effectively remove noisy features while preserving the edges of the object by minimizing its gradient magnitude. In this work, we consider a constrained optimization problem defined as:
where and denote the TV and norms, and is a data-tolerance parameter that controls the trade-off between regularization (smoothness) and data fidelity. Here, a tomography experiment is formulated as an inverse problem, , where is the projection matrix that models the measurement physics. The constrained optimization problem can be solved with a combination of adaptive steepest-descent (ASD) to minimize TV and projection onto convex sets (POCS) to enforce data constraints—commonly referred to as (ASD-POCS) Sidky and Pan (2008). This constrained optimization provides physical meaning to the tunable regularization parameter,
, which can be initially estimated from the data qualityLiu et al. (2017). Across all ‘flavors’ of compressed sensing tomography, the optimization process begins after all data has been collected. The iterative process can take thousands of iterations and runs from hours to a full day before converging to the designed solution. Moreover, the regularization parameter (here, ) is often task-dependent and needs be to adjusted to produce the best image quality, further increasing the computation time for the reconstruction process.
Iii Dynamic Compressed Sensing
Dynamic reconstruction during data collection allows researchers early insight into 3D structure throughout a tomography experiment. Figure 1 highlights the overall framework for the dynamic CS algorithm. Instead of starting at the end of an experiment, the reconstruction task begins immediately when the first projection is available. As more projections are experimentally acquired, the new information is accommodated as additional constraints in the optimization process (Fig. 1a) and improves the reconstructed tomogram quality (Fig. 1b). Because the reconstruction process is continuously running throughout the entire data acquisition, which typically take several hours, dynamic compressed sensing is able to produce a high-quality reconstruction before or upon arrival of the final projection (Fig. 1c).
For quantitative assessment, a 3D synthetic tomogram ( voxels) derived from an experimental gold decorated strontium titanate (Au/SrTiO) Padgett et al. (2017) nanoparticle was used to simulate projections along a single axis of rotation. A tilt range of with a two-degree tilt increment () approximates the angular range commonly collected in an electron microscope. Samples cannot normally tilt to high angles due to beam shadowing caused by the specimen or holder geometry Banjak et al. (2018)
. Poisson noise was added to the projections to give a signal to noise ratio (SNR) of 100. New projections were added every 3 minutes during the reconstruction to mimic typical experimental acquisition rates.
Dynamic compressed sensing produces a 3D tomogram with increasing quality throughout the experiment and by the end is nearly identical to the true SrTiO specimen. We can verify the performance by comparing the reconstruction to the ground truth (Fig. 1a vs c) and measuring the tomogram’s root mean square error (RMSE) (Fig. 1
b). The reconstruction error reduces with increasing number of iterations and projections acquired. The RMSE curve has a non-smooth, staircase structure when new projections are incorporated into the data vector. Early in the data acquisition process (hr, 14 projections), the limited tilt range results in expected streaking and blurring of the tomogram. As more projections are gathered, additional features are imposed on the 3D volume. Internal voids within cubic SrTiO nanoparticles become visible and small Au nanoparticles on the surface become less elongated. Towards the end of an experiment ( hrs, 45 projections), the SrTiO reconstruction qualitatively matches the true object.
Iv Evaluating Convergence
As dynamic CS progresses, in both iterations and the number of projections, it reduces RMSE but eventually diverges from the optimal solution with minimal error and approaches a solution defined by the optimization problem (Eq. 1). Iterative algorithms (e.g. Kaczmarz, Landweber, or Cimmino Method) typically deviate from solutions with minimum RMSE when applied to noisy data Elfving et al. (2014). While RMSE is a useful quantitative measure to assess reconstruction performance, it does not match the visually desired solution Jiang et al. (2017). Moreover, computing RMSE requires knowledge of the true 3D specimen structure.
In real experiments the true solution is unknown and a tomogram’s RMSE cannot be measured. Instead, the progression of a reconstruction’s TV and data distance (DD) can be utilized to assess the convergence towards an optimal solution. Figure 2 plots RMSE, DD and TV vs. time for the Au-SrTiO phantom nanoparticle during a dynamic compressed sensing reconstruction. Throughout an experiment (shaded green), the data distance trends downward to the specified data tolerance, (red line) indicating stable convergence. The incorporation of new data creates sharp discontinuities in DD and TV. Unlike RMSE which drops with the addition of new projections, DD and TV momentarily rise sharply because ASD-POCS is attempting to minimize the distance between DD and by iteratively adjusting the weights between data fidelity and regularization. After the arrival of new data, the algorithm will sufficiently converge to a solution within 125 iterations. Dynamic CS performs best when there are enough iterations to satisfy its data tolerance constraint (DD ) before new projections are introduced. If additional projections are added too quickly, the overall convergence may drift (See Supplemental Fig. 8) and the algorithm will be unable to reach its optimal solution by experimental completion.
Dynamic CS reliably produces high reconstruction quality with stable convergence—even before the experiment is complete (e.g. 80%, 3.5 hrs). After the last projection is collected (shaded red) the RMSE, DD, and TV for dynamic compressed sensing converges 50% faster than traditional compressed sensing (blue curve) that starts reconstruction after all data is collected. Dynamic CS benefits from a significantly closer optimal solution that allows for faster convergence when new data arrives. The final reconstruction produced by dynamic CS is indistinguishable from traditional CS and the solution converges to the true object’s total variation (TV) and .
V Dynamic CS for Experimental Electron Tomography
Dynamic CS can accurately reconstruct experimental electron tomography data. Figure 3 shows dynamic CS reconstruction of platinum (Pt) nanoparticles decorated on a carbon nanofiber. The experimental data had a limited a tilt range of and a two-degree increment (). The data, first acquired by Padgett et al. Padgett et al. (2017); Levin et al. (2016), was introduced into the reconstruction every 3 minutes to simulate experimental acquisition rates. The true object is unknown for experimental data. However, the convergence of DD and TV exhibit similar behaviors to simulation studies (Fig 2b). That is, DD is progressively driven down to and TV is iteratively converging to a minima.
Internal cavities in the carbon nanofiber and the platinum nanoparticles become visible within half the time to complete the experiment (e.g. 2 hrs). Details of the full physical structure becomes clearly visible within 62% of the complete acquisition (47 projections). The small Pt-nanoparticles are well resolved and show minimum elongation due to the missing wedge. The full structure is clearly visible about an hour prior to completion. 3D visualizations of the final reconstructions are highlighted in Fig. 3a.
Vi Dynamic Parameter Tuning
Selecting often requires computing several reconstructions and ultimately relies on the scientist’s judgement. Here we show dynamic compressed sensing allows to be tightened (decreased) or loosened (increased) mid-reconstruction. Furthermore, the data-constraint can be reversibly adjusted—reflecting stable and convex convergence. Generally speaking, selection of the data consistency constraint, , depends on the SNR as it accommodates all sources of data inconsistency (e.g. noise) and ensures re-projections are within a given distance from the actual (experimental) data Zhang et al. (2011). Dynamic parameter tuning allows researchers to more efficiently dial in the optimal parameter value. If is too low the reconstruction appears noisy; if is too high, detail and resolution is degraded.
Figure 4 demonstrates that manipulating throughout the reconstruction consistently converges in stable results. Here, we reconstructed a fully sampled Au/SrTiO tilt series and waited until all of the optimization parameters (RMSE, DD, and TV) were fully converged prior to perturbing the data constraint. During iterations: 0–4,500 we reduced eight times by before increasing by +0.025 eight times during iterations 4,500–8,500. For this dataset, we found small perturbations () guarantee convergence within 100 iterations. Dynamic CS reliably converges to solutions defined by the final chosen. Whether incrementally increased or decreased, the final value of determines a nearly unique solution, which can be seen both visually (Fig 4a,b) and quantitatively in the RMSE 4c. Note that the minimal RMSE () retains grainy artifacts and does not produce a desirable reconstruction due to Au particle’s high intensities. As discussed by Jiang, et. al., the visually appealing result is generally obtained from a slightly larger Jiang et al. (2017). We observed similar data-tolerance properties, the visually desirable reconstruction (highlighted in green Fig. 4) occurs at .
Vii Parallelization & Performance
Computational efficiency is key to the success of dynamic compressed sensing. As discussed previously (Section IV), the reconstruction process should reach stable solutions before incorporating new projections. We emphasize that the overall computation should always be faster than the data acquisition, so that by performing the image reconstruction on the fly, we can obtain the reconstructed image almost right after the experiment finishes. This is seen in all the previous cases we have presented, in which RMSE converges to a plateau value before the arrival of a new projection. However, as the size of the object increases, the computational complexity grows as , where is the size of the object in each dimension; the experimental time, however, grows only as . Therefore, for large physical systems (), single laptop/desktop or workstation is not powerful enough for dynamic compressed sensing. To overcome this problem, we deployed high performance computing (HPC) resources at Theta, a Cray XC40 11.69 petaflops supercomputer at Argonne Leadership Computing Facility.
We use Massage Passing Interface (MPI) Clarke et al. (1994) to parallelize our code across different nodes. In particular, we distribute the 2D slices to different MPI processes. There are two dominant computations involved in our algorithm: (1) ART – for minimizing ; (2) TV – for calculating the TV gradient of the object, . For ART, the computation is independent for different slices and no data exchange is needed among the processes. For TV, each process only needs to send the first and last slices owned by that process to the two nearby processes respectively. The communication overhead thus is minimal. We expect our algorithm to scale efficiently in HPC supercomputers. Besides using MPI for the inter-node parallelism, for the intranode parallelism, we use OpenMP Dagum and Menon (1998) to parallelize the computation intensive loops.
Viii Discussion & Conclusion
Across all varieties of compressed sensing tomography, the optimization process begins after all data has been collected. In this study we demonstrated these methods are compatible with a dynamic, continuously running reconstruction framework which automatically accommodates newly measured data. The algorithm allows for real-time analysis of 3D specimens as an experiment progresses and produces a high-quality reconstruction upon completion of data collection. Although conventional reconstruction algorithms, such as weighted back projection already provide real-time reconstructions due to their simplicity, the quality is significantly inferior to compressed sensing for electron tomography—especially whenever there is limited to moderate number of projections.
Dynamic CS leverages the time it takes to experimentally acquire projection data in electron tomography to compute high-quality 3D specimen structure. Evaluating dynamic CS against a fully-sampled conventional CS reconstruction, tomogram error is reduced by 27% in a Au/STiO dataset. The successive intermediate updates promote early convergence and it was shown to reconstruct fine features within 62% of the entire experimental acquisition time. The data tolerance parameter () can be interactively tightened or loosened mid-process to mitigate the complexities associated with selecting optimal regularization. Our numerical results illustrate that the proposed dynamic approach practically accelerates the convergence rate of conventional CS algorithms up to a factor of 3. Large-scale reconstructions are achieved with a massively parallelized implementation that benefits from the use of a high-performance computing cluster. The proposed framework builds a foundation for developing more sophisticated algorithms that incorporate projection alignment Odstrcil et al. (2019)
and machine learning techniquesDing et al. (2019). This provides new opportunities for other applications such as hyperspectral or ptychographic tomography—currently active areas of research. Combined with real-time 3D visualization tools, dynamic CS can enable real-time electron tomography. Scientists can visualize intermediate results with the high-fidelity of compressed sensing as projections are gathered to directly assess detailed specimen structure and optimize experimental parameters.
This work was supported by the Department of Energy (DOE) Office of Science (DE-SC0011385) and Argonne Data Science Program (ADSP). Simulations made use of the Advanced Research Computing Technology Services’ shared high-performance computing at the University of Michigan and Argonne Leadership Computing Facility (ALCF) at the Argonne National Laboratory. ALCF is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. Hardware provided by NVidia GPU Science Grant.
Author contributions statement
Computational experiments and tomographic reconstruction was conducted by J.S., Y.J., and R.H. Code was optimized and parallelized by H.Z. Data analysis and interpretation was carried out by J.S., H.Z., M.H., Y.J, and R.H. All authors reviewed the manuscript.
- David and Klug (1968) R. David and A. Klug, Nature 217, 130 (1968).
- Fernandez et al. (2010) R. Fernandez, B. Zuber, U. Elisabeth, M. Cyrklaff, W. Baumeister, and V. Lucic, J. Cell Biol. 188, 145 (2010).
- Yu et al. (2011) Y. Yu, H. Xin, R. Hovden, D. Wang, E. Rus, J. Mundy, D. Muller, and H. Abruña, Nano Letters 12, 4417 (2011).
- Sai et al. (2013) H. Sai, K. W. Tan, K. Hur, E. Asenath-Smith, R. Hovden, Y. Jiang, M. Riccio, D. Muller, V. Elser, L. Estroff, S. Gruner, and U. Wiesner, Science 341, 530 (2013).
- Ercius et al. (2009) P. Ercius, L. Gignac, and D. Muller, Microscopy and Microanalysis 15, 244 (2009).
- Jiang et al. (2017) Y. Jiang, E. Padgett, R. Hovden, and D. Muller, Ultramicroscopy 186, 94 (2017).
- Donoho (2006) D. Donoho, IEEE Trans. Inf. Thoer. 52, 1289 (2006).
- Candés et al. (2006) E. Candés, J. Romber, and T. Tao, IEEE Trans. Inf. Thoer. 52, 489 (2006).
- Saghi et al. (2011) Z. Saghi, D. Holland, R. Leary, A. Falqui, G. Bertoni, A. Sederman, L. Gladden, and P. Midgley, Nano Lett. 11, 4666 (2011).
- Bicer et al. (2015) T. Bicer, D. Gursoy, R. Kettimuthu, F. De Carlo, G. Agarwal, and I. Foster, Euro-Par 2015: Parallel Processing 9233, 289 (2015).
- Midgley and Weyland (2003) P. Midgley and M. Weyland, Ultramicroscopy 96, 413 (2003).
- Sidky et al. (2011) E. Sidky, Y. Duchin, and X. Pan, Med. Phys. 38, S117 (2011).
- Rudin et al. (1992) L. Rudin, S. Osher, and E. Fatemi, Physica D: Nonlinear Phenomena 60, 259 (1992).
- Sidky et al. (2006) E. Sidky, K. Chien-Min, and P. Xiaochuan, J.X-Ray Sci. Technol. 14, 119 (2006).
- Leary et al. (2013) R. Leary, Z. Saghi, P. Midgley, and D. Holland, Ultramicroscopy 131, 70 (2013).
- Guay et al. (2016) M. Guay, W. Czaja, M. Aronova, and R. Leapman, Scientific Reports 6, 1 (2016).
- Sidky and Pan (2008) E. Sidky and X. Pan, Pys. Med. Biol. 53, 4777 (2008).
- Liu et al. (2017) L. Liu, Y. Han, and M. Jin, PLoS One 12, 1 (2017).
- Padgett et al. (2017) E. Padgett, R. Hovden, J. DaSilva, B. Levin, J. Grazul, T. Hanrath, and D. Muller, Microscopy and Microanalysis 23, 1150 (2017).
- Banjak et al. (2018) H. Banjak, T. Grenier, T. Epicier, S. Koneti, L. Roiban, A. Gay, I. Magnin, F. Peyrin, and V. Maxim, Ultramicroscopy 189, 109 (2018).
- Elfving et al. (2014) T. Elfving, P. C. Hansen, and T. Nikazad, Inverse Problems 30, 1 (2014).
- Levin et al. (2016) B. Levin, E. Padgett, C.-C. Chen, M. Scott, R. Xu, W. Theis, Y. Yang, C. Ophus, H. Zhang, D.-H. Ha, D. Wang, Y. Yu, H. Abruña, R. Robinson, P. Ercius, L. Kourkoutis, J. Miao, D. Muller, and R. Hovden, Scientific Data 3, 1 (2016).
- Zhang et al. (2011) X. Zhang, J. Wang, and L. Xing, Med. Phys. 38, 701 (2011).
- Clarke et al. (1994) L. Clarke, I. Glendinning, and R. Hempel, Programming Environments for Massively Parallel Dirstributed Systems , 213 (1994).
- Dagum and Menon (1998) L. Dagum and R. Menon, IEEE Comput. Sci. Eng. 5, 46–55 (1998).
- Odstrcil et al. (2019) M. Odstrcil, M. Holler, J. Raabe, and M. Guizar-Sicairos, Optics Express 27, 36637 (2019).
- Ding et al. (2019) G. Ding, R. Zhang, and H. Xin, Scientific Reports 9, 1 (2019).
- Gordon et al. (1970) R. Gordon, R. Bender, and G. Herman, J. Theor. Biol. 29, 471 (1970).
- Kaczmarz (1993) S. Kaczmarz, Int. J. Control 57, 1269 (1993).
- Strohmer and Vershynin (2009) T. Strohmer and R. Vershynin, J. Fourier Anal. Appl. 15, 262 (2009).
- Sørensen and Hansen (2014) H. Sørensen and P. C. Hansen, SIAM J. Sci. Comput. 36, C524 (2014).
- Schwartz et al. (2019) J. Schwartz, Y. Jiang, Y. Wang, A. Aiello, P. Bhattacharya, H. Yuan, Z. Mi, N. Bassim, and R. Hovden, Microscopy and Microanalysis 25, 705 (2019).
Appendix A Supplemental Figures
Appendix B Dynamic Compressed Sensing Algorithm Parameters
The pseudo-code for our proposed algorithm is listed in Supplemental Algorithm 1. Its basic steps are summarized as follows: The outer loop controls the distribution of projections to the data matrix (). For each projection acquisition, the imaging data is appended to and utilized to further reduce the tomogram error (illustrated in Fig. 1a). The 3D reconstruction on lines 9 - 27 is similar to the pseudo-code in Sidky and Pan (2008) with some variation in the relaxation parameter reduction and POCS step. It is essential to allow the 3D reconstruction converge prior to appending more data to . Terminating the construction process will restrict the overall performance and cause convergence parameters to diverge from the true (optimal) solution. (illustrated in Fig. 8).
One can enforce data fidelity (POCS) by utilizing Kaczmarz algorithm, also known as the Algebraic Reconstruction Technique (ART). ART is an iterative method that solves inverse problems Gordon et al. (1970); Kaczmarz (1993) by sequentially using each row of the measurement matrix () to minimize the tomogram’s error. Each ART update follows:
where index is the iteration number and is a relaxation parameter. In this work, instead of sequentially cycling through the rows of ,we randomly select rows during ART update to achieve faster convergence Strohmer and Vershynin (2009) (demonstrated in Fig. 5). According to literature one can guarantee convergence when Sørensen and Hansen (2014). For the simulated and experimental objects studied here, we found optimal performance occurs when . If too large of a relaxation parameter is chosen, the ASD-POCS algorithm requires significantly more iterations to drive the data distance (DD) to reach (shown in Fig. 6).
The tomogram’s TV has a tendency to expand linearly with increasing projection number. As highlighted in the blue curve in Figure 8, the range between the minimal and maximal value grows consistently. To suppress this behavior, we introduced a envelope dampening calculation on line 8. This calculation ensures that when the relaxation parameter is reset, its magnitude decays linearly. Empirically we found that the ratio works consistently across all of our test data sets.
Outside of the influence from , there are additional parameters that control the dynamic CS algorithm. Complimentary to , there is reduction parameter that decays the step-size for ART. Similarly, there is also relaxation and reduction parameter for TV minimization (). should be less than so that the tomogram can maintain the data consistency constraint. The ratio of change in the 3D volume from the ART and TV stage is evaluated against (usually ). There also two variants total iteration number: minimum iteration to achieve convergence () and TV gradient descent (). Previous results show that 150 iterations is sufficient to achieve proper convergence Schwartz et al. (2019). Overall these parameters rarely need to tuned too drastically and the reconstruction quality will mostly depend on selection.