Introduction to Blob Detection (BD) Algorithm

We have employed blob detection (BD) algorithms (Lindeberg, Reference Lindeberg1993), specifically, those available in the Python package skimage (van der Walt et al., Reference van der Walt, Schönberger, Nunez-Iglesias, Boulogne, Warner, Yager, Gouillart and Yu2014). We find that standard BD works to locate the rough positions of atomic columns, but to further improve precision we introduce additional steps, illustrated in Figure 1, to refine the position and column intensity estimation. Indeed, blob detection routines, which are popular and well-developed methods in the field of computer vision (Lowe, Reference Lowe2004; Kong et al., Reference Kong, Akakin and Sarma2013; Xu et al., Reference Xu, Wu, Gao, Charlton and Bennett2020), are generally used to find regions within images that differ in certain properties, for example, in brightness, relative to neighboring regions. There are different forms of BD methods, and herein we will refer to two specific methods that outperform in both simulations and experiments: the “Laplacian of Gaussian” (LoG) method (Lindeberg, Reference Lindeberg1993, Reference Lindeberg1998) and the “Difference of Gaussian” (DoG) method (Lindeberg, Reference Lindeberg2015). A general sketch illustrating the workflow of BD algorithm is shown in Supplementary Figure S1.

Fig. 1. General scheme on how the coordinates and the intensity $\hat{I}$ of atomic columns are calculated starting from the blob detection algorithm output. The five steps followed on the procedure are further explain with details in the methodological section “Tailored adjustments for our specific experimental dataset”. Images on Steps 2, 3, and 5 show a simulated atomic column in TEM mode.

Before proceeding to the algorithm, we first introduce some notation and highlight relevant background information. The size of the images we consider is denoted as m × n, and the pixel lattice coordinates are given by the set Ω = {(i, j): i = 1, …, m, j = 1, …, n}. We introduce the mapping f(x, y) such that the pixel values of a TEM image are denoted as f _i,j =  f(x = i, y = j) > 0 for pixel indices (i, j) ∈ Ω. Note that in TEM images, areas with no sample are considered as measurements of the vacuum background, for example, the top-right corner of Figure 2 is unstructured noise; we denote f ^B to be the average pixel values in the vacuum background. In addition, we denote $f^\alpha$ to be α = 5th-percentile from all the pixel values{f _i,j:(i, j) ∈ Ω}. We define a two-dimensional, mean zero, spherical Gaussian density function as:

(1)$$G( {x, \;y, \;\sigma^2} ) = \Phi ( {( {x, \;y} ) , \;( {\mu_X, \;\mu_Y} ) = ( {0, \;0} ) , \;{\rm \Sigma } = \sigma^2I_2} ) = {\rm \;}\displaystyle{1 \over {2\pi {\rm \sigma }^2}}\exp \left({-\displaystyle{{x^2 + y^2} \over {2{\rm \sigma }^2}}} \right)$$

Fig. 2. Fluxional behavior in time-resolved in situ TEM experimental data. Images have been contrast-reversed following equation (2). (a) 2.5 s summed image of a CeO₂ nanoparticle, presenting a good signal-to-noise ratio which enables to efficiently resolve the position and the intensity of each present atomic column. (b) and (c) illustrate a 37.5 ms summed image at two different periods of time, where the red arrows indicate the disappearance of two atomic columns, this indicating that the summed image (a) is not representative for all the frames in the time series. (d) Time-resolved single frame of the experimental data, 2.5 ms, pointing out that the high content of noise (SNR = ~4) hinders the determination of the time the evolution of each atomic column. (e1) and (e2) represent a surface and bulk column, respectively, showing the difficulty of determining the central position and intensity of the atomic columns, especially on those at the surface where more dynamics are taking place.

for example, with mean (μ _X, μ _Y) = (0, 0) and the single-scale parameter σ ², where Φ denotes the standard bivariate Gaussian density function. We now introduce our full BD algorithm through the definitions and steps below, highlighting our additional tailoring to this specific problem.

Conventionally, the BD function of the Python's skimage package seeks bright pixels surrounded by darker ones (e.g., Supplementary Fig. S1a). However, the particular bright-field experimental data presented in this work was recorded initially with a bright background and black contrast of atomic columns, such that the regions of interest (the atomic columns) are darker than their surroundings. Hence, we introduce the following preprocessing step for experimental TEM images.

Definition 0. Define the function F(x, y) that generates a contrast reversal of the original image values f at (i, j) ∈ Ω as:

(2)$$F_{i, j} = F( {X = i, \;Y = j} ) = {\rm \;}\displaystyle{{\,f^B-f_{i, j}} \over {\,f^B-f^\alpha }}, \;$$

resulting in bright atomic columns (high pixel values) surrounded by dark background.

Remark on Definition 0. In equation (2), the numerator is reversing the contrast, while the denominator, $f^B-f^\alpha$, is used to normalize the observations. Additional subtraction by $f^\alpha$ in the denominator is performed such that the majority of the mapped values {F _i,j:(i, j) ∈ Ω} will be less than one, especially at the locations {(i, j)} of primary interest within the black contrast of an atomic column (where f _i,j ≤ f ^B usually holds). Such preprocessed values F _i,j will also tend to be distributed within [0, 1] quite densely, that is, 0 ≤ F _i,j ≤ 1 regardless of the raw observations f _i,j. As an added benefit, the BD parameters, as estimated below, under this normalized transformation may also be applied with similar detection power on new images with comparable SNRs, but independent of the new image's specific scaling. In practice, we have used the 5th-percentile for $f^\alpha$, as noted above; however, other percentiles in the range 1–10 may also be used depending on how variable or stable low percentile signal values behave.

Definition 1. The BD algorithm begins by convolving (discretely and finitely) the contrast reversed image {F _i,j:(i, j) ∈ Ω} with the Gaussian function G(x, y, σ ²) at the pixel locations Ω, illustrated in Supplementary Figure S1. Specifically, we define the convolution function C(x, y, σ ²) as:

(3)$$C( {x, \;y, \;\sigma^2} ) = G( {x, \;y, \;\sigma^2} ) {\rm \ast }F( {x, \;y} ) = \mathop \sum \limits_{( {i, j} ) \in {\rm \Omega }} G( {x-i, \;y-j, \;\sigma^2} ) \cdot F_{i, j}.$$

With this definition for the convolution function C(x, y, σ ²), we next introduce the two objective functions that are optimized within the LoG and DoG methods, respectively.

Remark on Definition 1. In addition to smoothing out the features arising from noise, the discrete Gaussian convolution is used to match our model Assumption 1 given later, and help determine the Gaussian parameters (both coordinates and variance) of the atomic columns (Supplementary Fig. S1).

Definition 2. We define L(x, y, σ ²) as the (scale-normalized) Laplacian, or second derivative of function C(x, y, σ ²), as:

(4)$$L( {x, \;y, \;\sigma^2} ) = {\rm \;}\sigma ^2\Delta C( {x, \;y, \;\sigma^2} ) = {\rm \;}\sigma ^2\mathop \sum \limits_{( {i, j} ) \in {\rm \Omega }} {\rm \Delta }G( {x-i, \;y-j, \;\sigma^2} ) \cdot F_{i, j}, \;$$

in which ΔC(x, y, σ ²) = C _xx(x, y, σ ²) + C _yy(x, y, σ ²) and ΔG(x, y, σ ²) = G _xx(x, y, σ ²) + G _yy(x, y, σ ²) denote the Laplacian operators of the function C(x, y, σ ²) in equation (3) and G(x, y, σ ²) in equation (1) with respect to x and y, respectively.

We further define Q(x, y, σ ²) as the (scale-normalized) first derivative of C(x, y, σ ²) with respect to σ ²:

(5)$$Q( {x, \;y, \;\sigma^2} ) = \sigma ^2C_{\sigma ^2}( {x, \;y, \;\sigma^2} ).$$

Remark on Definition 2. The scale parameter σ ² is included in both functions L(x, y, σ ²) and Q(x, y, σ ²) so that in addition to estimating the center locations of blobs, the minimization of the objective functions can also automatically match the exact scale/shape of the underlying blobs (see Lindeberg, Reference Lindeberg1998).

In practice, the BD algorithm evaluates the functions L(x, y, σ ²) or Q(x, y, σ ²) discretely at (x, y) = (i, j) ∈ Ω and $\sigma ^2 = \sigma _k^2 \in {\rm \Psi } = { {\sigma_1^2 , \;\cdots , \;\sigma_K^2 } }$, where $\sigma _1^2 \le \cdots \le \sigma _K^2$ are user-specified equally spaced BD variance parameters for discrete and approximate optimization. The Python function takes σ ₁, σ _K, and K as the input and automatically generates the corresponding arithmetical sequence Ψ. These sets generate a three-dimensional array of evaluation points for the final objective function. In the literature (Lindeberg, Reference Lindeberg2015), the partial derivatives appearing in equation (5) is usually approximated using discrete finite difference. Specifically, denoting Λ = Ω × {1, …, K} and $C[ {i, \;j, \;k} ] = C( {x = i, \;y = j, \;\sigma^2 = \sigma_k^2 } )$, for (i, j, k) ∈ Λ, we have

(6)$$C_{\sigma ^2}( {x = i, \;y = j, \;\sigma^2 = \sigma_k^2 } ) \approx \tilde{C}_{\sigma ^2}[ {i, \;j, \;k} ] \triangleq ( {C[ {i, \;j, \;k + 1} ] {-C} [{i, \;j, \;k} ] } ) /( {\sigma_{k + 1}^2 -\sigma_k^2 } ).$$

Hence, for (i, j, k) ∈ Λ, we have

(7)$$L( {x = i, \;y = j, \;\sigma^2 = \sigma_k^2 } ) = \tilde{L}[ {i, \;j, \;k} ] $$

which can be directly evaluated according to equation (4), and

(8)$$Q( {x = i, \;y = j, \;\sigma^2 = \sigma_k^2 } ) = \tilde{Q}[ {i, \;j, \;k} ] \triangleq \sigma _k^2 \tilde{C}_{\sigma ^2}[ {i, \;j, \;k} ].$$

The LoG method of the BD algorithm approximately minimizes the objective function L(x, y, σ ²) by searching for the local minimums across the 3D array $\tilde{L} = { {\tilde{L}[ {i, \;j, \;k} ] \colon ( {i, \;j, \;k} ) \in {\rm \Lambda }} }$, that is, those points that have values below a given threshold (δ) as compared with their immediate (3 × 3 × 3) neighbors in this array (see Supplementary Fig. S1e). The DoG branch does the same (local) minimization search, except for the array $\tilde{Q} = { {\tilde{Q}[ {i, \;j, \;k} ] \colon ( {i, \;j, \;k} ) \in {\rm \Lambda }} }$ instead.

Definition 3. The local minimum indices ${ {( {{\hat{x}}_r, \;{\hat{y}}_r, \;{\hat{k}}_r} ) \in {\rm \Lambda }\colon r = 1, \;\ldots , \;R} }$ of the 3D array $\tilde{L}$ given a threshold δ satisfy:

(9)$$ \eqalignb{& \tilde{L}[ {{\hat{x}}_r, \;{\hat{y}}_r, \;{\hat{k}}_r} ] < \tilde{L}[ {i, \;j, \;k} ] - \delta , \; \cr & {\rm \;for\;} \left\{{\matrix{ {i = {\hat{x}}_r-1, \;{\hat{x}}_r, \;{\hat{x}}_r + 1, \;} \cr {\,j = {\hat{y}}_r-1, \;{\hat{y}}_r, \;{\hat{y}}_r + 1, \;} \cr {k = {\hat{k}}_r-1, \;{\hat{k}}_r, \;{\hat{k}}_r + 1, \;} \cr } } \right.{\rm \;and\;}( {i, \;j, \;k} ) \ne ( {{\hat{x}}_r, \;{\hat{y}}_r, \;{\hat{k}}_r} ).$$

And an image's candidate blob set is hence defined as:

(10)$$B_\delta ^{\rm LoG} = { {b_r = ( {{\hat{x}}_r, \;{\hat{y}}_r, \;\sigma_{{\hat{k}}_r}^2 } ) \in {\rm \Omega } \times {\rm \Psi }\colon r = 1, \;\ldots , \;R} }.$$

Remark on Definition 3. The blob set $B_\delta ^{\rm DoG}$ from the other array $\tilde{Q}$ can be analogously defined given a threshold δ. With no further notice, the blob-like triplet b = (x, y, σ ²) will always be regarded as the circle centered at (x, y) with radius 2σ, including its edge. Instead of using radius $\sqrt 2 \sigma _{{\hat{k}}_r}$ as in the generic case described by skimage package (Lindeberg, Reference Lindeberg1993), we use the expanded radius $2\sigma _{{\hat{k}}_r}$ in order to have a larger blob area to compute the center of mass (as detailed below). Rather than outputting all local minimums, the function in the Python package preprocesses and removes largely overlapping blobs. To be specific, if the overlapping area of two of these blobs exceeds a given percentage (γ), then the smaller circle is eliminated from the set of blobs. Additionally, the generic version of BD only allows concentrically symmetric blobs, otherwise irregularly shaped atomic columns may be associated with multiple connected blobs. Ideally speaking, if the input image {f _i,j:(i, j) ∈ Ω} is well-conditioned, that is, with high SNR and well-separated nanoparticle structures, the number of detected blobs R should be equal (or close) to the number of atomic columns.

To summarize, for the optimal performance of a BD algorithm, the parameters that require careful initializations include the prefixed BD variance parameters Ψ (minimum and maximum sigma, and the number of steps between them), the threshold for local minimums δ, and the overlapping blob percentage γ (value between 0 and 1). As an example, Table 1 shows the parameters we have used for our simulated and experimental datasets. We have implemented LoG and DoG to our dataset and it turned out that both methods return excellent initial blobs for our following steps. Notice that the blobs are determined to within one pixel accuracy, sub-pixel refinement is performed at a later stage [see discussion of equation (14) below].

Table 1. Initialization Parameters Used in Our Noise-Free, Noisy, and Experimental Datasets.

Atomic Information Extraction and Refinements

Based on the general BD approach detailed above, we further introduce our tailored proposal for atom identification in time-resolved TEM image series. In particular, we utilize the output from a BD algorithm as our initial information regarding the atomic columns, we compute the blob area-based center of mass to refine estimates of column positions, and fit a weighted linear regression model to estimate the blob-based column intensities using our modeling assumptions detailed below. These steps, that are described in detail in the following paragraphs, are shown in Figure 1.

To eliminate the influence between different atomic columns due to the convolution equation (3) across the whole images, and to avoid the challenge of labeling detected blobs from the BD function with atomic columns’ indices, we choose not to implement BD globally on the whole image field of view. Instead, we implement separately within user-specified sub-region that roughly correspond to the approximate position of specific atomic columns (one sub-region per atomic column). Indeed, such idea also saves the computational power since the sizes of the arrays $\tilde{L}$ and $\tilde{Q}$ decrease significantly. Typically, a time series of image frames is considered for analysis. First, by summing all or a subset of the frames pixel-wise across the time series, the size and location of each search region can be chosen to guarantee that the sub-image is large enough to contain an atomic column, while sufficiently small to ensure that neighboring columns are not included. Then, the following tailored steps are applied sequentially to all atomic columns to recover the underlying nanoparticle structure within an image. Hence from now we will update the notation Ω as the indices set of the sub-images.

To begin, we first propose the following assumption, similar to Levin et al. (Reference Levin, Lawrence and Crozier2020), to model the pixel behavior of a specific atomic column.

Assumption 1. The pixel values within the atomic column can be exactly parametrized by the blob-like triplet $b_0 = ( {x_0, \;{\rm \;}y_0, \;{\rm \;}\sigma_0^2 } )$. In addition, every pixel (i, j) ∈ Ω within or near (as defined below) the circle b ₀ satisfies:

(11)$$f_{i, j} = \mu -A\cdot G( {x = x_0-i, \;y = y_0-j, \;\sigma^2 = \sigma_0^2 } ) , \;$$

where A is a scale for the Gaussian function, which is varying relative to the number of atoms present, and μ is the background intensity level (both defined below); and G(x, y, σ ²) is the Gaussian density function as defined in equation (1). Equivalently, the contrast behavior for the atomic column compared to its neighborhood can be derived directly from equation (11) as:

(12)$$A\cdot G( {x = x_0-i, \;y = y_0-j; \;\sigma^2 = \sigma_0^2 } ) = \mu -f_{i, j}.$$

Remark on Assumption 1. The parameters A and μ appearing in equations (11) and (12) are unknown and to be estimated as $\hat{A}$ and $\hat{\mu }$ in the steps below. Note that in TEM imaging, the background level μ is atomic-column specific since the intensity values in the surrounding region of the atomic column depend on the number of atoms present, and it differs from the value in the vacuum (Gonnissen et al., Reference Gonnissen, De Backer, den Dekker, Sijbers and Van Aert2017), that is, μ will differ from f ^B, in general.

Further note that in the following discussion, we derive our procedure with respect to {f _i,j:(i, j) ∈ Ω} not the contrast reversed {F _i,j:(i, j) ∈ Ω} to make later more straightforward comparisons with competing methods for column intensity estimation.

Recall in Definition 3 for the general introduction of BD approach, we define equation (10) as the set of candidate blobs from the LoG method (DoG analogously). For simplicity, in the following definition, Λ now represents the collection of sub-image pixels, and we continue using the same notation of equation (10) but for the BD output targeted at a specific atomic column region within an image.

Definition 4. Given the extracted sub-image's size, calibrated by the summed image sequence, the output of BD for a specific atomic column region is:

(13)$$B = { {b_1 = ( {{\hat{x}}_1, \;{\hat{y}}_1, \;\sigma_{{\hat{k}}_1}^2 } ) , \;\cdots , \;b_R = ( {{\hat{x}}_R, \;{\hat{y}}_R, \;\sigma_{{\hat{k}}_R}^2 } ) } } \subset {\rm \Lambda }.$$

Remark of Definition 4. In practice, we can implement both LoG and DoG methods, and combine the returned blobs together to obtain the set B as given in equation (13) after removing repeated and overlapping blobs. Unlike BD applied to the entire image, the number of blobs in the output set B for a single atomic column region is generally one (R = 1), however, in the current experimental dataset, there are some very noisy frames in which we find two (or even more) blobs within the image sub-region. When more than one blob is returned for a particular atomic column, reference coordinates are required to help pinpoint the optimal blob. With the summed image sequence mentioned above, which has a relatively high SNR, the average position of the atomic column can more easily be estimated, however, it then lacks all dynamics taking place across the sequence of images.

Definition 5. Define reference coordinates of an atomic column as b* = (x*, y*), that is, those from a summed sequence of images.

Whenever the set B has multiple blobs, the following proposition details how to proceed and denotes the final notation for the unique blob mapped to that atomic column.

Step 1. The blob b _r in candidate set B whose center $( {{\hat{x}}_r, \;{\hat{y}}_r} )$ has the shortest Euclidean distance to reference b* is selected to represent the specific atomic column when there is more than one returned. Denote the selected blob as $\hat{b} = ( {\hat{x}, \;\hat{y}, \;{\hat{\sigma }}^2} )$.

Remark on Step 1. In expectation, among our search grid Λ, the selected output $\hat{b} = ( {\hat{x}, \;\hat{y}, \;{\hat{\sigma }}^2} )$ will be the closest discrete location to the ground truth $b_0 = ( {x_0, \;{\rm \;}y_0, \;\;\sigma_0^2 } )$ defined in Assumption 1.

With the final blob estimate $\hat{b} = ( {\hat{x}, \;\hat{y}, \;{\hat{\sigma }}^2} )$ given, we now detail our tailoring for recovering the coefficients A and μ in model Assumption 1 (equation (12)), followed by our approach for more precisely estimating an (integrated) intensity for that atomic column.

We estimate the atom-specific background level μ to be the average of the pixel values near the edge of the blob circle $\hat{b} = ( {\hat{x}, \;\hat{y}, \;{\hat{\sigma }}^2} )$ which approximates $b_0 = ( {x_0, \;y_0, \;{\rm \;}\sigma_0^2 } )$.

Step 2. The atom-specific background estimator $\hat{\mu }$ for μ is computed as the average pixel values inside the annulus region with center $( {\hat{x}, \;\hat{y}} )$ and radii $1.8\hat{\sigma }$ and $2.4\hat{\sigma }$, respectively.

Given the lack of inherent precision for the discrete location coordinates $( {\hat{x}, \;\hat{y}} ) \in {\rm \Omega }$, we propose a center of mass analysis that converts integer-valued coordinates $( {\hat{x}, \;\hat{y}} ) \in {\rm \Omega }$ to the real-valued pair $( {\bar{x}, \;\bar{y}} )$, which provides sub-pixel level precision.

Step 3. Given the estimated background level $\hat{\mu }$, the center of mass is then computed within the circle defined by $\hat{b} = ( {\hat{x}, \;\hat{y}, \;{\hat{\sigma }}^2} )$ as

(14)$$\bar{x} = \displaystyle{{\sum\nolimits_{( {i, j} ) \in \hat{b}} {i\cdot {( {\hat{\mu }-f_{i, j}} ) }_ + } } \over {\sum\nolimits_{( {i, j} ) \in \hat{b}} {{( {\hat{\mu }-f_{i, j}} ) }_ + } }}, \;{\rm \;}\bar{y} = \displaystyle{{\sum\nolimits_{( {i, j} ) \in \hat{b}} {\,j\cdot {( {\hat{\mu }-f_{i, j}} ) }_ + } } \over {\sum\nolimits_{( {i, j} ) \in \hat{b}} {{( {\hat{\mu }-f_{i, j}} ) }_ + } }}, \;$$

where $( {\hat{\mu }-f_{i, j}} ) _ + { = } {\rm max}( {0, \;{\rm \;}\hat{\mu }-f_{i, j}} )$ is the non-negative weight for pixel $( {i, \;j} ) \in \hat{b}$.

Note that in equation (12), the scaling or amplitude A is now the only parameter that remains un-estimated if we further assume that $b_0\approx \bar{b} = ( {\bar{x}, \;\bar{y}, \;{\hat{\sigma }}^2} )$ and $G( {x = x_0-i, \;y = y_0-j; $ $\sigma^2 = \sigma_0^2 } ) \approx G( {x = \bar{x}-i, \;y = \bar{y}-j; \;\sigma^2 = {\hat{\sigma }}^2} )$ for all pixel coordinates (i, j) in or around circle $\bar{b}$. In fact, the same equation (12) holds for every pixel (i, j) within the blob $\bar{b}$ except the parameters $( {x_0, \;y_0, \;\sigma_0^2 } )$ replaced by $( {\bar{x}, \;\bar{y}, \;{\hat{\sigma }}^2} )$; hence, a weighted least square regression is proposed to estimate the amplitude A.

Step 4. The amplitude A can be estimated by the following regression problem:

(15)$$\hat{A} = \mathop {{\rm argmin}}\limits_{a > 0} \mathop \sum \limits_{( {i, j} ) \in \bar{b}} \bar{G}[ {i, \;j} ] \cdot ( {\hat{\mu }-f_{i, j}-a\cdot \bar{G}[ {i, \;j} ] } ) ^2, \;$$

where $\bar{G}[ {i, \;j} ] = G( {x = \bar{x}-i, \;y = \bar{y}-j; \;\sigma^2 = {\hat{\sigma }}^2} )$, and G is the Gaussian function defined in equation (1).

Our final step is to compute an integrated intensity estimate for the atomic column, and this can be done analytically via integration, based on the estimated model parameters derived above. For better comparisons with alternative methods for estimating the integrated intensity for an atomic column, we narrow the integration region to be $\tilde{b} = ( {\bar{x}, \;\bar{y}, \;( {\hat{\sigma }}^2/2) } )$, that is, the circle with center $( {\bar{x}, \;\bar{y}} )$ and radius $\sqrt 2 \hat{\sigma }$.

Step 5. With $\bar{b} = ( {\bar{x}, \;\bar{y}, \;{\hat{\sigma }}^2} )$ refined from BD estimate $\hat{b}$ and the estimated amplitude $\hat{A}$, we define the $( {\bar{x}, \;\bar{y}} )$-centered scaled Gaussian density function $\bar{g}( {u, \;v} )$ as $\bar{g}( {u, \;v} ) = \hat{A}\cdot G( {x =&EnLnBrk; \bar{x}-u, \;y = \bar{y}-v; \;\sigma^2 = {\hat{\sigma }}^2} )$. The corresponding estimated integrated intensity $\hat{I}$ is hence computed by integrating $\bar{g}( {u, \;v} )$ within circle $\tilde{b} = ( {\bar{x}, \;\bar{y}, \;( {\hat{\sigma }}^2/2) } )$, and it simplifies to:

(16)$$\hat{I} = \int_{( {u, v} ) \in \tilde{b}} {\bar{g}( {u, \;v} ) dudv = ( {1-{\rm \;}e^{{-}1}} ) \hat{A}, \;} $$

where e is the base of the natural logarithm. Note that this intensity is effectively a contrast reversed intensity, it initially increases as the number of atoms in the column increases. The same procedure must be applied to the simulated images employed to generate look-up tables if the user wants to convert the intensity $\hat{I}$ to the number of atoms in the column.

To summarize, based on the blob $\hat{b} = ( {\hat{x}, \;\hat{y}, \;{\hat{\sigma }}^2} )$ returned by the BD function, our proposed algorithm calculates the set of estimators ${ {\hat{\mu }, \;( {\bar{x}, \;\bar{y}} ) , \;\hat{A}, \;\hat{I}} }$, and finally outputs $( {\bar{x}, \;\bar{y}, \;\hat{I}} )$, that is, the refined coordinates and the integrated intensity, as the featured information of a specific atomic column.

Tailored Adjustments for Our Specific Experimental Dataset

For our experimental dataset, 25 atomic columns in total are expected (see Supplementary Fig. S2) within the region of interest, and every column has unique behaviors with respect to the variability of its center coordinates and intensity level, especially those columns at the surface of the nanoparticle where most of the fluxional behaviors take place. To address such atomic column variations, we propose some final adjustments to guarantee that we get reliable estimation using only a single pass through a sequence of images.

Step 6. The procedure for each image consists of two rounds of detection. The first considers a relatively large threshold δ and a wide column-wise sub-image selection for the initial input of BD, and the second narrows the sub-image region(s) while also decreasing the threshold level δ.

Remark on Step 6. In most image frames and atomic columns, it is clear by visual inspection that the output from the initial round is usually sufficient. However, in some extreme frames and column regions, due to low SNR and fluxionality within the system, the output may be more irregular and occasionally lead to poor performance. For these specific extraordinary cases, our second round helps optimize results since it is more concentrated at the average location and aims to detect the fainter blob that is very likely missing in the first round. By tuning the corresponding parameters, we are able to get excellent estimates uniformly over all image frames for all atomic columns.

Generation of the Simulated Dataset

To compare the performance of different atomic column tracking algorithms, an absolute reference is required. The absence of known ground truth in experimental data has led us to create a dataset based on simulations, where the location and occupancy of each atomic column can be extracted from a 3D model. In particular, our 3D model has been built using the freely available Rhodius software (Bernal et al., Reference Bernal, Botana, Calvino, López-Cartes, Pérez-Omil and Rodrı́guez-Izquierdo1998); the model consists of 25 Ce atomic columns (see Supplementary Fig. S2 for more information about the atomic model), and it has shape and size similar to the experimentally characterized CeO₂ nanoparticle shown in Figure 2 (and described in Section “Experimental data acquisition”). In addition, to determine the performance of algorithms under challenging conditions, the upper surface of the CeO₂ model exhibit columns with occupancies of 1, 2, or 3 atoms, situations that may occur at the surface of our experimental nanoparticle undergoing fast rearrangements.

The TEM simulations were performed using the multi-slice method, implemented in the Dr. Probe software package (Barthel, Reference Barthel2018). Each slice was of 0.19 Å in thickness and the electron-optical and imaging parameters were set to match those encountered for the experimental data of Figure 2. Specifically, we choose a third-order spherical aberration coefficient (C _s) of −30 μm, a fifth-order spherical aberration coefficient (C5) of 5 mm, and a defocus (C1) of –14 nm. All other aberrations (e.g., 2-fold and 3-fold astigmatism, coma, star aberration, etc.) were considered to be negligible. The image was simulated at 300 kV accelerating voltage with a beam convergence angle of 0.2 mrad and a focal spread of 4 nm, including partial temporal coherence and partial spatial coherence. An isotropic vibration envelope of 50 pm was applied during the image calculation, and the effect of the modulation transfer function (MTF) was turned off. The size of the simulated images was 512 × 512, leading to a pixel size of 0.097 Å/pixel (slightly smaller than the experimental one). Final simulated images are contrast-reversed, F _i,j, to show bright atomic columns.

To simulate the high degree of noise present in experimental time-resolved data, Poisson noise was added to the simulated images to match the SNR present in the experiment data. In the experimental images (see Section “Experimental data acquisition”), the SNR was estimated by taking the mean intensity in vacuum divided by the standard deviation in vacuum giving a value of ~4.00 for images with a time exposure of 1/400th s. To quantify the robustness of the algorithms to noise, ten different Poisson noise realizations created a series of ten different images. Supplementary Figure S3 shows the noise-free and noise-degraded simulated images.

The positions and intensities of the columns in the clean simulation (reference) and noise degraded simulations (degraded) was determined and compared. The real position/intensities (from the ground truth) and the outputs from these algorithms, that is, the error, was determined from the Euclidean distance (i.e., square root of the difference between the reference and degraded outputs).

Experimental Data Acquisition

To investigate the performance of the blob approach on real experimental data, we used time resolved data from a previously studied CeO₂ nanoparticle (Lawrence et al., Reference Lawrence, Levin, Miller and Crozier2019). An in situ TEM image series of a CeO₂ nanoparticle in a [110]-zone axis was recorded on an image corrected Thermo Fisher Titan 80-300 environmental TEM. The microscope was operated at 300 kV, running in ETEM mode with a pressure below 10⁻⁶ Torr. A Gatan K2 IS direct electron detector running in integration mode, allowed a high spatiotemporal resolution, yielding a fast acquisition of 400 frames per second (2.5 ms/frame) and pixel size of ~0.125 Å/pixel. In total, the CeO₂ nanoparticle was imaged for ~22 s, that is, ~9,000 frames. In order to drive strong fluxional behavior, a high electron fluence rate of 120,000 e⁻ Å⁻² s⁻¹ was employed (300 e⁻Å⁻² per frame) which resulted in electron beam reduction of CeO₂ and associated creation, migration, and annihilation of oxygen vacancies. Previously, it has been reported that this current is enough to induce structural rearrangements at the surfaces of CeO₂ nanoparticles (Möbus et al., Reference Möbus, Saghi, Sayle, Bhatta, Stringfellow and Sayle2011; Bhatta et al., Reference Bhatta, Ross, Sayle, Sayle, Parker, Reid, Seal, Kumar and Möbus2012; Bugnet et al., Reference Bugnet, Overbury, Wu and Epicier2017; Sinclair et al., Reference Sinclair, Lee, Shi and Chueh2017; Lawrence et al., Reference Lawrence, Levin, Boland, Chang and Crozier2021). The high concentrations of oxygen vacancies destabilize the nanoparticle surface leading to cation migration providing an ideal system for quantifying structural dynamics. For the current project, to simplify the analysis, imaging conditions and analysis procedures were developed to track cation columns only since they show up with much higher SNR in the noisy datasets. However due to the high temporal resolution, each individual image is significantly degraded by noise, with a SNR of ~4.