mpfit: a robust method for fitting atomic resolution images with multiple Gaussian peaks

The standard technique for sub-pixel estimation of atom positions from atomic resolution scanning transmission electron microscopy images relies on fitting intensity maxima or minima with a two-dimensional Gaussian function. While this is a widespread method of measurement, it can be error prone in images with non-zero aberrations, strong intensity differences between adjacent atoms or in situations where the neighboring atom positions approach the resolution limit of the microscope. Here we demonstrate mpfit, an atom finding algorithm that iteratively calculates a series of overlapping two-dimensional Gaussian functions to fit the experimental dataset and then subsequently uses a subset of the calculated Gaussian functions to perform sub-pixel refinement of atom positions. Based on both simulated and experimental datasets presented in this work, this approach gives lower errors when compared to the commonly used single Gaussian peak fitting approach and demonstrates increased robustness over a wider range of experimental conditions.


Introduction
The development of spherical aberration-correction for Scanning Transmission Electron Microscopy (STEM) imaging has been one of the biggest triumphs of electron microscopy over the past several decades, allowing the subångström resolution imaging of crystal structures [1,2,3]. Several pioneering STEM experiments have demonstrated the feasibility of this technique for the direct visualization of atom positions from aberration-corrected STEM images and has proved itself an invaluable tool for sub-ångström resolution structural measurements [4,5,6,7,8]. While the typical aberration-corrected STEM electron beam has a probe diameter approximately between 0.5Å and 1Å, supersampling scanning positions below the Nyquist-Shannon sampling limit and the subsequent tting of the probe image with a two-dimensional Gaussian function allows the sub-pixel precision assignment of atom column positions from aberration corrected STEM datasets [5,9,10,11,12,13,14]. This technique has been used for quantitative atomic displacement measurements across thin lms, 2D crystals, domain boundaries and has allowed the experimental observation of novel structural phenomena such as polar vortices [15,16,17,18,19,20].
While the Gaussian function tting approach is an extraordinarily powerful technique, one noted shortcoming is that it assumes well separated atoms with no overlap, or negligible aberrations in the beam itself -conditions that are only available under a certain limited set of imaging conditions [16,17]. Typically, such an imaging setup uses a ring shaped annular detector with the outer and inner detector collection circles centered along the microscope optic axis. Such a con guration will have an inner collection angle of approximately 85-90 mrad to capture only the incoherently scattered electrons, and is conventionally referred to as High Angle Annular Dark Field STEM (HAADF-STEM) imaging [5,21]. This mode of imaging is referred to as dark eld imaging since the vacuum is dark while the atom columns themselves are bright due to electrons preferentially scattering from atomic nuclei. Incoherent electron scattering at high angles is a consequence of Rutherford scattering where the electron is scattered due to nuclear Coulombic forces [22,23]. Since the Coulombic force experienced by the electron probe is directly proportional to the number of protons in the nucleus (Z), atom column images in HAADF-STEM datasets generate peaks with an almost linear relationship of intensity ∝ Z 2 with the atomic number and is also referred to as Z contrast imaging [24,25,26].
Z-contrast imaging however is generally considered unsuitable for imaging lighter elements such as oxygen, boron or carbon [19,21,20]. However, structural metrology for  Fig.1(a) and Fig.1(b) with the blue arrows referring to the intensity minima and the green arrows referring to the single peak fits.
many scienti cally important material systems such as ferroelectrics needs the imaging and quanti cation of lighter atoms as well as heavier elements [27,28]. This problem can be signi cantly mitigated in Bright Field STEM (BF-STEM) imaging, where rather than annular detectors a circular detector is used with the detector center coinciding with the optic axis of the microscope [19,29]. The conventional collection angle ranges in BF-STEM imaging extend up to 15mrad, signi cantly lower than even the inner collection angle for HAADF-STEM [29]. Because unscattered electron beams are imaged by this technique, in contrast to HAADF-STEM vacuum is bright, while the atom positions have comparatively lower intensity. The ideal BF-STEM image would thus have an intensity pro le complementary to the images obtained from HAADF-STEM imaging. However, in reality owing to coherent scattering e ects due to the low collection angles, atom positions are more blurred from aberrations that are more prominent in BF-STEM images [30]. Additionally, since BF-STEM images capture both light and heavy atom positions the inter-atomic distances are substantially smaller. These e ects result in atom positions that are non-Gaussian in shape, and often have intensity overlaps and tails coming from their neighbors making position metrology challenging in BF-STEM images.

Fi ing Atom Positions with Gaussians
The best modern aberration corrected microscopes can generate electron probes that are free of aberrations up to 30mrad, which corresponds to beam diameters that are of the order of 0.5 Å, or 50pm at 200kV [8,10]. Super-sampling the beam by a factor of ve results in scan positions that are spaced approximately 10pm apart from each other. For HAADF-STEM images where oxygen atoms are not observed, inter-atomic distances from the low index zones are mostly of the order of 1.5 Å, allowing enough distance between atoms so that they are well separated and thus an atom position can be reasonably approximated with a two-dimensional Gaussian intensity pro le. Since the FWHM of this Gaussian is around 50-75pm, this allows the determination of the peak of the Gaussian intensity distribution with accuracies approaching 0.5pm [14,17]. It is this combination of aberration-corrected imaging and Gaussian peak tting that has enabled modern electron microscopy to reliably measure domain walls, grain boundaries, defects, and strain with single picometer precision, making STEM imaging so powerful.
However, this approach runs into problems when applied to BF-STEM imaging. In Fig. 1(a), we show a typical BF-STEM image of LiNbO 3 with 4.9 pm scanning pixel sizes. The bright regions correspond to vacuum, while the darker regions are the niobium and oxygen atom columns with the blue dots corresponding to the intensity minima. While the intensity minima can be used as an initial estimate of atom positions, the error in such a measurement is at least of the order of the pixel size, which is 5pm in our case. This makes the error of measurement in BF-STEM an order of magnitude worse than the best HAADF-STEM results. Fig.1(b) demonstrates the same section of the BF-STEM image with the re ned atom positions obtained from tting the intensity distribution with a single Gaussian peak with green dots next to the intensity minima (blue dots). A visual estimation shows that the tted Gaussians do not reliably converge on the atom positions, and are often tens of picometers away when the intensity minima is weak, and the neighboring atom is close. In some cases, the re ned atom position is in the middle of the two neighboring atom columns with no de nite atomic intensity.
This can be quantitatively demonstrated by pro ling the summed intensity distribution (Fig. 1(c)) from the region shown along the white arrows in Fig.1(a) and Fig.1(b). The blue arrows in Fig.1(c) correspond to the intensity minima, while the green arrows correspond to the Gaussian re ned atom positions. The presence of an intense neighboring atom's intensity tail gives rise to a dip in the intensity away from the original minima, right in the middle of two atom columns and the Gaussian peak tting technique converges to that local minima rather than the original position. Previous BF-STEM imaging has attempted in circumventing such issues by using a multi-parameter Gaussian peak, or performing image metrology through multivariate statistics rather than tting each individual atoms [19,20]. Both approaches require an initial knowledge of the crystal structure being imaged. Multi-parameter Gaussian ts need an estimation of the number and location of the nearest neighbors, and thus cannot be applied as a robust technique as it necessitates custom tting equations for individual crystal structures. In particular, this restriction limits the application of this method where more than one crystal orientation may be present. Here, we propose a novel multi-Gaussian re nement routinemp t -that does not require prior knowledge of the crystal structures being imaged and can robustly re ne a wider variety of images by deconvolution of a subsection of the image into multiple overlapping two-dimensional Gaussians. Since HAADF-STEM image renement requires less stringent conditions, our algorithm extends equally well to such systems too.

The mpfit algorithm
The Gaussian curve is a centrosymmetric curve with wide uses in single processing for approximating symmetric impulse functions [31,32]. Moreover, it has been demonstrated that given a su ciently large number of Gaussians, any non-in nite signal can be approximated as a sum of overlapping Gaussians [31,32]. We use this insight and extend it into two dimensions by rst modeling our observed atom intensity as a sum of overlapping Gaussians. The second key idea is to recognize that not all Gaussian functions that are approximating the region of interest are in fact originating from the atom whose position we are trying to re ne. Thus the Gaussian functions are subsequently sorted and only a subset of them that approximate the atom position are used to re ne the atom. The ow chart of our algorithm is illustrated in Fig.2. The steps of the mp t algorithm can be described as: 1 Get intensity minima/maxima: The initial starting point of this algorithm is the calculation of intensity maxima for inverted contrast BF-STEM images or ADF-STEM images. 2 Calculate median inter-neighbor distance: Following the identi cation of intensity minima, the median inter-peak distance is calculated. 3 Get region of interest: The region of interest is cuto as a square with the intensity minima as the central pixel, with the sides of the square given by = + , where s is the side of the square and is the nearest even number to the median inter-peak distance. Thus the + , + pixel in the square is the intensity minima that was the original starting point. 4 Fit iteratively with Gaussians: The region of interest is then t by a single 2D Gaussian function with a user determined tolerance factor. The tolerance factor refers to the mean absolute di erence in intensity between the tted gaussian and the original data. The tted Gaussian function is then subtracted from the original region of interest, and the residual is subsequently tted again. This process continues for a pre-determined number of iterations, with the sum of all the Gaussians then subsequently representing the original box. In the authors' experience, the tolerance factor is less important than the number of iterative Gaussians used, with reasonable accuracy and speed being obtained with a tolerance of 10 −12 and 12 to 16 iterations for the mp t example presented here. 5 Sort peaks and get the re ned position: The Gaussian peaks are then subsequently sorted based on their distance from the original minima with only peak positions whose distances are less than half the nearest neighbor distance used for re nement. The re ned atom position is then the amplitude normalized average peak position of all the selected Gaussians.

Results and Discussions
Results on Simulated BF-STEM Images The e ciency and accuracy of the mp t algorithm was tested on simulated BF-STEM images of LiNbO 3 . The advantage of simulated data is that the accurate atom positions are already known and can be compared with mp t results. This allows the estimation of the relative errors of the single Gaussian and the multiple Gaussian mp t approaches, with the simulation parameters outlined in Table 1. Following the steps of the algorithm outlined in Fig.2, the intensity minima were rst calculated for the simulated image, with Fig.3(a) demonstrating the simulated BF-STEM image of LiNbO 3 with the intensity minima overlaid as blue dots. These intensity minima are subsequently used to calculate the median nearest neighboring distance between the minima, ( ), which is 18 pixels when rounded to the nearest even number. Based on the calculated value, the region of interest for this image thus corresponds to a square of 19 × 19 pixels, which is demonstrated for one of the atoms as a red square in Fig.3(a). The region of interest for that atom is shown in Fig.3(b) with the contrast inverted and the intensity minima for the atom in question overlaid as a blue dot. As could be ascertained from Fig. 3(b), the intensity distribution from the bottom left atom partially overlaps with the atom position we are aiming to re ne, precisely indicating the scenario where single peak Gaussian tting approaches often give erroneous results.
Twelve iteration steps were chosen to represent this section of the image, as per step four of the algorithm. The iteration steps and the evolution of the Gaussian summation is shown in Fig.3(c). The calculation of the Gaussian is performed by taking in the entire image, and calculating a two-dimensional Gaussian peak with the smallest absolute di erence with the initial region of interest. Multiple di erent Gaussian tting approaches can be used, with the tting equation used in this approach expanded in Eq. 1. As could be observed from Fig.3(c) the summation of the Gaussian peaks starts to approximate the region of interest within only a few iterations. This demonstrates that the iteration number chosen was su cient enough to capture the complexities of the data being tted. It is even more interesting to look at the result of the rst iteration, which is mathematically equivalent to the single Gaussian peak tting approach. As the rst iteration in Fig. 3(c) shows, the single Gaussian peak tting approach is a special case of the mp t algorithm, where the number of iterations is one. According to this image, the rst Gaussian peak does not exist near the center of the image, and is extracted towards the bottom left corner. The central peak related to the atomic column of interest is captured in the second iteration rather than the rst, thus visually demonstrating why the single Gaussian approach fails in some cases.
While it may be possible to adjust the calculation of the region of interest to capture the atom position accurately, this approach necessitates tinkering with multiple di erent collection areas and a non-uniform solution for all the atoms in the image. mp t on the other hand, removes the necessity for such complicated user modi cations, allowing the estimation of all the Gaussian peaks that contribute to the nal image. The individual peak positions are also visually represented as a function of the iteration number in Fig.3(d) where the central black line corresponds to the known atom position. The white sphere, corresponding to the zeroth iteration, is the intensity minima, which is the starting point of the calculations. The colors of the spheres themselves are scaled to their relative contributions. The spheres with a red border represent Gaussian peaks that were assigned to the neighboring atoms. As could be observed, the peak obtained from the rst iteration is assigned to the neighboring atom, with the second iteration being used for the atom position calculation. Summing all the Gaussian functions together we obtain Fig.3(e), which shows close delity to the input data ( Fig.3(b)).
Based on the nal step of the algorithm, the Gaussian peaks are assigned either to neighboring atoms or the central atom depending upon the distance of the peak center from the initial intensity minima. As can be seen in Fig.3(e), the intensity minima is not always a reliable estimator of the actual atom position, but the mp t algorithm converges extraordinarily close to the actual atom position (green and red dots overlapping), demonstrating its superiority. The representative summation of the Gaussian summation can thus be broken down into two components -the Gaussian peaks that were used for atom position renement, the sum of which is visualized in Fig. 3(f) and the Gaussian peaks that were further o , and assigned as contributions of neighboring atom intensity tails -represented in Fig.3(g). Thus, the combination of the main atom and the neighboring contributions gives rise to the total intensity pro le that was observed.
We further evaluated the accuracy of the mp t algorithm for an entire image rather than a single atom. Fig.4 shows and compares the three di erent atom position metrology techniques -intensity minima/maxima, single Gaussian peak tting and the mp t algorithm with each other respectively. Fig.4(a) demonstrates the intensity minima itself may not be coincident with the ideal atom positions due to minute intensity variations that are not accurately captured given a limited detector dynamic range, with the errors of the order of a single pixel. As a result, the intensity minima atom positions are clustered at several di erent clustering values, which can be understood based on the fact that results from the intensity minima are always on the order of a pixel. Thus compared to position re nement algorithms, just the minima itself is incapable of sub-pixel precision metrology. Fig.4(b) demonstrates the di erence of the single peak approach from the ideal atom positions, with the results being clustered into three distinct lobes. This can be understood based on the fact that there are three separate types of intensity distributions in the simulated data. For well-separated atoms, the single peak and the atom positions show close agreement, which generates the central lobes. However, there are also atom columns, where the neighboring atoms are either on the top left or the top right, giving rise to the two extra lobes -demonstrating the shortcomings of this approach when the intensity distributions of neighboring atoms approach the resolution limit of the electron microscope. The results from the mp t algorithm, demonstrated in Fig.4(c), are on the other hand clustered in a region less than 0.5pm across from the known atom positions -demonstrating it's accuracy. However in the authors' experience, the mp t technique fails to converge for the edge atoms, which shows up as atom positions that are not clustered and have a higher error. For the rest of the atoms in the image, however, mp t is signi cantly superior to the other approaches.

Results on Experimental BF-STEM Images
Along with simulated datasets, we additionally performed position metrology on experimental BF-STEM im- ages of LiNbO 3 viewed from the 1100 zone axis [15]. The results were obtained through STEM imaging in a spherical aberration corrected FEI Titan 3 transmission electron microscope, corrected for upto third order spherical aberrations. Imaging was performed at a camera length of 145mm and the BF-STEM images were collected using Gatan detectors with an outer collection semi-angle of 15mrad, using scanning pixel step sizes of 9.8pm. In contrast to simulated datasets, the exact ideal atom positions are not known owing to specimen drift, thermal vibrations, signal to noise ratio, and localized imperfections in the crystal lattice. Fig.5(a) demonstrates a region of interest in an experimental dataset, with the intensity reversed, with Fig.5(b) showing a section of the image marked by the red box in Fig.5(a). The blue dot in Fig.5(b) in corresponds to the intensity minima, while the green dot represents the position calculated by the single Gaussian peak tting approach. As can be visually ascertained, the calculated atom position does not correspond to the atom position, and thus is an inaccurate representation. Following step 4 of the algorithm, and similar to the procedure outlined in Fig. 3(c), the region of interest is represented by a succession of closely spaced two-dimensional Gaussian peaks over twelve iteration steps,with the contribution from the steps shown in Fig.5(c).
The individual Gaussian peaks that contribute to thenal representation of the region of interest are pictorially represented in Fig.5(d), with the zeroth iteration peak in white corresponding to the intensity minima, which is the starting re nement step. The higher order iteration peaks are colored based on the absolute magnitude of their contribution to the nal representation. As can be seen it is the rst 2-3 peaks that have the strongest contribution, demonstrating that twelve iteration steps are su cient. Peaks that are further from the original intensity minima by more than twice the median inter-peak distance are indicated by a red border and are assigned to the intensity tails from neighboring atoms and are not assigned to the main atom that is being re ned. respondence between the experimental and represented data. Extending the number of iterations would allow progressively smaller Gaussian peaks resulting in better correspondence, but would also increase the demand for computational resources without a correspondingly signi cant increase in precision. The intensity minima are overlaid on the images in blue, with the results from the single peak t approach in green and the mp t results in yellow respectively. Thus the mp t algorithm accurately determines the atom location rather than converging to saddle points created from intensity tails from neighboring atoms. Fig.5(f) represents the sum of the Gaussians that represents the atom being re ned and Fig.5(g) represents the contribution from the intensity tails from the neighboring atoms, and is calculated from the Gaussian peaks represented with red borders in Fig.5(d).
Returning back to the original experimental BF-STEM image in Fig.1(a), we revisit that experimental data in Fig.6(a), comparing the results obtained with the mp t approach. As could be visually ascertained, while the single peak t approach fails in some of the cases, the mp t approach reliably re nes to the atom position, which can also be ascertained by the intensity pro le demonstrated in Fig.6(b).

Comparisons with other algorithms
Several other specialized algorithms have been designed to quantify atom positions in electron microscope datasets, such as Atomap [33], StatSTEM [34] and oxygen octahedra picker [35]. The Atomap algorithm uses principal component analysis to obtain denoised STEM images and nds the center of mass based on the initial guess of local intensity maxima or minima. Using the center of mass as the starting estimate, it then subsequently approximates a twodimensional gaussian to locate the estimated position of atoms. Atomap can additionally sort the di erent species of atom columns in the image and analyze them individu- ally. StatSTEM on the other hand is a model-based tting algorithm for extraction of the atom position information from STEM images. StatSTEM models the atoms in the images as the superposition of two-dimensional gaussian peaks, and since this is a model based technique it requires prior knowledge of the crystal structure of the sample being imaged to give a better estimation of the initial guess. After obtaining the initial guess, the algorithm will go through iterations to reach the least squares estimation of tting parameters, and then determines the position. oxygen octahedra picker is a software specialized in identifying the octahedra rotations in the ABO 3 perovskite oxides. It sorts out the oxygen and B atom positions and provides users the option of selecting a fast center of mass estimation or a slower peak tting with two-dimensional gaussians. It exhibits an impressive accuracy of as small as 3 pm in simulated HAADF images. However, the existing methods still possess limitations in practical cases -with neither the oxygen octahedra picker and the Atomap software being able to process STEM images where atomic columns being measured will have intensity contributions from their neighbors. Thus both these approaches work well for wellseparated atom columns in HAADF images, but face accuracy penalties with BF-STEM images. While StatSTEM's model based algorithm is able to solve the overlapping issue by assuming atom columns as overlapping 2D Gaussian peaks, its iterative model tting process is computationally intensive, and requires prior knowledge of the crystal structure being imaged. As demonstrated in Fig. 7(a), visually there is almost no di erence between the tting results of StatSTEM versus mp t, with StatSTEM's results being slightly o -centered from mp t's estimation. Comparison of the results in Fig.7(b) demonstrates that both technique give results that are less than a pixel apart from each other, with mp t outperforming StatSTEM. The standard deviation ( ) of mp t's estimation from known atom positions is 1.49pm compared to a of 3.31pm from StatSTEM.

Conclusions
While it may be possible to assume from the results presented here that the single Gaussian peak tting approach fails to converge to atom solutions and gives erroneous results, it actually performs perfectly adequately for the majority of STEM experiments. However, for certain nonideal imaging conditions, the single Gaussian peak tting approach fails, while mp t accurately obtains precise atom positions. For well-separated atoms, the results from mpt and a single Gaussian re nement are in fact identical. Additionally, it has to be kept in mind, that even with parallelization implemented, the mp t algorithm solves for over ten Gaussian peaks in a batch process. On the other hand, the single Gaussian approach solves for just one peak, thus making the single Gaussian approach faster by at around an order of magnitude.
Future planned improvements include solving for neighboring peaks simultaneously using the tail functions to deconvolve the full obtained image as an independent set of impulse functions originating from individual atoms. Additionally, atom columns whose separation distances are below the resolution limit of the microscope may be particularly suited for this approach, by the deconvolution of the observed impulse function into two closely separated Gaussians and enabling the super-resolution metrology of atom positions from STEM datasets.
Thus, our results demonstrate that the mp t algorithm can reliably and robustly re ne the sub-pixel precision of atoms even without a priori knowledge of the underlying crystal structure. Additionally, since the single Gaussian approach is a special case of the mp t approach with the total number of iterations as one, this approach will also work for ADF-STEM images, enabling a single approach to the metrology of a wide variety of STEM data. The results are superior to existing algorithms, and exceeds the state of the art -StatSTEM in accuracy, with the added advantage of being agnostic to the crystal structure being imaged. where ( , ) is the Gaussian output as a function of and , and are the two normal distributions in the and directions, 0 and 0 are the position of the Gaussian peak, is the amplitude of the Gaussian peak and is the rotation in the counter-clockwise direction of the two-dimensional Gaussian peak. Thus given a set of , , and values from the experimental region of interest, a Gaussian curve is estimated from Eq.1 such that: where is the tolerance, which was 10 −8 for our implementation. The Equation itself is calculated through the least-squares approach using the trust region reflective algorithm. Trust-region algorithms are an evolution of Levenberg-Marquardt (LM) algorithms. However, compared to the LM algorithms, this algorithm is curvature independent and is thus computationally significantly faster [36,37,38].

Simulation Parameters
The LiNbO 3 images were simulated using the MacTempasX so ware, with the simulation parameters enumerated in Table 1 [39]. Author details