Determination of three-dimensional atomic positions from tomographic reconstruction using ensemble empirical mode decomposition

Tomographic reconstruction from a tilt series of electron micrographs has raised great interest in materials, chemical, and condensed matters science because of its capability of revealing 3D local atomic structures within nanomaterials. Previous breakthroughs have demonstrated that the positions of individual atoms can not only be visualized but also determined by combining a scanning transmission electron microscope with a high-angle annular dark-field detector, equally sloped tomography, and the filtering/denoising method. However, the filtering/denoising approach—whether imposed on 2D projections or 3D reconstruction—raises concerns regarding the robustness of image processing, the accuracy of atomic positions, and the artificial atoms introduced during filtering. In this article, we report a general method that overcomes these limitations. By removing unphysical oscillations in 2D projections through ensemble empirical mode decomposition and applying a standard Wiener filter to the 3D reconstruction, we are able to determine atomic structures with higher accuracy.


Introduction
Through considerable achievements in electron microscopy, researchers have been able to observe atomic structures in 2D images. Acquiring more high-quality projective images at different tilted angles is widely accepted as facilitating tomographic reconstruction for faithfully revealing the 3D distribution of atomic structures [1][2][3][4][5]. Recently, the 3D visualization of individual atoms has been realized by different groups, including in special conditions such as those in crystal structures and where sparse distribution is assumed [6][7][8][9][10][11][12][13].
One major problem in such visualization is that the signal-to-noise ratio of the acquired projections dominates the quality of the reconstruction. Poor reconstructions prohibit researchers from qualitatively extracting atomic positions. One solution is to assume crystalline structures as the given model and then generate a new atomic model according to few experimental projections. The advantage of this solution is obvious: with much less data acquisition, the rough 3D structures of radiation-sensitive samples can be obtained. However, if 3D local structures, such as defects and dislocations embedded within nanomaterials, are of interest, a tilted series of high-quality projections for 3D tomography is required.
To enhance the signal-to-noise ratio, denoising is required even when using a state-of-the-art instrument. For experimental projections with low signal-to-noise ratios, a typical approach is to use a filtering method [14]. However, the noise level must ususally be overestimated to obtain high-quality images [15]. When the noise is overestimated, a stronger filter is applied, thus generating artefacts such as atomic shifts and periodic atomic columns outside the particle, caused by the removal of a number of Fourier components [16,17]. Although Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence.
Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. those removed Fourier components with much smaller amplitudes are not likely to affect the appearance of large-scale dislocations, they cause ambiguity regarding the existence of point defects. Because the instrumentation is now sufficiently advanced for acquiring high-quality projections, performing gentle 2D denoising prior to or during 3D reconstruction has been reported to be feasible.
Several effective denoising techniques-from spatial to transform domain filtering-such as Wiener filtering [10,18], total variation minimization [19][20][21], and block-matching and 3D filtering [11,22], have been thoroughly developed for image processing and widely applied to tomographic reconstruction. Notably, denoising parameters are usually subjectively selected according to the visualization expected by individual researchers. The ensemble empirical mode decomposition (EEMD) method with physical understanding is a more general approach to addressing this problem.
Empirical mode decomposition (EMD) [23], based on Hilbert-Huang transform (HHT), is an adaptive signal-processing algorithm that decomposes a set of 1D data points into several intrinsic mode functions (IMFs). In contrast to Fourier-based transformation used in the reciprocal space, the decomposed IMFs generated by HHT in real space are not limited to linear and stationary data. Therefore, the obtained mixed data, including pixel-wise noise, subtle fluctuations, meaningful signals, and the background, could be decomposed in the order of different length scales. Although the mode mixing problem of EMD was reported previously [24], EEMD, in which white noise is imposed on multiple copies and modes are decomposed independently, has been proven to be an effective approach to the problem [25].
In consideration of the robustness of EEMD and the prior knowledge of the d-spacing of atomic structures, any IMFs with conformational oscillations lower than the minimal d-spacing or higher than the resolution limit can be further excluded. A projection, the collection of remaining modes, with a reasonable length scale is then regenerated while maintaining the density distribution of the original one.

Image analysis using ensemble empirical mode decomposition
In this section, we detail the procedures of the simulated projection-based EEMD method. An image is considered as a collection of column or row data, and EMD starts with the recognition of all extremes (i.e., local maxima and minima) within a column or row. Upper and lower enveloped lines are generated by connecting all the maximal and minimal points, respectively, by using cubic spline interpolation. The input data are then subtracted using a median envelope (i.e., the mean curve of upper and lower envelopes), thereby yielding new input data for the next sifting iteration. Such a sifting process goes on until certain stoppage criteria are met, at which point the first IMF is formed. Typical termination criteria require the difference of the number between the extremes and zero-crossing points in the computed IMF be zero or one, and that the median envelope equals zero at any point. Removing the first IMF component from the original data yields a residue that has a larger length scale, and the iteration continues until a specified IMF number is acquired or the residue is a monotonic function that cannot be sifted further.
A fivefold decahedral Pt nanoparticle with a size of approximately 7.3×7.0×4.5 nm 3 is designed for comparing the performance of different approaches. In addition to combining five tetrahedrons, we generate an edge dislocation embedded within the nanoparticle and several dents on the surface by removing several atoms. It is widely known that when a sample is sufficiently thin along the beam direction (i.e., the z-direction) and when exact zone-axis directions are avoided, a scanning transmission electron microscope (STEM) image with appropriate inner-outer angles of the annular dark-field detector can be considered a satisfactory approximation of linear projection [4,5,[26][27][28]. In such an instance, the contrast of the image is dominated by scattering from the 1s states. This effect can be approximately simulated by convolving a Gaussian function with the probe [29], whereas the standard deviation can be determined by comparing the projection with that calculated using the multisliced method [30]. Figure 1 shows the calculated projection at 0°. In figure 1(a), a model with 4,020 atoms is convolved with Coulomb potential and then the 3D potential is added up along the beam direction to generate a projection. We also performed a multisliced simulation according to the parameters of an FEI Titan 80-300 microscope (energy: 200 keV; spherical aberration: 1.2 mm; illumination semiangle: 10.7 mrad; HAADF detector inner angle: 35.2 mrad; HAADF detector outer angle: 212.3 mrad; pixel resolution: 0.352 Å) [10]. The simulated image is displayed in figure 1(b). To simulate the blurry effect caused by the probe, we convolve the projection with a Gaussian function. The standard deviation of this Gaussian function is determined by searching for the closest match with the atomic profile from the multisliced image. Figure 1(c) presents the convolved projection that satisfies the criteria, and figure 1(d) shows the line scans across the central row of atoms in figures 1(b) and (c). The agreement indicates that the linear approximation of the STEM image under the aforementioned circumstances is valid. Since the linearity can be optimized by adjusting the angles of the HAADF detector [31] and the experimental projections are usually more linear than simulated ones [9], this approximate model will be used to compare the performance of our approach with those of other methods.
To apply the Poisson noise, we have to convert the incident electron dose to the detected electron dose per pixel. In the experiment, the probe current is measured from the brightness of the phosphor screen. As a result, it can be treated as the incident electron dose on the specimen. The detected electron dose should be estimated as a fraction of the incident electron dose because the annual dark field detector only collects a small portion of incident electrons. A detected electron dose of´Å e 2.5 10 5 2 is used to generate 10% Poisson noise on the convolved projection, and this is our model. The level of Poisson noise is given by Now we can start applying 2D-EEMD to decompose the IMFs from the noisy projection. To avoid mode mixing, we add extra Gaussian white noise to the projection with an ensemble of 10,000 copies. The variance of the white noise is normalized to that of the input data. For each copy, 2D-EEMD is applied to decompose the noisy projection into three IMFs and the residue. After performing the same procedures for all 10,000 copies, we average each IMF and the residue. Figure 2 shows the final decomposition result. Unlike reciprocal space decomposition methods such as wavelet and Fourier transformation, EEMD decomposes signals in real space. Therefore, the IMFs in figures 2(a)-(c) correspond to small-to-large length scales, and figure 2(d), which has less obvious fluctuations, is the residue. Apparently, the second IMF contributes to most of the {111} and {200} lattice fringes. Because the resolution is limited by the probe size and spherical aberration, higher-resolution oscillations beyond the {200} lattice fringes are unlikely to be generated by true atomic structures; hence, the first IMF (i.e., the smaller length scale) can be removed because it is unrelated to the atomic structures of interest. The treated projection we use for the reconstruction is the add-up image of all other modes except the first IMF.
Notably, this approach is fundamentally different from conventional denoising techniques. Typical denoising methods are based on the noise model; they assume that noise is additive or multiplicative and attempt to identify mathematical or statistical means of separating the noise from the signal. This always raises concerns regarding whether the noise has been appropriately removed from the signal. In addition, the denoising parameters are subjectively selected; thus, denoising results from different individuals vary widely. In our approach, instead of performing denoising, we utilize physical understanding to remove nonstructural fluctuations from projections. Using a large ensemble with the EEMD method guarantees excellent robustness such that different researchers obtain the same results from the same data.

Compare tomographic reconstructions from different approaches
We repeat the aforementioned procedures to generate a tilt series of 63 projections with an angular range from −72.6°to 72.6°and equally sloped increments. All the projections are convolved with the same Gaussian function determined at 0°, and 10% Poisson noise is applied.
Next, we perform tomographic reconstruction using equally sloped tomography (EST) and then compare different denoising approaches. In the first approach, the EST reconstruction is obtained from raw projections alone without any treatments. A one-pixel thick central slice on the XZ plane is depicted in figure 3(a). The atomic fringes are vaguely visualized, and the white arrow indicates a dent on the surface. In the second approach, a standard 3D Wiener filter (i.e., λ=1) is applied to enhance image quality by minimizing the L2norm. The same slice is shown in figure 3(b) after filter application. Although the edge dislocation in the lower grain is clearly visible, a weak-density fake atom is generated in the dent region. In the third approach, to further enhance the signal-to-noise ratio, a stronger filter is usually applied to obtain a higher-quality image. In figure 3(c), the noise level is overestimated by a factor of two in the Wiener filter (i.e., λ=2). The image becomes much less noisy, yet the density of the fake atom becomes comparable to that of a real atom.
In our new approach, we apply 2D-EEMD to decompose the same 63 noisy projections and then remove the first IMF. After obtaining the 3D reconstruction, we perform standard 3D Wiener filtering. Figure 3(d) presents the same slice as that in figures 3(a)-(c), revealing that not only all the atoms are clearly visible, but no fake atoms are present.
To further quantitatively compare our method with other approaches, we perform atom tracing from 3D reconstructions. The atom-tracing procedure is as follows. (i) Interpolate twice on the voxels of the reconstruction (pixel resolution: 0.176 Å). (ii) Select a voxel as the central voxel of a spherical region with a ninepixel diameter (i.e., the peak size). A peak is considered present if the central voxel is the maximal value within its spherical region. (iii) After a peak is recognized, a spherical region with a 15-pixel diameter is then selected for calculating the centre of mass by The total density of a peak is the add-up value of all the voxels within the spherical region.
(v) Finally, verify every voxel in the interpolated 3D reconstruction, and repeat steps (ii)-(iv). Two considerations are carefully examined in the comparison. One is the accuracy of atomic positions and the other is the correctness of the number of atoms traced. The root-mean-square distance (RMSD) is a typical criterion for evaluating the accuracy of traced atomic coordinates given by where N represents the number of pairs of equivalent atoms between the tracing result and the model, and x i traced  and x i model  are the coordinates of the ith traced atom and of its corresponding atom in the model, respectively. In addition to the RMSD, we also determine the number of atoms according to the traced peaks. The 3D reconstruction contains both true atoms and density fluctuations; therefore, if any fluctuation meets the tracing algorithm criteria, then it is also recognized as a successful peak. A simple means of distinguishing true atoms from noisy peaks is to examine the total density of the peak. If the total densities are sorted from large to small, then the distribution is ideally similar to a step function, with values for true atoms and zeros for other peaks. When we take the derivative of the distribution, the index of the minimal derivative indicates the number of atoms.
According to the aforementioned ideas, we trace all the peaks in four 3D reconstructions: (i) the raw reconstruction, (ii) the raw reconstruction with a standard Wiener filter (λ=1), (iii) the raw reconstruction with a strong Wiener filter (λ=2), and (iv) the reconstruction after 2D-EEMD analysis with a standard Wiener filter (λ=1). Figures 4(a)-(d) show the distribution of cases (i-iv), with the traced peaks sorted from large-tosmall total density (in blue) and their corresponding derivatives (in green). Because of the noise and missing data, the distribution would not produce a perfect step function. The sharper the drop at the boundary is, the more precisely the position of the minimal derivative can be determined. Figures 4(a) and (b) appear to suggest that the correct numbers cannot be determined from the distributions because the drops are smooth around the 15,000th-18,000th and 5,000th-7,000th atoms, which are far from the correct number. After a stronger Wiener filter (i.e., λ=2) is applied, the number becomes reasonable. Figure 4(c) reveals that the number of atoms would be approximately 4,000-4,500; however, the smoothness around the edge increases the difficulty of determining the number precisely. After calculating the minimal derivative, the determined number is 4,301, which is still 7 % different than the correct number. Through our method, the curve drops almost vertically at the boundary of true atoms and noisy peaks. The distribution shown in figure 4(d) obtains the exact number of 4,020 atoms.
It is usually difficult to select a suitable number for the peak size when performing atom tracing. When the quality of reconstruction is poor, selecting a small peak size generates a large number of peaks, thus producing interstitial atoms. Nevertheless, if a large peak size is selected, then numerous peaks are excluded, thus generating vacancies. To demonstrate the robustness of our approach, we list both the RMSD values and the numbers of traced atoms with different peak sizes in table 1, from a diameter of 9 to 15 pixels. The large RMSD and the incorrect number of atoms indicate that neither correct atomic structures nor a true atomic profile are obtained. Applying a denoising technique clearly enhances the image quality, but the result is suboptimal when we use a standard Wiener filter. When a stronger Wiener filter is applied, the RMSD is significantly reduced, but the number of traced atoms remains higher than the correct number. The excess traced atoms do not truly exist in the specimen. This implies that artificial atoms are generated when a strong filter is imposed. Finally, our approach demonstrates the stability of the traced results from different peak sizes, indicating that removing unphysical oscillations from raw projections indeed reduces fake peaks in 3D reconstructions. In addition, applying a standard filter to the improved 3D reconstruction solves the fake-atom problem and even provides superior traced coordinate accuracy.

Apply ensemble empirical mode decomposition to experimental data
We apply the EEMD approach to experimental projections. Although no experimental 3D reconstruction is demonstrated in this paper, a specimen of AuPd bimetallic nanostructures grown using a sputter is employed to study the feasibility of our method. Figure 5(a) displays the raw image acquired using an FEI Titan at 200 keV. After the EEMD treatment, the first and second IMFs of the experimental image are removed because the third IMF corresponds to the main feature of lattice structures. The add-up image of the remaining IMFs and residue is shown in figure 5(b). Figure 5(c), the add-up image of the first and second IMFs, represents the difference between figures 5(a) and (b). The subtraction apparently removes unphysical oscillations from the real-lattice . The x and y axes indicate the index of the atoms and the corresponding total density, respectively. The y axis at the right indicates the strength of the derivative. Ideally, the distribution should form a step function with a sharp edge at the boundary of true and fake atoms. In reality, the curve is smooth. A sharper edge yields a more precise determination of the number of atoms. Table 1. Number of traced atoms (N atoms ) and the RMSD values from different reconstructions. The results obtained using different peak sizes for atom tracing are listed. The Wiener filter is labelled as WNF(λ), and D peak represents the peak diameter (in pixels). Our approach (i.e., EEMD + WNF (λ=1)) is the most accurate method in both the number and positions of traced atoms, even though the peak size is changed. structures. Close-ups of the insets (i.e., the yellow boxes) from top to bottom are provided in figure 5(d) for more detailed comparison. Notably, the mismatches of successive rows resulted from small errors of initial points during the scanning process are extracted from the original image, thus forming the pattern of dashed lines in figure 5(c). Clearly, those dashed lines are not true atomic structures and deteriorate the 3D reconstruction if they are not subtracted.

Conclusion
In this article, we first present the new idea of analysing acquired projections by removing unphysical oscillations of ultrashort-length scale rather than performing denoising artificially. We then employ the 2D-EEMD method to realize the decomposition of modes with different length scales. By removing some IMFs from the projection, those ultrashort-length oscillations no longer affect 3D reconstruction quality, thus making the atom tracing easier and more precise. In addition, we visually and quantitatively compare the results obtained using our method with those obtained using other methods. The results showed that our method is superior to the other methods regarding fake atoms, precise atomic position, and correct number of atoms. These crucial aspects are directly related to whether the 3D experimental reconstruction is sufficiently informative for the subsequent theoretical calculations. Finally, no parameters must be adjusted subjectively, thus guaranteeing consistent results.