Background

Electron tomography has been a powerful tool in determining 3-D structures and characterization of nanoparticles in the biological, medical, and materials sciences [1-3]. The method is carried out by acquiring a series of 2-D projection images of an object and then using these 2-D projections to reconstruct the 3-D object. Using the transmission electron microscope, these projections are collected at a number of different orientations, typically by tilting the sample about a fixed tilt axis [4], while other dual axis tilting schemes also exist [5]. A demonstration of the projection scheme is shown for a 2-D object in Figure 1. We will focus only on the case of a single fixed tilt axis in this paper, although our methods can easily be translated to dual axis schemes.

Figure 1
figure 1

1-D projections are taken of a 2-D object. The small ball along the edge is not projected in the 0° projection straight down due to the limited projection range. However, at the higher angles, this mass is now projected, which will affect an alignment based on the center of mass of these projections.

Full size image

Ideally, between two consecutive projections acquired at nearby tilts of the sample, one would observe only a small rotation of the projected image. However, due to unavoidable mechanical limitations, significant translation shifts are present. Therefore, the projections must be aligned into a common coordinate system to be properly interpreted. Once the projections are aligned, they can then be merged to approximate the 3-D structure of the sample. The alignment is a crucial part of the process, for the resolution of the reconstructed 3-D structures are limited to the accuracy in the alignment. In this paper, we demonstrate a new mathematically justified method for the alignment based on the apparent motion of the center of mass of many 2-D cross-sections of the sample.

Over the years, many traditional alignment techniques have been developed by the biological sciences [6]. The most commonly practiced are correlation techniques, feature tracking, and fiducial marker tracking. Correlation techniques are performed by selecting one of the projections as a reference image and aligning each pair neighboring images by selecting the cross-correlation peak between the images for the shift [7]. This method has been proven useful but can yield poor results, as small cumulative errors may result in a serious drift of the sample [8]. As we will show, cross-correlation will not recover the correct alignment even for noise-free data subjected to random shifts. The current work finds a solution without this deficiency.

Fiducial marker tracking is done by decorating the sample with small high-density particles that create high contrast in the projection images [9-12]. Individual markers are then identified in all projections. The alignment is determined based on tracking of the path of each marker through the projections. This method can be very accurate but requires a lot of manual interaction to properly locate and center the markers. The main drawback of marker tracking is that the markers will be present in the reconstruction and must be removed for accurate characterization of the sample. Since the markers are of such high density, the reconstruction of the markers will inevitably mix with the reconstruction, making the task of removal nontrivial and possibly inaccurate.

Feature tracking uses regions of high contrast or intensity as fiducial markers [13,14]. It requires the identification of suitable regions of high contrast that remain visible throughout the tilt series.

Others have begun to perform alignment techniques based on a refinement approach [6]. After a coarse alignment from cross-correlation, one proceeds in computing an initial 3-D reconstruction. This 3-D reconstruction is then reprojected and compared with the original projections. A new alignment arises from aligning the reprojected reconstruction with the original projections, and this process is iterated until convergence is met. In our experience with this method, the reconstruction always satisfies the projections, even if they’re misaligned, so that insignificant refinement occurs from updating.

Most recently, Scott et al. [15] introduced a technique based on the observation that as the sample is tilted about a fixed axis, the center of mass of the sample will spin in a circle, and if the center of mass is on the tilt axis, then it remains fixed. In this way, it was suggested to shift each projection so that the center of mass in each projection is fixed on a point and taking the line through this point parallel to the axis of rotation as the tilt axis. We believe this is not always applicable and can yield poor results in many settings. First, it requires a tilt series in which the total projected volume is fixed for each projection. However, in most practical settings, some mass will move in and out of the projection range as the sample is tilted, which will then significantly affect the location of the center of mass within the projection along both axes of the projection images. This transition of mass must be accounted for, as this transition will be along the edges of the projections, far from the center, and will thus weigh heavily on the calculated center of mass. Figure 1 demonstrates this transition of mass, with the small ball located on the left edge of the object that has only been projected at certain angles. An additional drawback is that using only the single center of mass point in each projection removes the use of any local structure of the projections as criteria for alignment.

In this paper, we give an alignment method that makes more detailed use of the path of the projected center of mass along many cross-sections of the object, perpendicular to the axis of rotation. In an ideal experiment, points on the sample move in circular trajectories. We define a viable path as the projection of such a circular orbit. By simple calculation, we derive an equation which describes all such viable paths of the projected centers of mass, as opposed to the one trivial path of a single point. From here, we show how one can determine a shift for each projection so that the center of mass of all cross-sections perpendicular to the axis of rotation nearly follows a viable path. In this way, since all cross-sections are considered in our alignment method, we will be able to avoid problems involved with error in the calculated centers of mass due to transition of volume in and out of the projections, and we maintain local analysis of the projections as means for the alignment. Additionally, our model aligns the projections based on the rotation about a chosen axis, so that manual interaction for determining the positioning of the tilt axis is avoided. In general, our method can be considered more statistically accurate, and we will show that it provides very dependable alignment and definitively improves the resolution of the reconstruction.

Methods

Notation

The 3-D density function for reconstruction will be denoted f(x,y,z)=f(x,(y z)), with (y z) a 2-D row vector. The data generated are the projections of f in the z-axis, about rotations around the x-axis. A rotation of f through θ about the x-axis can be written as:

$$ f\left(x,\left(y\kern1em z\right){Q}_{\theta}\right),\kern1em \mathrm{where}\kern1em {Q}_{\theta }=\left(\begin{array}{cc} \cos \theta & - \sin \theta \\ {} \sin \theta & \cos \theta \end{array}\right). $$

A projection about the rotation θ is then defined as:

$$ {P}_{\theta }(f)\left(x,y\right)=\underset{\mathbb{R}}{\int }f\left(x,\left(y\kern1em z\right){Q}_{\theta}\right)\kern1em dz. $$

We note that for each fixed x=x 0, P θ (f)(x 0,y) only contains information from f(x 0,y,z), and therefore, many of the alignment and reconstruction processes can be considered as 2-D rather than 3-D. Therefore, for convenience, we will sometimes denote:

$$ {f}_x\left(y,z\right)=f\left(x,y,z\right)\kern1em \mathrm{and}\kern1em {P}_{\theta}\left({f}_x\right)(y)={P}_{\theta }(f)\left(x,y\right). $$

In practice, we are given the unaligned data; therefore, we will regularly refer to the misaligned projections, denoted by \( {\overset{\sim }{P}}_{\theta }(f) \). We define these projections as:

$$ {\overset{\sim }{P}}_{\theta }(f)\left(x,y\right)={P}_{\theta }(f)\left(x-{x}_{\theta },y-{y}_{\theta}\right), $$

where the coordinates (x θ ,y θ ) are the shifts to be determined for the alignment. Similarly, we will denote:

$$ {\overset{\sim }{P}}_{\theta}\left({f}_x\right)(y)={P}_{\theta}\left({f}_x\right)\left(y-{y}_{\theta}\right), $$

where in this instance the shift x θ is not included. We do not include it, for determining the shifts x θ is a much more trivial task, so that most of our work here focuses on determining y θ after the x-axis alignment is completed.

We will denote the total mass about a cross-section x by \( {M}_x=\underset{{\mathbb{R}}^2}{\int }{f}_x\left(y,z\right)\kern1em dy\kern1em dz \). Then, the coordinates for the center of mass of a cross-section are denoted as:

$$ {c}_x^y=\frac{1}{M_x}\underset{{\mathbb{R}}^2}{\int }{f}_x\left(y,z\right)y\kern1em dy\kern1em dz,\kern1em {c}_x^z=\frac{1}{M_x}\underset{{\mathbb{R}}^2}{\int }{f}_x\left(y,z\right)z\kern1em dy\kern1em dz $$

We will denote the center of mass of a projected cross-section of f by:

$$ {t}_x^{\theta}\kern0.3em =\kern0.3em \frac{1}{M_x}\underset{\mathbb{R}}{\int }{P}_{\theta}\left({f}_x\right)(y)y\kern1em dy,\kern1.30em \mathrm{and}\kern1.30em {\overset{\sim }{t}}_x^{\theta}\kern0.3em =\kern0.3em \frac{1}{M_x}\underset{\mathbb{R}}{\int }{\overset{\sim }{P}}_{\theta}\left({f}_x\right)(y)y\kern1em dy $$

We take the conventional L p norm (denoted by ∥·∥ p ) of a function, say g, defined over ℝ n to be:

$$ \parallel g{\parallel}_p^p=\underset{{\mathbb{R}}^n}{\int}\Big|g(x){\Big|}^p\kern1em dx. $$

Similarly, for a vector x∈ℝ n, we take the ℓ p norm (denoted ∥·∥ p ) to be:

$$ \parallel x{\parallel}_p^p=\sum_{i=1}^n\Big|{x}_i{\Big|}^p. $$

Theoretical model

In practice, we are given the set of misaligned angular projections:

$$ {\overset{\sim }{P}}_{\theta_i}(f)\left(x,y\right),\kern1em \mathrm{f}\mathrm{o}\mathrm{r}\kern1em i=1,2,\dots, k. $$

Typically, the number of projections, k, can be from 50 to 200, with maximum tilts of ± 70°. The domain is of course limited, but for theoretical purposes, we will assume that the domain for y is all of ℝ. The problem is then to approximate the set of shifts \( \left({x}_{\theta_i},{y}_{\theta_i}\right) \) for alignment, so that \( {\left\{{\overset{\sim }{P}}_{\theta_i}(f)\left(x,y\right)\right\}}_{i=1}^k \) correspond to the aligned projections \( {\left\{{P}_{\theta_i}(f)\left(x,y\right)\right\}}_{i=1}^k \). Determining the shifts for the x-axis is much simpler, since the x-axis is the axis of rotation. We simply observe that the total mass in each cross-section should remain fixed, so that:

$$ {M}_x=\underset{{\mathbb{R}}^2}{\int }{f}_x\left(y,z\right)\kern1em dy\kern1em dz=\underset{\mathbb{R}}{\int }{P}_{\theta_i}\left({f}_x\right)(y)\kern1em dy. $$
((1))

Based on this simple observation, one should be able to approximate all shifts \( {x}_{\theta_i} \) based on a ‘conservation of mass’ approach. We design a ‘global’ alignment method for determining these shifts, by taking \( {x}_{\theta_i} \) to be the shift which minimizes the difference between the observed mass in each cross-section of \( {\overset{\sim }{P}}_{\theta_i}(f)\left(x-{x}_{\theta_i},y\right) \) and the average mass of all projections in each cross-section. More precisely, we let:

$$ \begin{array}{c}{x}_{\theta_i}\kern0.3em = \arg \underset{x^{\ast }}{ \min}\kern0.3em {\parallel \underset{\mathbb{R}}{\int }{\overset{\sim }{P}}_{\theta_i}\left(\kern0.60em f\right)\left(x\kern0.3em -\kern0.3em {x}^{\ast },y\right)\kern1em dy\kern0.3em -\kern0.3em \frac{1}{k}\sum_{l=1}^k\kern0.3em \left(\underset{\mathbb{R}}{\int }{\overset{\sim }{P}}_{\theta_l}(f)\left(x,y\right) dy\right)\kern0.3em \parallel}_1 .\end{array} $$
((2))

Of course, the averaged term, \( \frac{1}{k}\sum_{l=1}^k\left(\underset{\mathbb{R}}{\int }{\overset{\sim }{P}}_{\theta_l}\right.\left.(f)\left(x,y\right)\kern1em dy\right) \), is subject to error since the projections are not yet aligned, so the determination of each \( {x}_{\theta_i} \) is iterated a few times until there is no change. The number of iterations will depend on just how large the offset of the projections are, but we have typically observed no change in each \( {x}_{\theta_i} \) after just two iterations. A demonstration of this x-axis alignment is given in Figure 2.

Figure 2
figure 2

Images demonstrating the alignment along the x -axis. (a) 2-D projection image taken at a 30°tilt about the x-axis. (b) 1-D projection of (a) onto the x-axis. (c) 1-D projections onto the x-axis of all 2-D projections taken at different tilts about the x-axis. The misalignment is clearly shown in (c), as the 1-D projections should all be nearly the same. (d) Same 1-D projections in (c), shown after alignment is performed along the x-axis.

Full size image

One could also perform a similar ‘local’ method, by comparing the consecutive projections to each other instead of the average. This approach is subject to cumulative error in the alignment similar to cross-correlation; therefore, we avoid this approach.

From here forth, we will now assume that the \( {x}_{\theta_i} \) have been accurately determined, and consider each cross-section. For alignment along the y-axis, we again want to make use of physical properties. It has been noted, as f x (y,z) is rotated about the origin, the center of mass \( \left({c}_x^y,{c}_x^z\right) \) will spin in a circle around the origin. It is not immediately clear, however, how this property can be observed within the projections and used for alignment. Computing the center of mass of a projected slice, we obtain:

$$ \begin{array}{lll}{t}_x^{\theta_i}& =\frac{1}{M_x}\underset{\mathbb{R}}{\int }{P}_{\theta_i}\left({f}_x\right)(y)y\kern1em dy\kern2em & \kern2em \\ {}& =\frac{1}{M_x}\underset{\mathbb{R}}{\int}\left(\underset{\mathbb{R}}{\int }{f}_x\left(\left(y\kern1em z\right){Q}_{\theta_i}\right)y\kern1em dz\right)\kern1em dy\kern2em & \kern2em \\ {}& =\frac{1}{M_x}\underset{\mathbb{R}}{\int}\underset{\mathbb{R}}{\int }{f}_x\left(\alpha, \beta \right)\left(\alpha \cos {\theta}_i-\beta \sin {\theta}_i\right)\kern1em d\alpha \kern1em d\beta \kern2em & \kern2em \\ {}& =\kern0.60em \frac{ \cos {\theta}_i}{M_x}\kern0.3em \underset{\mathbb{R}}{\int}\kern0.3em \underset{\mathbb{R}}{\int }{f}_x\left(\alpha, \kern0.3em \beta \right)\alpha \kern1em d\alpha \kern1em d\beta \kern0.3em -\kern0.3em \frac{ \sin {\theta}_i}{M_x}\kern0.60em \underset{\mathbb{R}}{\int}\underset{\mathbb{R}}{\int }{f}_x\left(\alpha, \kern0.3em \beta \right)\beta \kern0.3em d\alpha \kern0.3em d\beta \kern2em & \kern2em \\ {}& ={c}_x^y \cos {\theta}_i-{c}_x^z \sin {\theta}_i,\kern2em & \kern2em \end{array} $$

where we applied the substitution \( \left(\alpha \kern1em \beta \right):=\left(y\kern1em z\right){Q}_{\theta_i} \). This tells us that the center of mass of each projected cross-section should follow the path given by:

$$ {t}_x^{\theta_i}={c}_x^y \cos {\theta}_i-{c}_x^z \sin {\theta}_i,\kern1em \mathrm{f}\mathrm{o}\mathrm{r}\kern1em i=1,2,\dots, k. $$
((3))

This equation gives us a local relationship between the relative positioning of all of the projections to use for the alignment. As discussed earlier, in [15], it was simply noted that if the center of mass is located at the origin on the tilt axis, then it does not move under rotations about that axis. This observation can be made through similar computations where the integrand is first taken over x, and then, the center of mass is computed for the total sum of the cross-sections, that is:

$$ {t}^{\theta_i}=\frac{1}{M}\underset{{\mathbb{R}}^2}{\int }{P}_{\theta_i}(f)\left(x,y\right)\kern0.3em dx\kern0.3em y\kern0.3em dy={c}^y \cos {\theta}_i-{c}^z \sin {\theta}_i, $$
((4))

where c y and c z here denote the center-of-mass coordinates along the y- and z-axes, respectively, independent of x, and M denotes the total mass of f. Therefore, it is suggested to shift each projection so that \( {t}^{\theta_i}=0 \) for all i, so that c y=c z=0. While this approach is theoretically sound in an ideal setting, summing over x immediately removes any consideration of local behavior of the projections of f. As we will show, in many settings, this simplification can be a major drawback.

Therefore, our approach is to determine a sequence of shifts so that for each cross-section there exists some deterministic center of mass \( \left({c}_x^y,{c}_x^z\right) \) so that Equation 3 is nearly satisfied. With this in mind, let us denote:

$$ \varTheta =\kern0.3em \left(\begin{array}{cc} \cos {\theta}_1& \kern0.3em - \sin {\theta}_1\\ {} \cos {\theta}_2& \kern0.3em - \sin {\theta}_2\\ {}\vdots & \vdots \\ {} \cos {\theta}_k& \kern0.3em - \sin {\theta}_k\end{array}\right),\kern1em {c}_x\kern0.3em =\kern0.3em \left(\begin{array}{c}{c}_x^y\\ {}{c}_x^z\\ {}\end{array}\right),\kern1em \mathrm{and}\kern1em {t}_x\kern0.3em =\kern0.3em \left(\begin{array}{c}{\overset{\sim }{t}}_x^{\theta_1}\\ {}{\overset{\sim }{t}}_x^{\theta_2}\\ {}\vdots \\ {}{\overset{\sim }{t}}_x^{\theta_k}\end{array}\right). $$

We note that from the acquired projection data we can compute both Θ and t x . Now from Equation 3, if our alignment is good, then for each cross-section x, there should exist some c x so that Θ c x ≈t x . Therefore, in order to yield a good alignment, we would like to determine:

$$ {y}_{\varTheta }=\left(\begin{array}{c}{y}_{\theta_1}\\ {}{y}_{\theta_2}\\ {}\vdots \\ {}{y}_{\theta_k}\end{array}\right), $$

so that there exist some c x satisfying:

$$ \varTheta {c}_x\approx {t}_x+{y}_{\varTheta },\kern1em \mathrm{f}\mathrm{o}\mathrm{r}\ \mathrm{all}x, $$

or equivalently:

$$ \underset{c_x}{ \min}\parallel \varTheta {c}_x-\left({t}_x+{y}_{\varTheta}\right)\underset{2}{\overset{2}{\parallel }}\approx 0\kern1em \mathrm{f}\mathrm{o}\mathrm{r}\ \mathrm{all}x. $$

In practice, we will have some finite number of cross-sections, say x j , for j=1,2,…n. Then, we would like solve the minimization problem:

$$ \underset{y_{\varTheta }}{ \min}\left(\sum_{j=1}^n\underset{c_{x_j}}{ \min}\parallel \varTheta {c}_{x_j}-\left({t}_{x_j}+{y}_{\varTheta}\right){\parallel}_2^2\right) $$
((5))

Now we can compute the minimization over c x directly. Given Θ and t x , the least square solution \( {c}_x^{\ast } \), to \( \parallel \varTheta {c}_x-\left({t}_x+{y}_{\varTheta}\right){\parallel}_2^2\kern0.3em \) :

$$ {c}_x^{\ast }= \arg \underset{c_x}{ \min}\parallel \varTheta {c}_x-\left({t}_x+{y}_{\varTheta}\right)\underset{2}{\overset{2}{\parallel }}, $$

can simply be found by differentiation so that:

$$ \begin{array}{lll}& \left(\frac{\partial }{\partial {c}_x^y}\parallel \varTheta {c}_x-\left({t}_x+{y}_{\varTheta}\right)\parallel {\kern1.60em }_2^2\right)\left|{\kern1.60em }_{c_x={c}_x^{\ast }}=0\kern1em \mathrm{and}\right.\kern2em & \kern2em \\ {}& \left.\left(\frac{\partial }{\partial {c}_x^z}\parallel \varTheta {c}_x-\left({t}_x+{y}_{\varTheta}\right)\parallel {\kern1.60em }_2^2\right)\right|{\kern1.60em }_{c_x={c}_x^{\ast }}=0.\kern2em & \kern2em \end{array} $$

Solving these equations, the solution can be found to be:

$$ {\varTheta}^{+}\left({t}_x+{y}_{\varTheta}\right), $$

where Θ + denotes the pseudo-inverse of Θ, given by Θ +=(Θ T Θ)−1 Θ. It should be noted that Θ T Θ is a 2×2 matrix with entries:

$$ \begin{array}{lll}{\left({\varTheta}^T\varTheta \right)}_{11}& =\sum_{i=1}^k\overset{2}{ \cos }{\theta}_i,\kern1em {\left({\varTheta}^T\varTheta \right)}_{21}={\left({\varTheta}^T\varTheta \right)}_{12}\kern2em & \kern2em \\ {}& =-\sum_{i=1}^k \cos {\theta}_i \sin {\theta}_i,\kern1em {\left({\varTheta}^T\varTheta \right)}_{22}=\sum \overset{2}{ \sin }{\theta}_i,\kern2em & \kern2em \end{array} $$

which is clearly invertible and without any notable computational cost.

Then, our minimization in Equation 5 becomes:

$$ \begin{array}{lll}\underset{y_{\varTheta }}{ \min }& \left(\sum_{j=1}^n\parallel \varTheta {\varTheta}^{+}\left({t}_{x_j}+{y}_{\varTheta}\right)-\left({t}_{x_j}+{y}_{\varTheta}\right){\parallel}_2^2\right)\kern2em & \kern2em \\ {}& =\underset{y_{\varTheta }}{ \min}\left(\sum_{j=1}^n\parallel \left(\varTheta {\varTheta}^{+}-I\right)\left({t}_{x_j}+{y}_{\varTheta}\right)\underset{2}{\overset{2}{\parallel }}\right).\kern2em \end{array} $$
((6))

If we let:

$$ A=\left(\begin{array}{c}\varTheta {\varTheta}^{+}-I\\ {}\varTheta {\varTheta}^{+}-I\\ {}\vdots \\ {}\varTheta {\varTheta}^{+}-I\end{array}\right),\kern1em \mathrm{and}\kern1em b=\left(\begin{array}{c}\left(\varTheta {\varTheta}^{+}-I\right){t}_{x_1}\\ {}\left(\varTheta {\varTheta}^{+}-I\right){t}_{x_2}\\ {}\vdots \\ {}\left(\varTheta {\varTheta}^{+}-I\right){t}_{x_n}\end{array}\right), $$

then the minimization problem in Equation 6 is equivalent to solving a standard least squares problem:

$$ \underset{y_{\varTheta }}{ \min}\parallel A{y}_{\varTheta }-b{\parallel}_2^2. $$
((7))

Practical implementation

The major consideration that we have ignored so far in the theoretical development but will handle in this section is that certainly the domain for y for \( {\overset{\sim }{P}}_{\theta_i}\left({f}_x\right)(y) \) is finite, say [−m,m]. As before with x, for all practical purposes, we will now additionally consider the y-axis to be discrete, and for each projection \( {P}_{\theta_i}(f)\left(x,y\right) \), the domain is given as:

$$ D=\left\{\left(x,y\right):\kern1em x=1,2,\dots, n,\kern1em y=-m,-m+1,\dots, m\right\}. $$

We chose the indexing for y symmetrically for convenience in the center-of-mass computations so that the center of the projections is along the modeled axis of rotation at y=0. Computing \( {t}_x^{\theta_i} \) now becomes:

$$ {t}_x^{\theta_i}=\frac{1}{M_x}\sum_{y=-m}^m{\overset{\sim }{P}}_{\theta_i}\left(x,y\right)y. $$

The first issue is that M x may vary through the tilt series for each cross-section; in particular, since the domain for y is limited, there may be some observable mass moving in and out of the field of view after rotation and projection, as we demonstrated in Figure 1. This is again why it’s important that we choose the alignment to be considered over many projected cross-sections.

To handle this issue, we multiply \( {\overset{\sim }{P}}_{\theta_i}(f)\left(x,y\right) \) by a window function, \( {\omega}_{\theta_i}\left(x,y\right) \), in the computation of \( {t}_x^{\theta_i} \) in order to alleviate some of this transition of mass in and out of the frame. The window function allows for the balance of the total mass within each projection. We choose our window functions to satisfy the following properties:

  1. (i)

    \( 0\le {\omega}_{\theta_i}\left(x,y\right)\le 1; \)

  2. (ii)

    \( M=\sum_{x=1}^n\sum_{y=-m}^m{P}_{\theta_i}(f)\left(x,y\right){\omega}_{\theta_i}\left(x,y\right) \), for i=1,2,…,k;

  3. (iii)

    \( {\omega}_{\theta_i}\left(x,y\right)\le {\omega}_{\theta_i}\left(x,y+1\right)\kern1em \mathrm{if}\kern1em y<0, \)

    \( {\omega}_{\theta_i}\left(x,y\right)\ge {\omega}_{\theta_i}\left(x,y+1\right)\kern1em \mathrm{if}\kern1em y\ge 0; \)

  4. (iiii)

    \( {\omega}_{\theta_i}\left(x,y\right)={\omega}_{\theta_i}\left(x+1,y\right) \), for x=1,2,…,n−1.

The first property simply emphasizes that multiplication of \( {\overset{\sim }{P}}_{\theta_i}(f) \) by \( {\omega}_{\theta_i} \) reweighs the projection values in order to dampen the introduction of new mass in to the frames. The second property then tells us that this dampening of the values of \( {P}_{\theta_i}(f) \) by multiplication of \( {\omega}_{\theta_i} \) yields the same total mass in each projection. Finally, properties (iii) and (iiii) describe how this dampening should be done. Property (iii) says that the window function decreases as the function moves away from the y-axis. This is because new mass would be introduced along the edge of the plane of view, so that we dampen these values more significantly. Property (iiii) is an additional property to help us better characterize \( {\omega}_{\theta_i} \) in a simple manner and simply says that we place the same weight for each cross-section x. One could remove property (iiii) and change property (ii) so that instead the mass M x is fixed for each cross-section of each projection. This could potentially cause bias in the alignment of the cross-sections, especially ones with considerable noise, and it would require much greater computational time to determine a window for each cross-section of each projection.

After the windowing function is determined, we then compute the center of mass for each projected cross-section \( {t}_{x_j} \), for j=1,2,…,n as:

$$ {\overset{\sim }{t}}_{x_j}^{\theta_i}=\frac{1}{M_{x_j}^{\theta_i}}\sum_{y=-m}^m{\overset{\sim }{P}}_{\theta_i}\left({f}_{x_j}\right)(y){\omega}_{\theta_i}(y)y\kern1em \mathrm{and}\kern1em {t}_{x_j}=\left(\begin{array}{c}{\overset{\sim }{t}}_{x_j}^{\theta_1}\\ {}{\overset{\sim }{t}}_{x_j}^{\theta_2}\\ {}\vdots \\ {}{\overset{\sim }{t}}_{x_j}^{\theta_k}\end{array}\right), $$

and solve a variant of Equation 6. The variation is that we only choose to minimize only a subset of the cross-sections, say T⊂{1,2,…,n}. This subset is chosen so that the selected cross-sections have a significant quantity of mass in each projection so that introduction of new mass along the edges has considerably less effect on the center of mass of this projected cross-section area. In addition, we only choose those in which the observable total mass within that cross-section varies little throughout all projections, to again avoid the cross-sections with large transition of mass.

More precisely, we pick the cross-sections in which the ratio of the average observed mass through the projections to the variance of the mass in the projections is above some specified tolerance. This tolerance can be chosen based upon quality of the data. Finally, the minimization for determining the shifts becomes:

$$ \underset{y_{\varTheta }}{ \min}\left(\sum_{j\in T}\parallel \left({\varTheta}^{+}\varTheta -I\right)\left({y}_{\varTheta }+{t}_{x_j}\right){\parallel}_2^2\right), $$
((8))

which can again be converted into a standard least squares minimization problem as done in Equation 7. We summarize the method with the simple schematic shown in Figure 3.

Figure 3
figure 3

The general workflow of our alignment approach.

Full size image

Reconstruction method

After the alignment, for the reconstruction, we use a compressed sensing approach by total variation (TV) minimization [16]. These methods have recently been gaining popularity for electron tomographic reconstructions [17-19]. In order to briefly describe the method, let us denote the 3-D reconstructed approximation of f by \( g={\left\{{g}_{x,y,z}\right\}}_{x,y,z=1}^N \), where for simplicity we now let our discrete 3-D domain be:

$$ D=\left\{\left(x,y,z\right):x,y,z\in \left\{1,2,\dots, N\right\}\right\}. $$

Most reconstruction methods are then designed so that numerical reprojection of g agrees with the experimental projections \( {P}_{\theta_i}(f) \), for i=1,2,…,k. In particular, reconstruction techniques typically minimize the distance between the projections of g and the experimental projections, sometimes called the projection error. This projection error can be expressed as:

$$ \sum_{i=1}^k\mathrm{dist}{\left({P}_{\theta_i}(f),{P}_{\theta_i}(g)\right)}^2\kern0.3em =\kern0.3em \sum_{i=1}^k\sum_{x,y=1}^N{\left({P}_{\theta_i}(f)\left(x,y\right)\kern0.3em -\kern0.3em {P}_{\theta_i}(g)\left(x,y\right)\right)}^2. $$
((9))

However, simple minimization of the projection error does not necessarily produce optimal results in the presence of noise. Therefore, methods, such as TV minimization, additionally apply regularization conditions on the reconstruction. In the case that our sample consists of homogeneous materials and relatively smooth surfaces, compressive-sensing theory allows us to assume that the reconstruction should have a small total variation norm, given by:

$$ \begin{array}{lll}\parallel g{\parallel}_{TV}& \kern0.3em =\kern0.60em \sum_{x,y=1}^N\sum_{z=1}^{N-1}\kern0.3em \left|{g}_{x,y,z+1}\kern0.3em -\kern0.3em {g}_{x,y,z}\right|\kern0.3em +\kern0.60em \sum_{x,z=1}^N\sum_{y=1}^{N-1}\left|{g}_{x,y+1,z}\kern0.3em -\kern0.3em {g}_{x,y,z}\right|\kern2em & \kern2em \\ {}& \kern1em +\sum_{y,z=1}^N\sum_{x=1}^{N-1}\left|{g}_{x+1,y,z}-{g}_{x,y,z}\right|.\kern2em & \kern2em \end{array} $$

With this in mind, we would like for Equation 9 to be relatively small, while also applying a penalty on ∥g∥ TV for noise reduction, so that our method solves:

$$ \underset{g}{ \min}\left\{\right.\parallel g{\parallel}_{TV}+\lambda \sum_{i=1}^k\mathrm{dist}{\left({P}_{\theta_i}(f),{P}_{\theta_i}(g)\right)}^2\left\}\right.. $$
((10))

Results and discussion

We will give the results for experimental and simulation data. We compare the reconstructions from alignment using cross-correlation and our center-of-mass technique, while also demonstrating the advantage of using many slices for the center-of-mass alignment, as opposed to just one center-of-mass calculation.

Experimental results

For the experimental data, we have an alumina particle sitting on a holey carbon grid. The sample was prepared by grinding the alumina spheres into powder. A suspension of the powder is prepared in ethanol and sonicated for 5 min. The suspension was then added drop-wise over the lacey carbon film supported on 200 mesh Cu TEM grids (Structure Probe, Inc., West Chester, PA, USA) and dried at room temperature. The sample is analyzed using the FEI Titan 80-300 Scanning Transmission Electron Microscope equipped with a spherical-aberration probe-corrector (CEOS GmbH, Heidelberg, Germany) operating at 200 kV. The images were collected using the high-angle annular detector with the camera length of 195 mm and at 80,000 X magnification. The acquisition time was set to 15 s over an image area of 1024 X 1024 pixels resulting in a pixel size of 0.2411 nm. The tilt series is collected using linear tilt scheme continuously from -70° to +70°with tilt increments of 2°. Dynamic STEM focus function is used to compensate for change in focus across the image. The projection of the sample at 30°degrees is shown in Figure 2, and the aligned projections are shown in a video in Additional file 1.

Total variation minimization is valid for this data set, as the alumina particle and the carbon grid are known to be uniform in density. In addition, regularization of the reconstruction with TV minimization is critical to the quality of the results due to the low-dose sampling conditions necessary for acquisition of the projections due to beam sensitivity of the material. The reconstructed images from cross-correlation and our alignment methods are shown in Figure 4. While the overall particle morphologies are similar, the reconstruction resulting from our alignment displays much more uniform densities and clearer particle structures. This will result in more confident segmentation and characterization of the reconstructed particle, which is crucial to the interpretations of the experiment. In the 3-D images (visualized using tomviz software [20]), the overall structures appear similar. However, less rigid particle structure is recovered with the cross-correlation alignment, as the red glow around the particle demonstrates blurring from the main particle structure to a lower gray level represented by red in the colormap. In Figure 5, we plotted the centers of mass, t x , for two cross-sections. Plotted together with t x are least squares solutions of the center of mass, \( \left({c}_x^y,{c}_x^z\right) \), based Equation 3 given the computed t x . It is evident that our method finds a nearly viable path for the motion of the center of mass, as we set out to do. On the other hand, the alignment from cross-correlation clearly fails to do so, resulting in low-resolution reconstructions.

Figure 4
figure 4

Reconstructions from cross-correlation and our alignment approach. (a-c) Cross-section images of the 3-D volume from cross-correlation alignment. (d-f) Same cross-sections shown as (a-c) resulting from implementing our alignment method. (g, h) 3-D volume renders of the two reconstructions from cross-correlation alignment (g) and our alignment method (h). The scale bar in (a) is valid for (a-f), and the scale bar in (g) is valid for (g) and (h). It is apparent from these images that more blurring is present from the cross-correlation as a result of misalignment.

Full size image
Figure 5
figure 5

Location of the centers of mass of single cross-sections for each projection angle (blue) and the least squares solutions to fit the viable paths (red) given by Equation 3. The results from cross-correlation for two cross-sections are given in (a, c), and the results from our alignment method for the same cross-sections are shown in (b, d).

Full size image

In Figure 6, additional results are given using the alignment method described in [15]. Again the 3-D visual comparison of the reconstructions show that our alignment has produced a more rigid structure, as there is less red glow from the main particle but less significant than the results from cross-correlation. Similarly, the images in Figure 6c,d,e,f of the 2-D cross-sections show a more rigid structure and less noisy artifacts due to misalignment. The plots in Figure 6 give a quantitative comparison of the alignment approaches. In Figure 6g,h, the location of the global projected center of mass along the y-axis is shown for the two methods. The plot in Figure 6g shows the only consideration for the originally proposed center-of-mass alignment, as the center of mass in the projections along the y-axis is shifted to the tilt axis. With pixelation of the images, there is still a small negligible distance (less than half a pixel) between the center of mass in each projection and the tilt axis. The location of this center of mass resulting from our approach is shown in Figure 6h and does not necessarily follow a viable path, because we choose a different minimization and allow our approach to avoid problematic cross-sections. In Figure 6i,j, the path of the projected center of mass is shown for a single cross-section for the two alignment methods, where, for this cross-section, our methods demonstrate a viable path and the approach based on the single global center of mass does not. Inevitably, our method produces better reconstruction results, demonstrating that a more sophisticated alignment approach should be taken for dependable results as we have done, taking into account not one single data point but rather all cross-sections as unique data points. The resulting segmentation of the alumina particle is shown in 3-D in a video in the Additional file 2.

Figure 6
figure 6

Results from alignment in [ 15 ] and our approach. (a, b) Images of 3-D volume rendering of the reconstructions from the [15] (a) and our method (b). (c, e) 2-D cross-sections images of the 3-D reconstruction shown in (a). (d, f) Images of corresponding 2-D cross-sections of the 3-D reconstruction shown in (b). (g, h) Plots of the path of the projected global center of mass along the y-axis for the two alignment methods. (i, j) Plots of the path of a center of mass along a single cross-section of the projections for the two alignment methods.

Full size image

Simulation results

As a numerical test, we reconstructed simulated data by projecting a discrete 3-D volume with binary intensities at the same tilt angles as the experimental data: a maximum tilt range of ± 70 °in 2°-angle increments. We align the projection images according to the various alignment methods, and each realigned set of projections is reconstructed again using TV minimization. The results from the simulations are shown in Figure 7. The total projected volume shows little variation depending on the tilt angle, with the exception of a small mass appearing in the projection range at high-tilt angles. This is indicated in the projection images shown in Figure 7a,b, where, in Figure 7a, the bundle of mass is located towards the upper right of the projection image, and in Figure 7b, this bundle of mass has nearly moved completely out of the projecting range. With the special example we have here, this small transition of mass will significantly affect the results of an alignment approach such as in [15]. This is very clear from the resulting blurry reconstruction in Figure 7e that does not resemble a binary reconstruction. In addition, it can be seen in Figure 7d that even in this noise-free simulation cross-correlation also produces very poor results simply because the model is not appropriate. In Figure 7c, it is seen that our center-of-mass approach still yields optimal results displaying a near binary reconstruction image that almost completely resembles the original phantom not presented in the figure. The adaptability of our method to choose only the appropriate cross-sections with little variability of mass is clearly advantageous as demonstrated in these simulations.

Figure 7
figure 7

Tomographic simulations with a binary 3-D phantom. (a, b) Projection images of the phantom tilted about the axis at -50° and -32°, respectively. (c-e) 2-D cross-section of the reconstructed phantom from registering the data with different alignment techniques. (c) Result from our center of alignment method. (d) Result from cross-correlation. (e) Result from originally proposed center-of-mass technique.

Full size image

Conclusions

Our method has a sound physical basis: the movement of the center of mass in each cross-section. By selecting shifts for individual tilt-series images that globally lead to physically plausible motions for the centers of mass of many cross-sections, our method effectively utilizes the assumption that the sample object is rigid to improve the alignment and the resolution of the final reconstruction. We have shown that conventional alignment procedures, which shift the global center of mass to the origin, may not produce physically plausible motions in other cross-sections. We have generalized these methods in a computationally feasible manner that can be easily be incorporated into electron tomography workflows. We have demonstrated the significance of such consistency between cross-sections and the effectiveness of the presented method by improving the resolution of 3-D reconstructions of simulated and actual data.