Materials property mapping from atomic scale imaging via machine learning based sub-pixel processing

Direct visualization of the atomic structure in scanning transmission electron microscopy has led to a comprehensive understanding of the structure-property relationship. However, a reliable characterization of the structural transition on a picometric scale is still challenging because of the limited spatial resolution and noise. Here, we demonstrate that the primary segmentation of atomic signals from background, succeeded by a denoising process, enables structural analysis in a sub-pixel accuracy. Poisson noise is eliminated using the block matching and three-dimensional filtering with Anscombe transformation, and remnant noise is removed via morphological filtering, which results in an increase of peak signal-to-noise ratio from 7 to 11 dB. Extracting the centroids of atomic columns segmented via K-means clustering, an unsupervised method for robust thresholding, achieves an average error of less than 0.7 pixel, which corresponds to 4.6 pm. This study will contribute to a profound understanding of the local structural dynamics in crystal structures.


INTRODUCTION
Over the past decade, the experimental resolution of scanning transmission electron microscopy (STEM) has been improved using the spherical aberration corrector 1 . This success has enabled the formation of electron beam probes at a size of 1 Å, and the observation of local crystallographic structures at the atomic scale. Beyond an ordinary visualization of atomic structures, quantitative analyses of atomic positions, such as strain [2][3][4] , polarization [5][6][7][8][9] , and oxygen octahedral tilt [10][11][12][13] , have garnered an increasing interest. In general, structural transitions occur at several picometers of shifts in the atomic positions. Accordingly, the precise acquisition of atomic positions is significant in the quantitative analysis of STEM images. The demand for picometer-scale accuracy in atomic positioning motivated us to develop PRESTem, a program for analyzing atomic positions of STEM images at the sub-pixel level.
In the STEM images, atomic positions are deduced from the intensity values enclosing the atomic columns. Thus, noise in the intensity level causes inaccurate acquisition of atomic positions; hence, it should be eliminated as the primary factor. In recent years, considerable signal-to-noise ratio (SNR) enhancement has been reported in the summation of multiple images via image registration [14][15][16][17] , correcting the positional mismatches between the serially obtained images, which is typically called image drifts. However, in the expense of high signal-to-noise ratio, a summation of multiple images inevitably increases the total electron dose, and it adversely affect the specimen condition by causing structural transition or contamination, which result in degradation of data reliability; hence, a single image enhancement can be a general application for analyzing atomic positions. For single image enhancement, optimization of the Weiner filter and band-pass filter in the frequency domain yielded a considerable noise reduction by exploiting the periodic nature of the crystal structure in atomic scale images [18][19][20] . In addition, advanced noise filtering techniques, such as principal component analysis (PCA) 15 , nonlocal means filtering 21,22 , and block matching and three dimensional (BM3D) filtering 23,24 have been gradually introduced to the field of electron microscopy from the classical image processing field. However, these techniques are generally effective for Gaussian noise, whereas Poisson noise is the principal noise in STEM images.
Advanced computational techniques achieved precise quantification in atomic positions beyond the spatial resolution limit of STEM. The model-based estimation method is the most widely used approach, which adapts the pixel intensities of an atomic column to a specific model function and defines the atomic position where maximal intensity exists 2,5,16,[25][26][27][28][29][30] . Simple twodimensional (2D) quadratic function 30 and 2D Gaussian function 2,5,16,31-33 have been introduced as the model functions. The 2D Gaussian function has been continuously modified with empirical parameters to enhance the accuracy in detecting atomic positions, however, resulted in increased computational burden [25][26][27][28][29]31 . On the other hand, the center of mass method, a conventional approach to locate multiple particles in experimental images, still has a potential to achieve the sub-pixel accuracy, with low computational costs. The fundamental idea of the center of mass method is the summation of the intensities and pixel positions belonging to the particles. An optimal use of the center of mass method can yield a precision of 0.017 pixels with a low computational burden 34 . The simplicity of the method has enabled its easy application in locating atomic columns in STEM images 10,[27][28][29] . However, in most applications, the pixels to be summed using the center of mass method are approximately estimated from the lattice spacings and atomic column radii. This rough selection of estimating region inevitably includes background signal outside the atomic columns, resulting in a degraded accuracy. Thus, a clear separation of atoms from the background was set to be the starting point in the application of the center of mass method to detect atomic positions in STEM images.
Several software tools with unique combinations of denoising techniques and algorithms for determining atomic locations have been presented [28][29][30] . Peak pairs analysis (PPA) 30 , originally developed for phase contrast images in conventional TEM, has been extensively used for STEM images 25 . As a pathfinder, PPA employs simple algorithms, such as the selected inverse Fourier transformation for noise reduction and 2D quadratic fitting for locating atomic columns, providing limited information on the crystal structure. CalAtom provides multiple choices for denoising and locating processes 28 . However, CalAtom focuses on the intensity of each heavy element in 2D materials rather than their positions. Consequently, it is not optimized to extract information from atomic peaks with weak signals, such as oxygen. Atomap is represented by a 3-step serial fitting process for independent chemical species 29 . Atomap initially finds atomic positions using a local maximum greater than the intensity threshold from a denoised image with 1D PCA, and then iteratively refines the positions using the center of mass and 2D Gaussian methods. However, this procedure significantly increases the computational cost. One of the main contributions of this study is the evaluation of the precision of atomic positions extracted using the proposed method and other reported approaches.
Here, we propose a strategy for the precise analysis of atomic structures based on the image enhancement and detection algorithm optimized for atomic scale STEM images. Moreover, image enhancement is achieved by removing the predominant Poisson noise from the combination of Anscombe transform and BM3D. Additionally, morphological filtering and histogram stretching were applied to enhance the weak signals of light elements. For the precise detection of atomic positions, the binarization of atomic columns and background is introduced, which effectively defines the fitting region for the center of mass method. By using K-means clustering, an unsupervised learning method, this binarization can be achieved without any user intervention. The remainder of this paper is organized as follows. First, we describe the proposed noise reduction process for image enhancement and atomic segmentation. Next, we evaluate and compare the reliability of the acquired atomic positions using the PRESTem and other available software tools. Subsequently, the structural analysis of the perovskite oxide is presented. Finally, we draw the conclusions and describe the outlook for the PRESTem.

RESULTS AND DISCUSSION
For structural analysis of STEM images, the peak position of each atom should be accurately determined. However, STEM images are difficult to analyze owing to their low resolution caused by noises during the acquisition process. To guarantee the analytical accuracy of STEM imaging, we propose the sequential steps as follows: (i) pre-processing for noise reduction, and (ii) segmented searching for atom positioning. To this end, both high-angle annular dark field (HAADF) and annular bright field (ABF) images were effectively utilized to segment the entire atoms; thus, the positions of the heavy and light atoms are detected from the HAADF and ABF images, respectively. Figure 1 shows a flowchart of the proposed methodology. In the pre-processing step, the noise in the STEM image was removed using BM3D. Morphological filtering was subsequently performed to remove the background stains. In addition, histogram stretching was performed to increase the contrast between the atoms and background. In the segmentation step, K-means clustering was employed, which separates atoms and background as a binarized form. The coordinates of each atom were obtained from the center of mass in the binarized image. During the entire process, we first obtained the coordinates of heavy atoms in the HAADF image using the above process followed by the coordinates of light atoms in the ABF, except those found in the HAADF image (more details in Supplementary Note 1).

Noise reduction
Noises in STEM images can be modeled using Poisson and Gaussian distributions. Unlike Gaussian noise, which appears independently of images, Poisson noise is difficult to handle because it appears in an image-dependent form. However, owing to the characteristics of STEM images, the Poisson component cannot be excluded 35 . Hence, we applied a mixed Poisson-Gaussian model for noise modeling. Poisson-Gaussian noise modeling is defined as follows: where x ∈ X is the pixel position in the image domain X, z is the observed signal, y is the original noise-free signal, σ is a function that represents the standard deviation of the entire noise component, and g represents an independent random noise with a zero mean and standard deviation of 1 36 . The noise elimination process consists of three steps: noise estimation, Anscombe transformation, and BM3D. The first step is the noise estimation of images, which estimates the variance of Gaussian and Poisson noises 36 . Noise estimation is a method of dividing an image into patches to locally estimate the standard deviations of the noise components in each patch and fit the estimated values of noise components in the entire image to equalize the noise distribution in each area of the image. We investigated the effects of exposure time on the acquisition process of experimental HAADF-STEM images (Supplementary Note 2). Our results confirmed that the standard deviations of Gaussian noise ranged from 0 to 13.54, whereas those of Poisson noise ranged from 41.70 to 67.63, which verifies the dominance of Poisson noise in STEM images.
Anscombe transformation step 37,38 stabilizes the variance of Poisson noise and makes it a constant value. This process can transform Poisson noise into Gaussian noise and eliminate it using the Gaussian noise removal method. The last step is BM3D, which is an effective denoising method for Gaussian noise 23 . BM3D consists of the two iterative processes, including the grouping and collaborative filtering, as shown in Fig. 2. In the initial grouping stage, the noisy image is partitioned into several fragments and grouped into blocks with similar ones, as shown by the blue boxes in Fig. 2. The image fragments are matched by the similarity index of each fragment (L2-distance). In the collaborative filtering stage, each block is linearly transformed into a 3D form, and the noise is attenuated using a hard threshold. The primary denoised image blocks are restored to the original 2D form and position. As initial collaborative filtering is completed, a second grouping process is conducted on the filtered image, achieving an enhanced block matching, as shown by the red boxes in Fig. 2. A final estimate of the true signals with minimized noise is generated owing to the additional collaborative filtering using the Wiener filtering on the image blocks.
After the denoising process, blur and stains with bias remain in the background. In addition, strengthening of atomic contrast is required for precise segmentation. To this end, background correction and histogram stretching are performed. Morphological filtering is an effective method to remove stains and bias in the background 39 . In gray scale image, morphological dilation is the process of adding the structural operator values to the image intensity value and then selecting the maximum value among them to perform the calculation. As a result of morphological dilation, the bright part of the image increases. Conversely, morphological erosion is a method that selects and calculates the minimum value among the difference between the intensity value of the image and the structural operator, and the dark part of the image increases. Morphological opening, a combination of these two methods, is used to remove stains and bias from the background. In the morphological opening process, the image is first eroded and subsequently dilated to estimate the trend in the background signal. By subtracting the estimated background from the denoised image, small noise and stains in the background are eliminated (Supplementary Note 3). Next, histogram stretching is performed to emphasize only the peak part of atoms. Unlike histogram equalization, histogram stretching is more appropriate for the segmentation of background and atoms because it can emphasize the desired intensity part.

Finding the atomic positions
The conventional methods [25][26][27][28][29]31 have implemented segmentation using the Gaussian fittings. However, because the actual atoms are not always Gaussian, we applied the K-means clustering to binarize the background and atoms in the proposed method. Figure 3 shows the entire segmentation process. Figure 3a demonstrates the simulated STEM image, which reveals that the target images are not in a perfect Gaussian form.
K-means clustering divides the data into clusters of specified k based on the intensity levels of the images 40 . In this method, the data are divided by minimizing the variance of the distance in each cluster through iteration. As K-means clustering is an unsupervised method, when only k clusters are set, other parameters are not needed because the method does not require the ground truth. Thus, we performed K-means clustering to segment the background and atoms of the STEM image (Fig. 3b).
To find the correct atomic coordinates after K-means clustering, the center of mass of the individual binarized column was set as an atom position. Figure 3c presents the results of applying the Kmeans clustering and center of mass method to determine atomic positions in the entire STEM image. After these peak detection process, all atoms are classified by inputting unit cell information through our program (see Supplementary Note 4).

Performance depending on the noise level
The performance of our proposed noise reduction method was evaluated by varying the noise parameters and comparing the obtained results with those of the conventional denoising techniques 18,19,23 . Test images for the evaluation were generated from the simulated HAADF-STEM image of SrTiO 3 , with an incremental Poisson noise σ P in the range of 20-80, fixed Gaussian noise of σ G = 10, and additional background noise equivalent to 10% of the maximum intensity. Except for image noise, we excluded other nonideal simulation parameters influencing the atomic positions in the STEM image, such as sample mistilt, aberrations, thermal vibrations, and defocus. Detailed information for simulation condition is in "Methods." The criterion for the performance is set to be the peak signal-to-noise ratio (PSNR) calculated using the mean squared error between the true and   Figure 4a shows that our proposed method exhibits a visually stable denoising performance regardless of the input noise level, which corresponds to the stable PSNR values presented in Fig. 4b (comparisons with actual denoised results in Supplementary Note 5). Our noise reduction method, represented by the Anscombe transform, BM3D, and morphological filtering enhances the PSNR by approximately 8-12 dB compared to the noisy input image. It is noted that our proposed method shows better performance compared to other denoising techniques, including the BM3Donly, at every noise level. This is because these denoising techniques cannot efficiently remove the background noise alone; however, they significantly affect the image quality (Supplementary Note 6).
Additionally, we investigated the reliability of the atomic positions depending on the input noise level. Figure 4c shows the average errors of the atomic positions acquired from the denoised HAADF-STEM images for each noise parameter. As the noise level increases, the accuracy of the peak positions deteriorates, whereas PSNR of the denoised images exhibits a slight change. However, our proposed method exhibits a sub-pixel precision under severe noise conditions.

Benchmark
We benchmarked our proposed method against other reported methods 28-30 by comparing the precision of the estimated atomic positions. Simulated HAADF-STEM images with identified atomic positions are employed again for the benchmarks. Here, Poisson-Gaussian noise was applied with a Poisson noise of σ P = 40 and Gaussian noise of σ G = 10, which are estimated using the experimental HAADF-STEM images and noise estimation algorithm 36 . Figure 5a shows a partial area of the final test image with an artificial noise. The Sr atomic columns exhibit a brighter contrast compared to the Ti atomic columns. Oxygens do not appear in the HAADF image because of their small scattering angles. Ground truth, the actual atomic positions of the simulated images, were derived from the fractional coordinates of the input crystal structure and pixel scale of the simulated image. The atomic positions of the Sr and Ti sites were extracted using PRESTem from the test image and compared with the ground truth. Figure 5b shows the distribution of errors in both the x-and y-axis. The average error and standard deviation are 0.71 and 0.51 pixels, respectively, which implies the achievement of sub-pixel precision. The average error and standard deviation were calculated using the absolute values of errors. Histograms in two side panels of Fig. 5b count the error points at intervals of 0.1 pixels along both axes. These histograms demonstrate a normal distribution centered at the zero-error point, revealing no biases in the acquired atomic positions.
To compare the obtained results with those of the other STEM analysis programs, such as Peak pairs analysis (PPA) 30 , CalAtom 28 , and Atomap 29 , the same evaluation tests were performed. These three programs were selected as they employ different denoising techniques and locating algorithms for atomic positions. All processes of atomic positioning were followed using the available manuals, except for PPA that we changed the denoising technique to that of PRESTem. This was because the Bragg filter, the denoising technique in PPA, generated false atomic intensities in the blank space. The information for computation time is given in Supplementary Note 6, which includes the dependency of computational time on the number of atomic columns in the STEM images. Figure 5c presents the error distribution of PPA. The average error and standard deviation are 1.88 and 1.17 pixels, respectively. The error points were randomly spread, and no bias was observed. CalAtom yields an average accuracy and a standard deviation of 1.82 and 1.54 pixels, respectively, using a Gaussian noise filter and multiple-ellipse fitting (MEF) method 41 . Although CalAtom exhibits normally distributed error points near the origin in Fig. 5d, it has a similar average error and larger standard deviation compared to PPA because of its larger maximum error. Atomap, combined with 1D PCA denoising and multistep fitting process represented by 2D elliptical Gaussian, presents an average accuracy of 0.83 pixels and a standard deviation of 0.35 pixels. These values were comparable to or even better than those of PRESTem. However, the error points are biased along the y-axis in Fig. 5e, and the iterative process excessively increases the computational time one to two orders of magnitude larger than that of PRESTem ( Supplementary Fig. 5).
The center of mass method used in PRESTem is one of the fastest algorithms for locating particles in experimental images. However, its accuracy is often undervalued compared to modelbased estimation methods. Dedicated denoising techniques for Poisson noise and accurate selection of the fitting region can increase the accuracy compared to model-based estimation methods, while maintaining its fast speed.

Demonstration
Based on the atomic coordinates, the local structural information of materials can be derived at the unit cell scale, the so-called peak analysis. In practice, the accuracy of structural information obtained using the peak analysis should be confirmed in addition to the precision of the acquired atomic coordinates. In this study, we adopted the ABO 3 perovskite system as a pilot target in PRESTem and implemented an automated peak analysis procedure (Supplementary Note 7). Perovskite oxide systems have been widely studied owing to their valuable functional properties, including superconductivity, ferroelectricity, piezo-, and metal-insulator transition. Despite the simple crystallographic structure of ABO 3 , the functional properties of perovskite oxides arise from various structural transitions, such as, polar distortion, lattice elongation, and lattice shear. Thus, investigation should be conducted on these structural changes in perovskite oxides to unveil the relationship between the functional property and its crystallographic structure. In this context, we evaluated the accuracy of the quantification by estimating the atomic displacement and the ratio between the lattice parameters c and a (c over a ratio) of the perovskite structure.
As for the atomic displacement, both the magnitude and angle of the displacement vector are important for the classification of the phase and estimation of polarity. Therefore, we defined an arbitrary system with SrTiO 3 , mimicking the polar vortex discovered in the perovskite oxide superlattice 42 . In this arbitrary structure, each Ti atom is displaced from the centroid of neighboring Sr atoms, while displacements are arranged in a vortex shape (Fig. 6a, b). The HAADF-STEM simulation was conducted along the [100] zone axis, and the output image was corrupted with a Poisson noise of σ P = 40 and Gaussian noise of σ G = 10, the same parameters used in benchmarks. Displacements of the Ti atoms were evaluated using PRESTem and compared with the true values provided from the simulated structure. The average and standard deviation of the error were 1.57 pm/4.25°a nd 1.13 pm/3.07°, respectively, calculated from 256 unit cells. A detailed comparison of the selected 8 unit cells is presented in Fig.  6c. Consequently, the proposed method guarantees an error of <3 pm. In addition, the tendency of a gradual decrease in the displacement angle was maintained in the evaluated results despite some errors.
To demonstrate the c over a ratio evaluation in the heterostructure, HAADF-STEM simulations were conducted on the (SrTiO 3 ) 10 /(PbTiO 3 ) 10 superlattice (Fig. 6d). In this structure, the lattice parameter c in PbTiO 3 is compressed and gradually changes between the interface 43 . This image was corrupted with the same Poisson-Gaussian noise, and subsequently processed with PRE-STem, which is comparable to the displacement test above. The evaluated c over a ratio in the PbTiO 3 layers was in the range of 1.00-1.05, which is in accordance with the truth value, as shown in Fig. 6e. Remarkably, the standard deviation was 0.014 in the SrTiO 3 region compared to 0.004 in the PbTiO 3 region. This is because the contrast of Sr diminishes owing to the strong signal of Pb in the atomic number (Z) contrast imaging of HAADF-STEM.
In conclusion, we proposed a powerful system for the quantitative analysis of local crystal structures from atomic-scale STEM images. The BM3D, morphological filtering, and histogram stretching processes were applied consecutively that effectively removed the STEM image noise. Segmentation using the unsupervised training method, K-means clustering, robustly separated the atomic signals from the background plane. Centroids of each atomic column were extracted and classified by user-input periodicity. We developed the PRESTem program using these methods with an automated structure analysis.
The accuracy of the atomic positions obtained using the PRESTem program was evaluated and compared with that of the existing programs represented by 2D quadratic fitting, MEF optimization, and 2D Gaussian fitting. With the binary segmentation process of K-means clustering, the commonly undervalued centroid method in PRESTem exhibited errors of less than 1 pixel, which was comparable to that of the 2D Gaussian fitting without high computational costs and biased errors.
Recent progress in STEM imaging has enabled rapid acquisition of atomic images with valuable information. The development of a precise and automated analytical method can bridge the gap between the two processes of image acquisition and image analysis. It is expected that our program will provide an opportunity for the statistical investigation of the material system overcoming the technical limitations of STEM, which cannot represent the overall characteristics of the specimen. In addition, our proposed methods for image enhancement and locating atomic positions can be robustly applied to atomic-scale STEM images of any other crystallographic structures, while structural quantification is preliminarily focused on perovskite.

STEM simulations
Except for the demonstration part, the input crystal structure for the simulation was built from cubic Pm3m SrTiO 3 , but Ti positions in each unit cell were arbitrarily modified to diversify the atomic positions. Within 20 × 20 × 1 cells of SrTiO 3 , all Ti atoms had a deviation of 20 pm from the body-centered position, along a random direction under the (100) plane. In all crystal structure for the STEM simulation, the decimal point of fractional coordinates was considered up to 6. STEM simulation was conducted by the multi-slice method using Dr. Probe software 44 . Sampling rate was set to be 6.6 pm per pixel, which approximately accords to the experimental imaging condition of 15 million times magnification, and the thickness of specimen was set to be approximately 5 nm by duplicating the input crystal structure along the projection axis. Detailed simulation parameters are in Table 1.

DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.

CODE AVAILABILITY
The codes used to calculate the results of this study are available from the corresponding author upon reasonable request. Also, our program is publicly available on GitHub (https://github.com/Cheque93/PRESTem).