Optimized imaging using non-rigid registration

The scanning procedure creates an array of intensities that we call a frame. It can be interpreted as an image reflecting the properties of the investigated specimen on some domain. However, due to distortions the locations corresponding to the measured intensities are not known exactly. Thus, the frame is not just the image but is also related to an unknown positioning system determining what exactly the image represents. The question is how the information received from several frames $f_{j},\,j=1,\dots,n$, can be effectively combined and how a specific knowledge about the distortions could be obtained.

To outline our approach we introduce some notation. It is convenient to represent a frame as a function $f$. The array values of the frame define the function values $f(x)=f(x_{1},x_{2})$ for certain integer coordinates $x_{1}$ and $x_{2}$ on a regular, rectangular grid. Then, using an appropriate approximation method, e. g. bilinear interpolation, one can extend this function to intermediate real positions $x=(x_{1},x_{2})$. Consider two frames $f_{j}$ and $f_{k}$. We want to establish a correspondence between their positioning systems by approximating the deformation $\phi_{j,k}$ that gives a position $\phi_{j,k}(x)$ in frame $f_{j}$ for a position $x$ in frame $f_{k}$ in the sense that $f_{j}\circ\phi_{j,k}\approx f_{k}$. The process of finding $\phi_{j,k}$ is often referred as registration of $f_{k}$ and $f_{j}$.

A particularly simple version of a registration map is a rigid motion, in general a combination of a translation and a rotation. This is typically determined by first manually aligning the respective frames, using specific image features, sometimes followed by adjustments based on numerical indicators such as the cross-correlation of the adjusted frames. In the case of STEM images, collected pixel-by-pixel, rigid registration is inherently limited due to the nature of the acquisition process.

Our nonrigid series registration is founded on a variational approach that transforms two frames $f$ and $g$ into a common coordinate system with a nonparametric, nonrigid transformation $\phi$ such that $f\circ\phi\approx g$. The objective functional consists of the sum of two terms, the normalized cross-correlation of $f\circ\phi$ and $g$ as fidelity term, and the Dirichlet energy of the displacement, i. e. $\phi$’s deviation from the identity as regularizer. While this nonparametric transformation model is very flexible, an effective numerical realization is rather difficult due to the large number of local minima our objective functional exhibits. The proposed hybrid minimization strategy rests on several important mathematical concepts:

• 1.

  a multilevel scheme that processes the frames at different levels of resolution, going from coarse to fine;

• 2.

  minimization on each level based on a gradient flow,a generalization of the gradient descent concept;

• 3.

  an explicit time discretization of the gradient flow ODE combined with an automatic step size control to ensure convergence

Since feature extraction from noisy data is a challenging task by itself, we employ a registration approach that does not extract features but works directly on the image intensities. For a comprehensive introduction to image registration in general we refer the reader to the book by Modersitzki [16].

Ideally, assuming no noise, no distortions, and a perfect extension of the frames to real valued functions, one should expect that at any $x=(x_{1},x_{2})$ from the scanned domain the composition $(f_{j}\circ\phi_{j,k})(x)=f_{j}(\phi_{j,k}(x))$ has the same value as $f_{k}(x)$. In reality, one attempts to find a deformation $\phi_{j,k}$ that provides a good fit of $f_{j}\circ\phi_{j,k}$ to $f_{k}$ by replacing the ideal pointwise condition $(f_{j}\circ\phi_{j,k})(x)=f_{k}(x)$ with a suitable similarity measure. Such a measure quantifies the similarity of two frames $f_{j}\circ\phi_{j,k}$ and $f_{k}$. An overview of various popular similarity measures can also be found in Modersitzki [16].

To describe such a similarity measure, as usual, we denote by $\overline{f}=\frac{1}{|\Omega|}\int_{\Omega}f\operatorname{d\mathit{x}}$ the mean value of $f$ over the image domain $\Omega$. The standard deviation of $f$ over $\Omega$ is denoted by $\sigma_{f}=\sqrt{\frac{1}{|\Omega|}\int_{\Omega}(f-\overline{f})^{2}\operatorname{d\mathit{x}}}$. The normalized cross-correlation of two functions $f$ and $g$ over the domain $\Omega$ is defined as

This quantity is between $-1$ and $1$, while $\mathrm{NCC}[f,g]=1$ if and only if $f=ag+b$ for some real constants $a>0$ and $b$. Adopting the tradition to formulate a variational approach for determining the deformation as a minimization problem, we define the data term of our objective functional as

Without further constraints, finding a deformation $\phi=(\phi^{1},\phi^{2})$ among all vector-valued piecewise bilinear functions that minimizes $S[\phi]$ is a severely ill-posed problem. Therefore, we add the regularization term

to penalize irregular displacements and receive the objective functional

to be minimized. The larger the regularization parameter $\lambda>0$, the smoother the resulting deformation $\phi$. This regularization term was used in Modersitzki [16] and its use is often referred to as diffusion registration.

While the non-rigid transformation model (2) is very flexible, its numerical realization is rather difficult. Due to the high beam sensitivity of zeolites, subsequent frames deteriorate rapidly and the periodic nature of typical input frames renders the estimation of the deformation $\phi$ that yields the global minimum of $E[\phi]$ particularly challenging. In fact, even after including a regularization term, the objective functional $E$ typically still exhibits a large number of local minima. Therefore, the numerical minimization requires sophisticated algorithms designed to handle a large number of local extrema. Hence, our proposed minimization strategy rests on several conceptual pillars that guide the development of the algorithmic procedures. The remainder of the current section provides the details of each component of the solution process. Primarily, it is based on a multilevel scheme that processes different levels of resolution of the data to be registered (see Subsection 2.4.5). The registration is first done on the coarsest level, the result obtained on this level is prolongated to the next finer level and used as the initialization for the minimization on this level. The process is repeated until a registration at the resolution of the input images is achieved. The minimization on each level is based on the gradient flow concept which is described in Subsections 2.4.6 and 2.4.7).

A multilevel scheme serves as the outermost building block for the registration of a pair of frames. As indicated by numerous authors, e. g. [17, 16, 18] just to name a few, a coarse to fine registration strategy is a powerful tool to help prevent a registration algorithm from getting stuck in undesired local minima. The basic ideas of multilevel schemes go back to multi-grid algorithms [19, 20] and here we use the standard terminology, such as prolongation and restriction mappings, from that subject area. In particular, a hierarchy of computational grids is built, ranging from a very coarse grid level, e. g. $2^{4}\times 2^{4}=16\times 16$ pixels, to a fine level that matches the resolution of the input images, e. g. $2^{10}\times 2^{10}=1024\times 1024$. The exponent of $2$ used to build a grid in the hierarchy is called the grid level of the multiscale grids. The registration is first done on a coarse staring level, the result obtained on this level is prolongated to the next finer level and used as the initialization for the minimization on this level. The process is repeated until a registration at the resolution of the input images is achieved. The coarser the level, the fewer structures of the input images are preserved. Hence, the minimization at a certain level avoids all local minima induced by the small structures not visible at this level. Typically grids of the size $2^{d}\times 2^{d}$ or $(2^{d}+1)\times(2^{d}+1)$ are used for multilevel schemes since there are canonical prolongation and restriction mappings to transfer data from one level to the other for these kind of grids.

In the case of $2^{d}\times 2^{d}$ grids, the grid is interpreted in a cell-centered manner. When going from one level to a finer level, each grid cell is split into two times two cells (two in each coordinate direction). The prolongation of a function copies the value from the coarse grid cell to the corresponding four grid cells in the finer grid, while the restriction uses the average value of four grid cells in the fine grid as the value for the corresponding coarse grid cell. For $(2^{d}+1)\times(2^{d}+1)$ grids a mesh-centered approach is used. Refining such a grid is done by adding a new grid node between every pair of neighboring nodes in the coarse grid. The prolongation copies the values of the coarse grid nodes to the nodes at the same positions in the fine grid and determines the values for the fine grid nodes that are not in the coarse grid by bilinear interpolation of the neighboring coarse grid values. The restriction is the transposition of the prolongation operation, but rescaled row-wise such that if we have a value of one on every fine grid node, we get a value of one at every coarse grid node.

The minimization on a given level is based on a so-called gradient flow [21], a generalization of the gradient descent concept. As with the latter, the basic principle of the minimization is to go in the direction of steepest descent and formalized with the ordinary differential equation (ODE)

Here, $\operatorname{grad}_{G}$ denotes the gradient with respect to a scalar product $G$. In this notation, the well-known classical gradient in finite dimensions is the gradient with respect to the Euclidean scalar product. Since a gradient flow, as a gradient descent, is attracted by the “nearest” local minimum, and the scalar product determines how distances are measured, a suitable choice of the scalar product avoids undesired local minima. As shown in [17], the ODE (3) is equivalent to

where $E^{\prime}$ denotes the first variation of $E$ and $A$ is the bijective representation of $G$, such that $G(\zeta_{1},\zeta_{2})=\left<A\zeta_{1},\zeta_{2}\right>$ for all $\zeta_{1},\zeta_{2}$. This formulation illustrates the influence of the scalar product and leads us to choose it such that $A^{-1}$ is smoothing its argument. This way the descent will avoid non-smooth minimizers. A Sobolev $H^{1}$ inner product scaled by $\sigma$, i. e.

has proved to be a suitable choice. Here, we denote $X:Y=\sum_{ij}X_{ij}Y_{ij}$ for matrices $X,Y\in\mathbb{R}^{2\times 2}$, i. e. the Euclidean scalar product of the matrices interpreted as long vectors. Note that for the corresponding metric resulting from the inner product (5), the operator $A^{-1}$ is equivalent to one implicit time step of the heat equation with a step size of $\textstyle{\frac{\sigma^{2}}{2}}$, an operation widely known for its smoothing properties.

From (4) we see that the first variation of the objective functional $E$ (i. e. $E_{f,g}$) is necessary in order to use the gradient flow. With a mostly straightforward but somewhat lengthy calculation, one obtains

For the spatial discretization of a numerical solution to the weak form (6) we employ piecewise bilinear Finite Elements (FE) where the computational domain $\Omega$ is a uniform rectangular mesh. Introducing the FE method is beyond the scope of this paper, for a detailed introduction to the FE concept we refer to the textbook by Braess [22]. Of course, it is also possible to use different discretization approaches, for instance Finite Differences, but Finite Elements offer a canonical way to conveniently handle $\left<E^{\prime}[\phi],\zeta\right>$, the weak formulation of the first variation of $E$. To compute the integrals in the energy, its variations and the metric, we use a numerical Gauss quadrature scheme of order 3 (see, for example, §3.6 of Stoer and Bulirsch [23]).

The gradient flow ODE (4) is solved numerically with the Euler method, an explicit time discretization approach. Thus, $\partial_{t}\phi$ is approximated by $\textstyle{\frac{\phi^{l+1}-\phi^{l}}{\tau}}$, where $\phi^{l}$ denotes the deformation at the $l$-th iteration of the gradient flow and $\tau$ is the still-to-be-chosen step size. This leads to the update formula

Since the minimization does not require the knowledge of the trajectory of $\phi$ given by the ODE but only the equilibrium point, it is sufficient to use this kind of simple explicit method to solve the ODE. It is very important though to carefully choose the step size $\tau$ since the iteration can diverge with a poorly chosen step size. To this end we employ the so-called Armijo rule with widening [24]. This rule is a line search method used on the function

and selects the step size such that the ratio of the secant slope $\textstyle{\frac{\Phi(\tau)-\Phi(0)}{\tau}}$ and the tangent slope $\Phi^{\prime}(\tau)$ always exceeds a fixed, positive threshold $\rho\in(0,1)$. This means that the actual achieved energy decay (secant slope) is at least $\rho$ times the expected energy decay (tangent slope). Note that the specific choice of the step size control is not important – there are several popular step size controls that could be used here. It is crucial though to use a step size control that guarantees the convergence of the time discrete gradient flow (7).

The resulting multilevel descent algorithm is summarized in Algorithm 1.

Our multilevel descent algorithm is well-suited for the beam sensitive materials we consider here. Beam damage starts as a loss of high frequency information in the image as the atomic structure of the sample is lost with increasing electron irradiation. At the coarser stages of the registration, we find the distortion map for the large scale differences between the image frames. As the registration proceeds, the distortions for the earlier distortion map are used as the initial values for the next level. In the case where either the noise level in the input frames is extremely high or the material has changed from frame to frame, the distortions from the coarser level are retained and the regularization of the algorithm smoothly adjusts between the coarser grid points. The step size control described above means that our registration proceeds only to the point where the comparison between the frames is still meaningful. When either the noise level of the input or a significant difference in the two frames occurs, the registration algorithm is terminated, minimizing the artifacts of our method.

Since consecutive images in our input series were acquired directly one after another, we can assume the differences between two consecutive images, both regarding displacement and beam damage, to be small enough to use the identity as initialization for the minimization on the coarsest level in Algorithm 1. In case this assumption is not fulfilled it may be necessary to add a fiducial mark to the specimen to facilitate the registration algorithm.

After having a suitable deformation for each consecutive image pair we can chain those to get estimates for non-consecutive pairs. Let $f_{1},\ldots f_{n}$ denote the images in our series ordered by their time of acquisition. As outlined above we can estimate deformations $\phi_{2,1}$ and $\phi_{3,2}$ that match the first and second image pair respectively, i. e. $f_{2}\circ\phi_{2,1}\approx f_{1}$ and $f_{3}\circ\phi_{3,2}\approx f_{2}$. Then $\phi_{3,2}\circ\phi_{2,1}$ is a good initial guess for the deformation matching $f_{1}$ and $f_{3}$ because $f_{3}\circ(\phi_{3,2}\circ\phi_{2,1})=(f_{3}\circ\phi_{3,2})\circ\phi_{2,1}\approx f_{2}\circ\phi_{2,1}\approx f_{1}$. Therefore, our proposed registration algorithm with $\phi_{3,2}\circ\phi_{2,1}$ as initial guess allows us to reliably register $f_{3}$ to $f_{1}$. Iterating this provides a way to accurately register any frame of the series to the first frame (cf. Figure 3).

After registering all frames to the first one, we can fuse the information contained in the series by a suitable averaging operator ${\cal A}$ acting on the deformed frames, $g_{i}=f_{i}\circ\phi_{i,1}$, where $\phi_{1,1}$ is the identity mapping. Given a sequence of frames $\{g_{i}\}_{i=1}^{n}$ the first example is the mean

which provides the best least squares fit to the sequence. The arithmetic mean is a natural choice to calculate an average and is very easy to compute, although, depending on the input data, it may not be the best choice. In particular, if the input data contains noticeable outliers, which due to the scattering nature of the electrons in the HAADF STEM measurements is expected to be the case, the median is preferable to the arithmetic mean. Recall the definition of the median

which is the best $\ell_{1}$ fit to the sequence and therefore less sensitive to outliers. This selection of ${\cal A}$ is used in all subsequent experiments.

The estimated reconstructed image $f={\cal A}\{f_{i}\circ\phi_{i,1}\}$ from the previous section can be considerably improved by iteration of this basic procedure. Choosing $f_{1}$ as reference frame was a canonical but arbitrary choice, whereas the newly obtained average $f$ is a much more suitable reference frame. Thus, it is reasonable to repeat the process using $f$ as reference frame. Since $f_{1}$ was the reference frame when calculating $f$ the identity mapping is a good initialization for the registration of $f_{1}$ to $f$ and we can simply prepend $f$ to the input series, calling it $f_{0}$, and repeat the process. Here however, we don’t include $f_{0}$ when calculating the average since the original images $f_{1},\ldots f_{n}$ obviously contain all of our measurement information. Thus, an average of the deformed frames may be calculated via $f={\cal A}[\{f_{i}\circ\phi_{i,0}\}_{i}]$. This results in a further improved reconstructed image $f$ and thus the process can be repeated again, leading to an iteration procedure to calculate $f$.

The $k$-th estimate of $f$ is denoted by $f^{k}$. Initially, $f^{0}=f_{1}$ since the specimen is least damaged during the acquisition of $f_{1}$. We use $f_{0}=f^{k-1}$ as our reference frame at this stage and calculate $f^{k}$ using ${\cal A}[\{f_{i}\circ\phi_{i,0}^{k}\}_{i}]$ where $\phi_{i,0}^{k}$ denotes the deformation from $f_{i}$ to $f_{0}=f^{k-1}$ at the end of each iteration.

A critical issue in finding the deformations $\phi_{{\color[rgb]{0,0,1}i},0}^{k}$ is how to determine the initial guess for the nonrigid registration process we have described. Since the specimen can move considerably during the acquisition of the image series, we use the identity as an initial guess only for the transforms $\phi_{j+1,j}$ between two consecutive frames (note that $\phi_{1,0}^{0}$ is the identity). Then the initial guess for $\phi_{1,0}^{k+1}$ is the transform $\phi_{1,0}^{k}$, while for $\phi_{j+1,0}^{k}$ it is the composition $\phi_{j+1,j}^{k}\circ\phi_{j,0}^{k}$. Note that the estimation of $\phi_{j+1,0}^{k}$ depends on $\phi_{j,0}^{k}$, so as indicated before, these deformations need to be determined one after another. As an example, in Figure 4 we present frame 9 registered to frame 1.

This core methodology opens an avenue to further algorithmic options which facilitate an even increased denoising effectiveness at the expense of additional computational effort (see Subsection 3.5.2). The point is that within some regime, in principle, additional computational work gains increased information quality. The full strategy is summarized in Algorithm 2.

Note that after the first time the average has been computed, i. e. when $k\geq 2$, the reference frame $f_{0}$ is considerably less noisy than the input data. Thus, it is then possible to reduce the regularization parameter $\lambda$ when calculating $\phi^{k}_{i,0}$ for $k\geq 2$. Here, we typically use $0.1\,\lambda$ as regularization parameter for $k\geq 2$, i. e. the same regularization parameter value is used for all $k\geq 2$, but this value is smaller than the one used for $k=1$.

The inevitable movement of the sample during the scanning process combined with how STEM rasters the electron probe on the sample to acquire an image leads to a spatial distortion in the STEM images. In particular, Figure 4 shows a significant translation of the field of view encountered during the acquisition of the series of frames. In order to retain as much information as possible and ease the registration process by minimizing frame-to-frame distortions we introduce next a simple strategy for capturing the most prominent “macroscale” translational parts of the distortion. To formulate this model we need to make assumptions on the movement of the sample and need to denote some quantities inherent to the STEM scanning process.

Let $t$ denote the time it takes to scan a single line in the image and let $t_{f}$ denote the flyback time, i. e. the amount of time it takes to move the STEM probe from the end of one line to the beginning of the next line. As first order approximation, we assume that the sample undergoes a constant translational movement $v=(v_{1},v_{2})\in\mathbb{R}^{2}$ while a STEM image is acquired. The movement vector may differ from frame to frame though. Furthermore, we assume that the STEM probe starts to scan the sample at position $(0,0)$. Thus, when we scan the first pixel, corresponding to the position $(0,0)$, we also scan position $(0,0)$ of the sample. When scanning the last pixel in the first line, corresponding to the position $(1,0)$, the time $t$, i. e. the time to scan a line, has passed and the sample has moved by $tv$. Hence, we are actually scanning the position $(1,0)-tv$ of the sample. Let $h$ denote the height of a STEM line as fraction of $1$, i. e. $h=\frac{1}{N-1}$ where $N$ denotes the number of lines in the STEM image. Then, when scanning the first pixel of the second line, corresponding to the position $(0,h)$, the time $t+t_{f}$ has passed and we are actually scanning the position $(0,h)-(t+t_{f})v$. Consequently, when scanning the last pixel of the second line, corresponding to the position $(1,h)$, the time $2t+t_{f}$ has passed and we are actually scanning the position $(1,h)-(2t+t_{f})v$.

In general, when scanning a position $x_{1}\in[0,1]$ in line $k$, corresponding to $(x_{1},kh)$, a time of $tx_{1}+(t+t_{f})k$ has passed and we are actually scanning the position $(x_{1},kh)-(tx_{1}+(t+t_{f})k)v$. Therefore, assuming that the distortion is linearly interpolated between the lines, the image domain is deformed by the linear mapping

In other words, denoting the matrix by $M$, instead of $f(x)$ the intensity acquired at position $x$ is actually showing $f(Mx)$. Thus, using the inverse of $M$, it would be possible to remove the distortion caused by a constant translation of the sample during the acquisition of the frame.

Algorithm 2 allows us to generate a reconstruction of the series that structurally closely resembles the first frame $f_{1}$ due to the fact that $f_{1}$ was used as initial reference frame in the algorithm. In particular, if we drop the outer loop, i. e. for $K=1$, the reconstruction actually is a denoised version of $f_{1}$ (note that this is not guaranteed for $K\geq 2$). Using a slight alteration of Algorithm 2, we can register all frames to the $j$-th frame. This is achieved by applying the algorithm with $K=1$ and, without the averaging step, once on the series $f_{j-1},\ldots,f_{1}$ and once on the series $f_{j+1},\ldots,f_{n}$, each time using $f_{j}$ as reference frame. Having the transformations of all frames to the $j$-th frame, we can create a denoised reconstruction $\tilde{f}_{j}$ of the $j$-th frame and thus of every input frame. Furthermore, we can use the algorithm to register all of the denoised frames to the first of these frames $\tilde{f}_{1}$. Contrary to the initial registration of the input frames to $f_{1}$, the registration of the denoised frames can be done with a considerably lower regularization parameter $\lambda$ and at a higher precision since this registration is not hampered by the large amount of noise found in the input frames. Due to their improved accuracy, the deformations determined based on $\tilde{f}_{1}\ldots,\tilde{f}_{n}$ can be used to generate an improved reconstruction from our original input $f_{1},\ldots,f_{n}$ by using these deformations when averaging the input frames with $f_{1}$ as reference frame. If desired, this process can even be iterated, i. e. instead of only generating an improved reconstruction of $f_{1}$ one can generate an improved reconstruction of every frame and again create an improved average of the original data by registering the improved reconstructions and so on, which can be seen as an extension of the $k$-loop idea in Algorithm 2. Let us emphasize that all reconstructions here are always obtained by averaging the original noisy input frames, only the deformations used for the averaging were determined based on reconstructions we calculated previously.

The main drawback of this variant is the vastly increased computational cost compared to our original algorithm. Since the original algorithm essentially needs to be run $n+1$ times, once for each input frame to get a denoised version of this frame and once to create the improved reconstruction based on the registration of the denoised frames, the computational cost is increased by a factor of $n+1$. Thus, it is only feasible to use this approach if the computational cost is not a concern or if the number $n$ of input frames is moderate.