Automatic alignment for three-dimensional tomographic reconstruction

TL;DR

This paper introduces a general algorithmic framework to automatically correct for misalignments in 3D tomographic reconstruction, improving image quality when calibration errors are present.

Contribution

It presents a novel combined approach for parameter retrieval and algebraic reconstruction, applicable to noisy and incomplete measurement data.

Findings

Effective correction of alignment errors demonstrated on simulated data.

Improved reconstruction quality shown on electron tomography dataset.

Framework adaptable to various sampling and error conditions.

Abstract

In tomographic reconstruction, the goal is to reconstruct an unknown object from a collection of line integrals. Given a complete sampling of such line integrals for various angles and directions, explicit inverse formulas exist to reconstruct the object. Given noisy and incomplete measurements, the inverse problem is typically solved through a regularized least-squares approach. A challenge for both approaches is that in practice the exact directions and offsets of the x-rays are only known approximately due to, e.g. calibration errors. Such errors lead to artifacts in the reconstructed image. In the case of sufficient sampling and geometrically simple misalignment, the measurements can be corrected by exploiting so-called consistency conditions. In other cases, such conditions may not apply and we have to solve an additional inverse problem to retrieve the angles and shifts. In this paper we propose a general algorithmic framework for retrieving these parameters in conjunction with an algebraic reconstruction technique. The proposed approach is illustrated by numerical examples for both simulated data and an electron tomography dataset.

Equations

$$p(\mathbf{s},\mathbf{\eta })={\int }_{\mathbb{R}}\mathrm{d}tu(\mathbf{s}+t\mathbf{\eta }),$$

$$\hat{p}(\mathbf{\xi },\mathbf{\eta })=(2\pi {)}^{-\frac{3}{2}}{\int }_{\mathbb{R}}\mathrm{d}t\hat{u}(\mathbf{\xi })\exp \left(-ıt(\mathbf{\xi }\cdot \mathbf{\eta })\right),$$

$$W\mathbf{u}=\mathbf{p},$$

$$W(\mathbf{a})\mathbf{u}=\mathbf{p},$$

$$\underset{\mathbf{a},\mathbf{u}}{\min }f(\mathbf{a},\mathbf{u})=\frac{1}{2}\Vert W(\mathbf{a})\mathbf{u}-\mathbf{p}{\Vert }_2^2+h(\mathbf{u}),$$

$$\nabla f(\mathbf{a},\mathbf{u})=\left(\begin{array}{c}G(\mathbf{a},\mathbf{u}{)}^T\left(W(\mathbf{a})\mathbf{u}-\mathbf{p}\right)\\ W(\mathbf{a}{)}^T\left(W(\mathbf{a})\mathbf{u}-\mathbf{p}\right)+\nabla h(\mathbf{u})\end{array}\right),$$

$${\nabla }^{2}f(\mathbf{a},\mathbf{u})=\left(\begin{array}{cc}G(\mathbf{a},\mathbf{u}{)}^{T}G(\mathbf{a},\mathbf{u})+R(\mathbf{a},\mathbf{u},\mathbf{r})& G(\mathbf{a},\mathbf{u}{)}^{T}W(\mathbf{a})+H(\mathbf{a},\mathbf{r}{)}^T\\ W(\mathbf{a}{)}^{T}G(\mathbf{a},\mathbf{u})+H(\mathbf{a},\mathbf{r})& W(\mathbf{a}{)}^{T}W(\mathbf{a})+{\nabla }^{2}h(\mathbf{u})\end{array}\right),$$

$$\underset{\mathbf{a}}{\min }\overline{f}(\mathbf{a}),$$

$$\nabla \overline{f}(\mathbf{a})$$

$${\nabla }^2\overline{f}(\mathbf{a})={\nabla }^{2}f(\mathbf{a},\overline{\mathbf{u}})\backslash (W^{T}W+{\nabla }^{2}h(\overline{\mathbf{u}})),$$

$$\kappa \left({\nabla }^2\overline{f}\right)\le \kappa \left({\nabla }^{2}f\right),$$

$$\Vert {\nabla }_{\mathbf{u}}f(\mathbf{a},{\overline{\mathbf{u}}}_{ϵ}){\Vert }_2\le ϵ,$$

$${\overline{f}}_{ϵ}(\mathbf{a})=f(\mathbf{a},{\overline{\mathbf{u}}}_{ϵ}).$$

$$\nabla {\overline{f}}_{ϵ}(\mathbf{a})={\nabla }_{\mathbf{a}}f(\mathbf{a},{\overline{\mathbf{u}}}_{ϵ}).$$

$$\Vert \nabla {\overline{f}}_{ϵ}(\mathbf{a})-\nabla \overline{f}(\mathbf{a}){\Vert }_2\le Cϵ,$$

$$\nabla {\overline{f}}_{ϵ}(\mathbf{a})=\nabla \overline{f}(\mathbf{a})+{\nabla }_{\mathbf{a},\mathbf{u}}^{2}f{\left({\nabla }_{\mathbf{u},\mathbf{u}}^{2}f\right)}^{-1}{\nabla }_{\mathbf{u}}f(\mathbf{a},{\overline{\mathbf{u}}}_{ϵ}),$$

$$\Vert \nabla {\overline{f}}_{ϵ}(\mathbf{a})-\nabla \overline{f}(\mathbf{a}){\Vert }_2\le Cϵ,$$

$$\underset{\mathbf{a}}{\min }\overline{f}(\mathbf{a})+g(\mathbf{a}),$$

$$\frac{\mu }{2}\Vert \mathbf{a}-{\mathbf{a}}^{\prime }{\Vert }_2^2\le \overline{f}(\mathbf{a})+g(\mathbf{a})-\overline{f}({\mathbf{a}}^{\prime })-g({\mathbf{a}}^{\prime })$$

$$\underset{\mathbf{u}}{\min }f(\mathbf{u},{\mathbf{a}}^{(k)})$$

$$\left(\overline{f}({\mathbf{a}}^{(k+1)})-\overline{f}(\overline{\mathbf{a}})\right)\le (1-\mu /L)\left(\overline{f}({\mathbf{a}}^{(k)})-\overline{f}(\overline{\mathbf{a}})\right)+\frac{C^2{ϵ}_k^2}{2L},$$

$${\mathsf{prox}}_{\gamma g}(\mathbf{x})=\underset{\mathbf{a}}{\mathrm{argmin}}\frac{1}{2}\Vert \mathbf{a}-\mathbf{x}{\Vert }_2^2+\gamma g(\mathbf{a})$$

$$W_{\text{roi}}(\mathbf{a}):=T_{\text{roi}}W(\mathbf{a}).$$

$$p(\mathbf{s},\mathbf{\eta })={\int }_{\mathbb{R}}\mathrm{d}tu(\mathbf{s}+t\mathbf{\eta }),$$

$$\mathbf{p}=W(\mathbf{a})\mathbf{u}.$$

$$W(\mathbf{a})=\left(\begin{array}{c}{W}_1({\mathbf{a}}_1)\\ W_2({\mathbf{a}}_2)\\ ⋮\\ W_{N_{\text{proj}}}({\mathbf{a}}_{N_{\text{proj}}})\end{array}\right).$$

$$G(\mathbf{a},\mathbf{u})=\frac{\partial W(\mathbf{a})\mathbf{u}}{\partial \mathbf{a}}=\left(\begin{array}{cccc}\frac{\partial W_1}{\partial {\mathbf{a}}_1}({\mathbf{a}}_1)\mathbf{u}& & & 0\\ & \frac{\partial W_2}{\partial {\mathbf{a}}_2}({\mathbf{a}}_2)\mathbf{u}& & \\ & & \ddots & \\ 0& & & \frac{\partial W_{N_{\text{proj}}}}{\partial {\mathbf{a}}_{N_{\text{proj}}}}({\mathbf{a}}_{N_{\text{proj}}})\mathbf{u}\end{array}\right).$$

$$W_i({\mathbf{a}}_i)=S_{z}R({\mathbf{a}}_i),$$

$$R({\mathbf{a}}_i)(u{)}_{{\mathbf{x}}_{\text{s}}}=\tilde{u}(A_R({\mathbf{a}}_i{)}^{-1}({\mathbf{x}}_{\text{s}}))$$

$$\tilde{u}(\mathbf{x})$$

Full text

Published as: Inverse Problems 34(2):024004, 2018

1 Introduction

In tomographic reconstruction, the aim is to recover an unknown object described by a function $u:\mathbb{R}^{3}\rightarrow\mathbb{R}_{+}$ with compact support from a finite number of samples of its X-ray transform

[TABLE]

where $\mathbf{s}\in\mathbb{R}^{3}$ and $\boldsymbol{\eta}\in\mathbb{S}^{2}$ determine the origin and direction of the X-ray (with $\mathbb{S}^{2}=\{\mathbf{x}\in\mathbb{R}^{3}\,\,|\,\,\|\mathbf{x}\|_{2}=1\}$ denoting the three-dimensional unit sphere). The samples of the X-ray transform are typically obtained by a scanning device, such as a CT scanner, where a source and detector traverse a certain trajectory around the object, which defines the set of sample points $\{(\mathbf{s}_{i},\boldsymbol{\eta}_{i})\}_{i=1}^{m}$. A few common setups are illustrated in figure 1. Here, a number of individiual ray-trajectories are illustrated as well as the axis and range of rotation of the object.

The inverse problem consists of reconstructing $u$ from samples $\{p(\mathbf{s}_{i},\boldsymbol{\eta}_{i})\}_{i=1}^{m}$ of $p$. Upon Fourier transforming $p$ with respect to its first argument we get

[TABLE]

which forms the basis for many analytic reconstruction formulas. Indeed, for the parallel-beam geometry we have $\boldsymbol{\xi}\cdot\boldsymbol{\eta}=0$ and we retrieve the well-known Fourier Slice Theorem, which states that the Fourier transform of $p$ along $\mathbf{s}$ gives the Fourier transform of $u$ along a slice perpendicular to $\boldsymbol{\eta}$ [25]. For some specific acquisition geometries, such as parallel-beam and fan-beam, explicit inverse formulas have been derived [29]. These typically involve a filtering of the projection data followed by a backprojection step (filtered backprojection) or interpolation followed by an inverse Fourier transform. We can think of these as invert-then-discretize approaches to solving the inverse problem. Being derived from a continuous model, these formulas only give an accurate reconstruction for fully and regularly sampled data.

On the contrary, various non-regular and under-sampled scenarios occur regularly in practice and can often be conveniently approached by first discretizing and then inverting the resulting system of linear equations. Such a discretize-then-invert approach is often referred to as algebraic reconstruction in the literature [18, 22].

Upon expressing $u$ in terms of an appropriate basis with coefficients $\{u_{i}\}_{i=1}^{n}$ and approximating the integral, the reconstruction problem can be posed as system of linear equations

[TABLE]

where $W\in\mathbb{R}^{m\times n}$ is the projection matrix whose elements depend on the acquisition parameters, $\mathbf{u}\in\mathbb{R}^{n}$ represents the object and $\mathbf{p}\in\mathbb{R}^{m}$ are the projection data for given rays characterized by $(\mathbf{s}_{i},\boldsymbol{\eta}_{i})$, i.e., $p_{i}=p(\mathbf{s}_{i},\boldsymbol{\eta}_{i})$. This moves the reconstruction problem into completely different realms of mathematics. Results and algorithms from the field of convex analysis and image processing can be brought to bear to include various regularization terms such as Total Variation to define a unique regularized solution. Although exact reconstruction can no longer be guaranteed with such approaches, they are often used in practice and may yield accurate results even in the case of highly corrupted or incomplete data [35, 19].

When dealing with tomographic problems in practice, a secondary type of inverse problem must often be solved as part of the image reconstruction process. In the standard models for tomographic reconstruction, the geometry of the imaging system (position and 3D orientation of sources and detectors) are assumed to be perfectly known and are a fixed part of the inverse problem that requires no further attention. In practice, however, the corresponding acquisition parameters, which determine the trajectories of the X-rays through the unknown object, may be known a-priori only with limited accuracy. Recovering these parameters then becomes an inverse problem of itself. An extreme case is single-particle analysis based on images obtained by cryo-electron microscopy, where the orientation of the object with respect to the incident X-rays is completely unknown [15]. Moreover, even if the acquisition parameters are prescribed by the setup, calibration errors or (rigid) movement of the object itself may result in a mismatch between the assumed ideal geometry and the effective geometry used to record the projections. Reconstruction from such distorted datasets introduces artifacts in the recovered object and leads to a loss of resolution. Estimating the acquisition parameters to correct such alignment errors is therefore of crucial importance in many practical applications [24, 31, 28].

A relatively straightforward and practical approach to deal with the estimation of the alignment parameters is the introduction of markers attached to the scanned object, each of which can be tracked through the series of projection images. The motion of these markers (often small beads made of dense material) can then be related to their 3D position, determined from the series of observed 2D positions in the projection images. The strength of this approach lies in its ability to accurately recover many different geometrical parameters, including not only shifts of the projection images, but also detector tilts and variations in the position of the X-ray source. A significant drawback is the need to prepare the object before the scan to apply the markers, and also the effects that the markers themselves can have on the reconstructed images: as marker particles have to be dense structures to show up in the individual projection images, they typically lead to streak artefacts in the reconstructions [37].

A common approach to the alignment problem that does not require the presence of markers is cross-correlation. The idea is that two projections with a small angular separation will be highly similar in structure, except for a geometric stretching of the data. By incorporating such geometric assumptions in the alignment procedure and correlating consecutive images to each other sequentially, some shift and rotation parameters can be estimated in a heuristic way, but with limited accuracy. Out-of-plane movements on the other hand, such as detector tilts and variations in the position of the X-ray source, cannot be retrieved with such approaches.

The category of alignment methods that do not require the use of markers or feature-detection in the projections is mainly based on exploiting the redundancy in the tomographic data. Similar to the development of analytical and algebraic methods for solving the reconstruction problem, the alignment problem can also be approached along these two different paths. On the analytical side, consistency conditions can be formulated for various source trajectories, including the well-known Helgason-Ludwig consistency conditions [21] for parallel-beam illumination. These prescribe sets of equations that must be satisfied by any valid set of tomographic projections. For example, the first order Helgason-Ludwig condition prescribes that the center of mass of the scanned object must project onto the center of mass of each of the measured projections, which can be used to estimate unwanted shifts of the detector that vary across the range of angles. A drawback of such methods is that they require full-field, i.e. non-truncated projection data. Hence, they are - strictly speaking - not applicable to region-of-interest tomography, although necessary conditions [11] and heuristic adaptations [33] have been proposed recently. Local consistency conditions can also be defined in some cases, e.g. based on John’s equation for cone-beam illumination [38]. Although these conditions can be evaluated in a region-of-interest tomography setting, they instead rely on the assumption that the angular sampling is sufficient to allow for the numerical computation of derivatives that approximate a continuum of projection angles. Recently, a new type of consistency conditions has been introduced in [1] based on epipolar geometry that is more general and may lift some of these constraints, but its full relation to the classical consistency conditions remains to be investigated.

Alternatively, estimation algorithms for the alignment parameters have been developed based on minimizing the inconsistencies between the forward model applied to the reconstructed image and the measured projection data. In [9], a quasi-newton algorithm is used to simultaneously compute the reconstructed image and the corresponding alignment parameters, using a limited-memory Broyden method to approximate the Hessian term. A similar approach is followed in [6], where a Levenberg-Marquardt method is used instead. In [31], a scheme for alignment of X-ray images is proposed that alternates between image reconstruction steps (using fixed alignment parameters) and alignment steps (using a fixed reconstruction). In addition to the consistency with the measured data, [23] employs a sharpness criterion in the image domain as part of the optimization process. This broad class of methods, known as projection matching, follow a discretize-then-invert approach and are much more general than the analytical approach. Similar to the algebraic reconstruction methods, they allow for the incorporation of prior knowledge about the object within the cost function for the parameter estimation and they can deal effectively with non-standard geometrical setups in principle, although existing implementations are often limited to correcting in-plane misalignment.

In this paper, we consider a rather general formulation of the projection matching approach and investigate various mathematical and numerical aspects of different solution strategies. Although a broad variety of projection matching algorithms have been described in the literature, they are typically proposed in relation to a specific application, where a single algorithm variant is proposed evaluated based on empirical results only. Instead, we focus on the underlying numerical mechanisms and outline the different choices that can be made in splitting the problem into subproblems.

We consider the setting where alignment errors can be modelled by three-dimensional rigid motion of the object. Accordingly, each projection has its own set of up to six alignment parameters (since any rigid motion in three dimensions can be described by three shifts and three rotations). We assume that for each projection, the corresponding alignment parameters are approximately known. This is true for most tomographic imaging setups, where projection images are acquired from known directions and only small deviations from the known geometry need to be determined, but not in the field of single-particle-analysis for example, where projections can be randomly distributed with completely unknown angles. In such cases, different algorithms have to be applied such as the approach in [12] based on graph laplacians. Note that the alignment parameters for each projection image are considered as independent, and thus can be used to model model rigid motion of the object relative to the fixed source-detector-system during the data acquistion.

The forward model for the tomographic alignment problem can be expressed as

[TABLE]

where $\mathbf{a}\in\mathbb{R}^{\ell}$ with $\ell=6N_{\textup{proj}}$ is a vector containing a suitable parametrization of a rigid object motion, $N_{\textup{proj}}\in\mathbb{N}$ denoting the number of projection images. We propose an algorithmic approach for joint alignment and tomographic reconstruction based on the formulation (1). This paper provides a mathematical framework that encompasses a large class of such methods, allowing for analysis in a common framework.

The paper is organized as follows. In section 2, we present the main ideas underlying our approach and provide a detailed analysis of the proposed algorithm. Some practical aspects of the implementation of the algorithm are discussed in section 3. Results on both phantom and real datasets are presented in section 4. Section 5 concludes the paper.

2 Non-linear least-squares

We pose the joint reconstruction and alignment problem as a non-linear least-squares problem

[TABLE]

where $h$ is an appropriate regularization term. This problem is related to the Total Least-Squares (TLS) problem in the sense that the matrix itself is (to some extent) unknown [40]. In the traditional TLS formulation, however, the perturbation $W$ would be modelled as a matrix-valued unknown, i.e., $W+\Delta W$. Although this provides an interesting viewpoint on the alignment problem, it is outside the scope of this paper to pursue this connection further.

Throughout the paper we will make the following assumptions on $f$:

• (A1)

$W(\mathbf{a})$ has full rank for all admissable $\mathbf{a}$ or $h$ is strongly convex. It follows from these assumptions that $f$ is strongly convex in $\mathbf{u}$ for any admissible $\mathbf{a}$ so that the solution $\overline{\mathbf{u}}(\mathbf{a})$ that minimizes $f(\mathbf{a},\mathbf{u})$ is unique.

• (A2)

There exists an optimal set of alignment parameters and image, $(\overline{\mathbf{a}},\overline{\mathbf{u}})$, that minimizes $f$ and $f$ is convex in the neighborhood of this minimizer.

We note that only the first assumption is absolutely critical, but can be fulfilled easily in practice. The second assumption is necessary for some of the results that follow. For the remainder of the section we will furthermore assume that

• (A3)

$f$ is twice continuously differentiable in $\mathbf{a}$. This follows immediately if $W$ is twice continuously differentiable in $\mathbf{a}$ which requires the use of higher-order quadrature rules for approximating the projection-integrals.

• (A4)

$h$ is twice continuously differentiable.

We will discuss to what extent these latter assumptions can be dropped in section 2.3.

The gradient and Hessian of $f$ are given by

[TABLE]

and

[TABLE]

where $G(\mathbf{a},\mathbf{u})$ is the Jacobian of $\mathbf{a}\mapsto W(\mathbf{a})\mathbf{u}$ with columns $(\partial W/\partial a_{i})\mathbf{u}$, $H(\mathbf{a},\mathbf{v})$ is the Jacobian of $\mathbf{a}\mapsto W(\mathbf{a})^{T}\mathbf{v}$ with columns $\left({\partial W}/{\partial a_{i}}\right)^{T}\mathbf{v}$, $\mathbf{r}:=W(\mathbf{a})\mathbf{u}-\mathbf{p}$ and $R(\mathbf{a},\mathbf{u},\mathbf{v})\in\mathbb{R}^{\ell\times\ell}$ has entries $\mathbf{u}^{T}(\partial^{2}W(\mathbf{a})/\partial a_{i}\partial a_{j})\mathbf{v}$. Note that by (A2), the full Hessian is positive semidefinite close to the minimizer.

We could in principle view (2) as an optimization problem in $(\mathbf{a},\mathbf{u})$ and apply a gradient-based method to solve it directly. Such a brute-force joint optimization approach introduces several problems, however. Firstly, the problem can be severely ill-conditioned due to differences in sensitivity of $\mathbf{u}$ and $\mathbf{a}$. Furthermore, a joint approach makes it hard to employ advanced reconstruction techniques that have been developed for linear tomographic problems.

2.1 Variable projection

The joint estimation problem (2) has a separable structure; it is strongly convex in $\mathbf{u}$ and differentiable in $\mathbf{a}$ (by assumption (A3)). Following the variable projection approach [17], we eliminate the reconstructed object by setting $\overline{\mathbf{u}}(\mathbf{a}):=\operatornamewithlimits{argmin}_{\mathbf{u}}f(\mathbf{a},\mathbf{u})$. The problem now reduces to

[TABLE]

where $\overline{f}(\mathbf{a})=f(\mathbf{a},\overline{\mathbf{u}}(\mathbf{a}))$.

We can now state the following results:

Theorem 2.1**.**

Let $\overline{f}(\mathbf{a})=\operatornamewithlimits{argmin}_{\mathbf{u}}f(\mathbf{a},\mathbf{u})$. The gradient of $\overline{f}$ is given by

[TABLE]

i.e. it is simply the gradient of the full objective, evaluated at the optimal $\mathbf{u}$. The Hessian of $\overline{f}(\mathbf{a})$ is given by

[TABLE]

i.e., the Schur complement of $(W^{T}W+\nabla^{2}h)$ in the full Hessian of $f$.

Proof.

The gradient-expression follows directly from [5, Thm. 2] which requires assumption A1, while the Schur-complement expression for the Hessian was presented by [2]. ∎

Remark 2.1**.**

By Theorem 2.1 and assumption (A2) we find that $\overline{f}$ is locally convex as well, according to the properties of the Schur complement. In addition, we find that a (local) minimum $\overline{\mathbf{a}}$ of $\overline{f}$ together with the corresponding $\overline{\mathbf{u}}$ is a (local) minimum of $f$.

The reduced problem (5) has a much smaller dimensionality and is expected to be less ill-conditioned. Indeed, we find the following relation between the condition numbers of the Hessians of $f$ and $\overline{f}$.

Theorem 2.2**.**

Under assumptions A1-A4 we find the following relation between the condition numbers of the full and reduced Hessians locally around a minimizer:

[TABLE]

where $\kappa$ denotes the condition number of the corresponding matrix.

Proof.

Following [36, Thm. 5], we find that the eigenvalues of $\nabla^{2}\overline{f}$ interlace those of $\nabla^{2}f$. The result immediately follows. ∎

A more detailed analysis of the expected convergence behaviour when minimizing separable functions like $f$ and its reduced variant, $\overline{f}$, with Newton-like methods is presented by [2].

2.2 Inexactness

In practice, it is not feasible to compute $\overline{\mathbf{u}}$ such that $\nabla_{\mathbf{u}}f(\mathbf{a},\overline{\mathbf{u}})=0$ exactly. Rather, let us assume we have an approximation $\overline{\mathbf{u}}_{\epsilon}$ such that

[TABLE]

and define

[TABLE]

Although we could in principle compute the exact gradient of $\overline{f}_{\epsilon}$ by taking into account the finite accuracy of the reconstruction, this would require us to differentiate the algorithm used to obtain $\overline{\mathbf{u}}_{\epsilon}$. This is not very attractive from a practical point-of-view. Instead, we approximate the gradient as $\nabla\overline{f}_{\epsilon}(\mathbf{a})\approx\nabla\overline{f}(\mathbf{a})=\nabla_{\mathbf{a}}\overline{f}(\mathbf{a},\overline{\mathbf{u}}_{\epsilon})$ and analyze the error we introduce. Under assumptions (A1)-(A4), we obtain the following result.

Theorem 2.3**.**

Given $\overline{\mathbf{u}}_{\epsilon}$ satisfying $\|\nabla_{\mathbf{u}}f(\mathbf{a},\overline{\mathbf{u}}_{\epsilon})\|_{2}\leq\epsilon$, define the approximate gradient

[TABLE]

We then have

[TABLE]

for some $C>0$.

Proof.

The approximate gradient is readily expressed as

[TABLE]

from which we find

[TABLE]

where $C=\|\nabla^{2}_{\mathbf{a},\mathbf{u}}f\left(\nabla^{2}_{\mathbf{u},\mathbf{u}}f\right)^{-1}\|_{2}$. ∎

Gradient-descent algorithms with approximate derivatives have been studied extensively, and there are convergence guarantees to a local minimum of $\overline{f}$ as long as $\epsilon\downarrow 0$ or to a point within distance $\mathcal{O}(\epsilon)$ of the local minimum for a fixed $\epsilon$ [13, 16].

2.3 Extensions to less-differentiable cases

We can drop the requirement for the existence of second derivatives (A3) if we instead assume that the gradient of $f$ w.r.t. $(\mathbf{a},\mathbf{u})$ is Lipschitz continuous. The resulting reduced objective $\overline{f}$ will have a Lipschitz continuous gradient as well. Results analogous to Theorems 2.2 and 2.3 can then be derived using the local convexity of $f$ (A2).

Of special interest are non-differentiable regularization terms on $\mathbf{u}$ such as Total Variation (TV) or bound-constraints. We can readily drop assumption (A4) because the assumption of strong convexity of $f$ in $\mathbf{u}$ (A1) is sufficient to guarantee differentiability of $\overline{f}(\mathbf{a})$, see [32, Thm. 9.18].

In many practical settings, the evaluation of $W$ is implemented via (multi-)linear interpolation [41], in which case the derivatives with respect to $\mathbf{a}$ are not continuous. Under assumption (A2) that $f$ is (locally) convex, the reduced function $\overline{f}$ is (locally) convex as well [32, Thm. 2.22] so that we can use a deritative-free method for minimizing $\overline{f}(\mathbf{a})$ in this setting. A promising candidate for such problems is the implicit filtering method, which approximates the derivative using finite-differences and adapts the stepsize according to its convergence [30].

3 Algorithms

Based on the ideas presented in the previous section, we can design and analyze various joint alignment and reconstruction algorithms that aim to solve

[TABLE]

where $g(\mathbf{a})$ is an appropriate regularization term. Below, we present a few prototype gradient-based algorithms for solving (7) and discuss their convergence behaviour. To facilitate this discussion, we assume that $\overline{f}$ has a Lipschitz continuous gradient (which requires Lipschitz continuity of the gradient of $f$ w.r.t $\mathbf{a}$) and that $\overline{f}+g$ is strongly convex, i.e.,

[TABLE]

for some $\mu>0$. As noted before, we do not expect $\overline{f}$ to be globally (strongly) convex, but it is reasonable to assume $\overline{f}$ to be locally (strongly) convex to facilitate the exposition. We can typically drop the assumption of strong convexity to mere convexity at the cost of a slower convergence rate. When dropping the requirement of convexity altogether, we are left with a non-linear optimization problem without any special structure and cannot state any convergence results.

3.1 Joint reconstruction and alignment

Based on the properties of the regularization term $g$, we can design a rule to update the alignment parameters. A basic first-order algorithm for a smooth regularization term $g$ is shown in Algorithm 1.

Here and in the following, the operation $\mathsf{reconstruct}(W(\mathbf{a}^{(k)}),\mathbf{p},\epsilon_{k})$ corresponds to solving the tomographic reconstruction problem at fixed alignment parameters

[TABLE]

up to a tolerance $\epsilon_{k}$ in the gradient $\|\nabla_{\mathbf{u}}f(\mathbf{a}^{(k)},\mathbf{u}^{(k+1)})\|_{2}\leq\epsilon_{k}$. For quadratic regularizers $h$, the minimization can be performed efficiently via a conjugate-gradient method. For non-smooth convex penalties involving for example total variation or bound-constraints, Chambolle-Pock’s algorithm [10] or FISTA [4] are suitable choices.

The following theorem establishes convergence of Algorithm 1 to the minimizer, $\overline{\mathbf{a}}$, of $\overline{f}$ for a fixed steplength.

Theorem 3.1**.**

Let $\gamma^{(k)}=1/L$ and $\|\nabla_{\mathbf{u}}f(\mathbf{a}^{(k)},\mathbf{u}^{(k+1)})\|_{2}\leq\epsilon_{k}$ hold in each iteration $k$ of Algorithm 1. Then we have for some constant $C>0$

[TABLE]

where $L$ denotes the Lipschitz constant of $\nabla\overline{f}+\nabla g$ and $\mu$ is the modulus of strong convexity of $\overline{f}+g$.

Proof.

From theorem 2.3 we immediately obtain a bound of the error in the gradient. Applying [16, Lemma 2.1] gives the desired result. ∎

According to theorem 3.1, Algorithm 1 is linearly convergent under the given Lipschitz-continuity- and strong-convexity-assumptions if we decrease the tolerance $\epsilon_{k}$ at a linear rate (e.g., $\epsilon_{k}=r^{k}$ with $r<1$) [16].

When the regularization term $g$ is convex but not differentiable, we can apply a proximal gradient algorithm, as stated in Algorithm 2. Here, the proximal operator is defined as

[TABLE]

This can be used, for example, to impose bound-constraints, i.e. admissible minimum and maximum values, on the alignment parameters.

Convergence of proximal-gradient algorithms with approximate gradients has been established by [34]. Provided that $\overline{f}$ is strongly convex and that the error in the gradient decreases at a linear rate, Algorithm 2 converges linearly. If $\overline{f}$ is only convex, we obtain a sublinear convergence rate.

A variant of the proposed algorithms is stated in Algorithm 3. The method $\mathsf{align}(\mathbf{u},\mathbf{p},\delta)$ computes minimizing alignment parameters $\operatornamewithlimits{argmin}_{\mathbf{a}}f(\mathbf{a},\mathbf{u}^{(k)})$ up to some tolerance $\delta$. Convergence of Algorithm 3 can be guaranteed in case $f(\mathbf{a},\mathbf{u})$ is convex in $\mathbf{a}$ as well as $\mathbf{u}$ [7, 8]. Convergence rates for this algorithm have, to the best of our knowledge, not been established.

4 Results

The following results have been obtained with a custom implementation of a three-dimensional projection operator for parallel-beam setups that allows for full three-dimensional alignment of all projection images independently. The alignment parameters are taken to be three rotations and two shifts relative to the direction of the beam, as shown in figure 2. We refer to the rotations as in-plane (around the $z$-axis), nod/pitch (around the $x$-axis) and tomographic (around the $y$-axis) and the shifts as lateral (along the $x$-axis), axial (along the $y$-axis). For the considered parallel beam setup, shifts along the optical axis (blue axis in figure 2) do not influence the projections. The tomographic projector is implemented in Matlab, exploiting its GPU capabilities to achieve computational efficiency. Implementation Details on the projection operators, the corresponding Jacobian, the reconstruction and alignment algorithms are described in the appendix. In all of the subsequent numerical experiments, we simulate the misaligned tomographic data using trilinear interpolation, but applied the bicubic scheme from A.3 in the joint alignment and reconstruction.

4.1 Convergence experiments

We test algorithms 1 and 3 on a phantom of size $128\times 128\times 128$, shown in figure 3. We consider a parallel-beam setup with 128 tomographic projections uniformly covering the full angular range $[0;\,180^{\circ}]$. We simulate noise-free tomographic data under small perturbations of the alignment parameters that are shown in figure 3. For tomographic reconstruction, we use the Conjugate Gradient (CG) method to solve the gradient-penalized Tikhonov functional in (17). The regularization paramater is fixed to $\alpha=10^{3}$. In the experiments, the reconstruction tolerance $\epsilon$ (cf. Algorithm 1) is a relative tolerance (relative to the initial norm of the gradient) and we control the alignment tolerance $\delta$ (cf. Algorithm 3) by choosing the number of inner alignment gradient-descent iterations $N_{\text{align}}$. A large $N_{\text{align}}$ thus corresponds to a smaller $\delta$. for $N_{\text{align}}=1$, algorithms 1 and 3 are obviously identical.

Figure 4 shows the un-aligned and aligned reconstruction for $\epsilon=10^{-4}$, $N_{\text{align}}=1$, as well as the reconstruction with the true alignment parameters. The un-aligned reconstruction clearly shows blurring and streaks, which disappear in the aligned reconstruction. Moreover, the aligned reconstruction is visually indistinguishable from the one using the true alignment parameters. The fitted alignment parameters are plotted in 5. We see that the shifts are fairly well-recovered in all cases. However, the fit of the reconstruction of the rotational misalignment breaks down for the lowest CG-accuracy $\epsilon=10^{-1}$. For smaller reconstruction-tolerances, all alignment parameters except for the tomographic angle are recovered reasonably well. The reconstruction of the latter remains corrupted by low-frequency errors, which may point to an inherent non-uniqueness in the problem.

In figure 6, we show the convergence of Algorithm 3 in terms of the optimality ($\|\nabla\overline{f}\|_{2}$), the object-reconstruction error and the alignment error (the difference between the true and reconstructed alignment parameters). We note that using a larger reconstruction tolerance may lead to a premature stalling of the convergence, or even to divergence, compare results for $\epsilon=10^{-1}$. Interestingly, a large $\epsilon$ also seems to result in a slightly faster decrease of the alignment error in the first few iterations. This might be explained by the fact that fewer CG-iterations lead to a smoother recontruction, which is easier to align as argued in A.4. Furthermore, we observe that the alignment error (right column) always exhibits a slight increase after an initial phase of convergence. This is due to the emerging low-frequency errors in the fitted tomographic angle observed in figure 5. We also see that increasing the number of alignment steps ($N_{\text{align}}$) leads to slightly more optimal results (left column) but not to a higher accuracy in terms of the reconstruction and even gives a slightly larger alignment error. We may thus conclude that performing more than a single alignment update per iteration (cf. Algorithm 3) is not necessarily better than Algorithm 1.

The final aspect we wish to investigate in these experiments is the sensitivity of the iterative approach to the magnitude of the alignment errors. To this end, we use the exact same alignment errors as in the previous experiments, except that we increase their magnitude by a factor of 2, 4, 8 and 16 respectively. For the alignment, we run Algorithm 1 with $\epsilon=10^{-4}$ for a total of 20 outer iterations. The regularization parameter $\alpha$ is chosen using a heuristic described in section A.4. This heuristic is based on an estimate of the misalignment and effectively assigns a higher $\alpha$ (ranging from roughly $10^{1}$ to $10^{4}$ in our experiments) for larger misalignment. In our experience, this aspect is quite crucial since a smoother reconstruction, as induced by a larger $\alpha$, improves the initial convergence of the alignment. On the contrary, choosing $\alpha$ too low may cause the alignment to stall prematurely. In practice, one could of course use a continuation, starting from large $\alpha$ and reducing it as the alignment improves. The resulting final reconstructions (computed using $\alpha=10^{2}$ fixed) are shown in figure 7 and the corresponding alignments are shown in figure 8. Even for the most severe initial misalignment, we observe that all alignment parameters are accurately recovered, except for the tomographic angle. If the misalignment exceeds the angular sampling, it becomes exceedingly hard to find the correct projection angles as nearby projections are similar. As long as enough of the projections have been aligned correctly, the effect on the final reconstruction is not dramatic. Overall, the experiment shows that the proposed alignment procedure deteriorates gracefully as the initial misalignment increases.

4.2 Alignment of truncated projection data

So far, we have only considered settings where the object is fully captured in all of the tomographic projections. In many applications, however, only a limited region-of-interest (RoI) of an extended sample is imaged due to setup limitations or practical considerations. This setting can be modelled by composing the projection matrix $W(\mathbf{a})$ with an additional truncation matrix $T_{\textup{roi}}\in\mathbb{R}^{m_{\textup{roi}}\times m}$ that restricts the (hypothetical) full-field tomographic data to the information captured by the $m_{\textup{roi}}$ detector pixels:

[TABLE]

The limited field-of-view allows for the integrated object mass in the detection domain to change while the sample is rotated. Most importantly, this effect rules out the theoretical exactness of center-of-mass-based alignment methods, although approximate schemes, where smooth window functions dampen the in- and outflow of object mass, have been proposed [33]. On the contrary, our alignment algorithms remain applicable in the RoI-setting. We simply have to exchange the projector $W(\mathbf{a})$ by $W_{\textup{roi}}(\mathbf{a})$ and otherwise follow exactly the same theoretical and numerical approach as in the full-field case. Note that the gradient penalty in the Tikhonov-regularized reconstruction (see A.4) results in a smooth extrapolation of the recovered object outside the region-of-interest so that padding of the truncated data (see e.g. [27]) is not necessary.

As a test case, we consider a $256^{3}$-sized object composed of randomly generated ellipsoids in a cylindrical support domain and simulate 128 projections in between $0^{\circ}$ and $160^{\circ}$, truncated to the central $128^{2}$ pixels. The tomographic data is simulated under both systematic and random lateral- and axial shifts of the object, compare figure 2. Rotational misalignment is omitted in this test case. The exact phantom (central slice), the simulated RoI-tomographic data (central sinogram) and the misalignment shifts (dotted lines) are shown in figures 9(a), (b) and (c), respectively. We apply Algorithm 1 to the data, where the iterations are automatically stopped after 39 alternating reconstruction- and alignment steps when the computed correction in the alignment parameters $\mathbf{a}^{(k)}-\mathbf{a}^{(k-1)}$ drops below $0.05$ pixels in maximum norm.

The central slice of the final object reconstruction is shown in figure 9(f), the difference between the reconstructed- and exact alignment parameters is plotted in (c) (solid lines). While rapid variations of shifts $\{(s_{i,x},s_{i,y})\}$ with the tomographic angles $\theta_{i,xz}$ are almost completely eliminated by the alignment method, some slowly varying misalignment components remain. This is too not surprising since region-of-interest tomography is non-unique even without an additional alignment problem to be solved [29]. Notably, however, figure 9(c) reveals that also the bulk components of the simulated object motion are significantly damped compared to the initial misalignment. Likewise, the motion artifacts in the recovered object are almost invisible in the final reconstruction whereas quite pronounced initially, compare figures 9(e) and (f). Visually, the achieved reconstruction using the fitted alignment parameters is indeed of comparable quality as a reconstruction using the exact alignment, which is shown in figure 9(d) for comparison. The convergence plot in figure 9(g) reveals the fit of the alignment parameters to improve only very slowly, which indicates that the alignment problem is highly ill-conditioned in this setting. Indeed, the curves suggest that further improvements might be achievable by performing significantly more iterations.

4.3 Application to an electron tomography data set

For evaluation of our alignment scheme in a real-world context, we apply it to an electron tomography data set (HAADF-STEM). The sample is given by an assembly of gold nanoparticles of $20$ nanometers mean diameter, embedded in a polymeric matrix. See [39] and references therein for experimental details and first tomographic reconstructions. Although actually measured as a dual-axis tilt series, i.e. for two perpendicular tomographic axes as sketched in figure 1(c), we consider only a single-axis sub-data set for simplicity. After background substraction and cropping to the relevant part of the field-of-view, the tomographic data consists of 75 projections of size $512\times 512$, uniformly covering an angular range of $150^{\circ}$. As our alternating minimization alignment method is not global, we prealign the data by computing the center-of-mass (COM) of every projection and shift these to center the COM. The central sinogram of the prealigned tomographic is shown in figure 10(d).

In order to limit the computational effort, we bin the data to size $256\times 256\times 75$ and reconstruct the object in a cylindrical support domain enclosed by a $256^{3}$ voxel cube. For each projection, lateral and axial shifts $s_{i,x},s_{i,y}$ as well as in-plane- and pitch-rotation angles $\theta_{i,xy},\theta_{i,yz}$ are fitted (compare figure 2), i.e. a fully 3D alignment is performed except for perturbations in the tomographic angles. The regularization parameter $\alpha$ in the Tikhonov cost functional is chosen according to an estimated misalignment of 10 pixels, see A.4.

Slices of an object recovered from the prealigned data (corresponding to the initial iterate of the alignment algorithm) are plotted in figures 10(a)-(c). Figures 10(e)-(g) show the same slices for a reconstruction based on the final alignment parameters obtained after 47 alternating reconstruction- and alignment steps of Algorithm 1. Note that the recovered objects are not the (smooth) intermediate reconstructions computed by the alignment method. Rather, they are calculated using a regularization parameter $\alpha$ that is reduced by a factor of 100 in order to facilitate the identification of misalignment artifacts and the comparison to the results in [39]. The fitted alignment parameters are plotted in figure 10(h).

Figure 10 reveals that the stripe-shaped misalignment artifacts are considerably reduced by the alternating-minimization alignment compared to the initial prealignment: figures 10(e)-(g) show well-defined compact particles as anticipated for the considered sample. However, in addition to shadow- and stripe artifacts associated with the incomplete tomographic angular range (missing wedge of $30^{\circ}$) and the coarse angular sampling, respectively, the plots show some residual crescent-shaped appendices of the particles along specific directions. This indicates a residual misalignment in the determined position of the tomographic axis.

The systematic occurence of negative values in the recovered objects suggests that the alignment might be improved by imposing a non-negativity constraint in the tomographic reconstruction. As this is not possible for the conjugate-gradients-based Tikhonov regularization, we run Algorithm 1 a second time using the alternative Kaczmarz-type tomographic reconstruction method described in A.4. The regularization parameter $\alpha$ and all other settings are chosen exactly as before.

Notably, convergence occurs much faster in this case, the stopping criterion of less than $0.05$ pixel increment in the alignment parameters being met already after 15 iterations. The final fit of alignment parameters, plotted in figure 10(l), deviates from the previous results (figure 10(h)) in the order of $\pm 1$ pixels. In order to compare the quality of the two fits, we compute a final (Tikhonov-)reconstruction with the same regularization as for the results shown in figures 10(e)-(g), yet using the alignment parameters obtained from the Kaczmacz-type approach. The recovered object slices are shown in figures 10(i)-(k). While the artifacts associated with a missing-wedge- and the angular sampling remain, the directed crescent-structures are almost completely eliminated. This indicates a significantly improved alignment of the tomographic axis.

5 Conclusions and discussion

In many applications of tomographic imaging, the experimental setup deviates geometrically from the idealized mathematical model. Such a misalignment can cause severe distortions in the reconstructed image. Therefore, an alignment procedure is needed to match the mathematical model to the true tomographic geometry. In some cases, this alignment can be done prior to reconstruction by simple transformations of the projection images. For general three-dimensional misaligment of the setup, reconstruction and alignment have to be done jointly. In this paper, we discuss the alignment problem from an algebraic point-of-view and propose an algorithmic framework for solving joint reconstruction and alignment problems. By posing the joint problem as a separable least-squares problem, we can exploit the special structure of such problems: they are quadratic in the image but generally non-convex in the alignment parameters. Since these optimization problems allow for a unique reconstructed image for a given set alignment parameters, we can eliminate the image by solving a (regularized) least-squares problem. The remaining problem is shown to be smooth in the alignment parameters (given a sufficiently differentiable interpolation scheme), even when using non-smooth regularization for the image reconstruction. The alignment problem may include additional regularization terms on the alignment parameters and can be solved using a gradient-based method. We present a general mathematical framework for designing and analyzing such gradient-based methods and discuss their convergence properties. Under the additional assumption of (local) convexity, we show that a linear convergence rate can be achieved even when performing inexact reconstructions.

We present numerical experiments on both phantom and real data to illustrate the proposed approach. For simulated full-field tomographic data, we find that small to moderate perturbations in almost all alignment parameters can be stably retrieved by Algorithm 1, fitting the alignment parameters by simple gradient-descent steps. The only exception is given by inaccuracies in the tomographic tilt angle, for which we observe a slight instability of the reconstruction to low-frequency errors. By its generic nature, our method can also be applied in the case of truncated projection data as measured in region-of-interest tomography. Indeed, our simulations show that it may correctly detect misalignment in the data even in this setting, for which already tomographic reconstruction without an alignment problem is non-unique. However, convergence turns out to be rather slow in this case. Aligning an experimental electron tomography data set by our algorithm, we demonstrate its practical applicability in more challenging, real situations, in particular involving a limited range of tomographic angles. We furthermore find that the alignment can be improved by additionally imposing a non-negativity constraint on the reconstructed images. This reveals that the efficiency of the presented alignment framework greatly benefits from the exploitation of available a priori knowledge on the imaged objects. Optimizing the incorporation of constraints is thus an important goal for further development.

It should be emphasized that, rather than developing a highly-taylored alignment method for specific types of datasets, we set out to describe a general algorithmic framework in which joint reconstruction-alignment methods can be formulated and analyzed. Our implementation contains certain features to make it generic and robust. This includes a linesearch procedure, estimation of the regularization parameter, and the ability to work with a wide range of acquisition setups. Alignment of challenging large-scale datasets, however, requires further tuning and a more efficient implementation, which is subject to future research.

Acknowledgements

This work was funded in part by Netherlands Organisation for Scientific Research (NWO) through projects 613.009.032 and 639.073.506 and by the Deutsche Forschungsgemeinschaft (DFG) through project C02 of SFB 755 - Nanoscale Photonic Imaging.

Appendix A Implementation

A.1 Discretization of the X-ray transform

We express the X-ray transform of a three-dimensional object $u$ as

[TABLE]

where $\mathbf{s}\in\mathbb{R}^{3}$ and $\boldsymbol{\eta}\in\mathbb{S}^{2}$, $\mathbb{S}^{2}=\{\mathbf{x}\in\mathbb{R}^{3}\,\,|\,\,\|\mathbf{x}\|_{2}=1\}$ denoting unit sphere, determine the origin and direction of the X-ray. The projection data consists of a collection of such projections $p_{i}=p(\mathbf{s}_{i},\boldsymbol{\eta}_{i})$, for $i=1,\ldots,m$. Sampling the object $u$ on a grid of $n$ voxels with values $u_{i}$ for $i=1,\ldots,n$ we discretize the X-ray transform as

[TABLE]

The exact form of the system matrix $W(\mathbf{a})\in\mathbb{R}^{m\times n}$ depends on the applied interpolation method, see A.3. The tomographic data can be grouped into the different projections $\mathbf{p}=(\mathbf{p}_{1},\mathbf{p}_{2},\ldots,\mathbf{p}_{N_{\textup{proj}}})^{\textup{T}}$, each acquired at a specific, yet possibly unknown set of alignment parameters $\mathbf{a}_{i}$ - a vector parametrizing the six degrees of freedom of a rigid object motion. The projection matrix then has a block-structure

[TABLE]

Importantly, this also induces a block-structure of the Jacobian to be evaluated in Algorithm 1:

[TABLE]

A.2 Parametrization of the alignment transforms

We represent the alignment parameters for one projection as $\mathbf{a}_{i}=(\mathbf{s}_{i},\boldsymbol{\theta}_{i})$, where the shift-vector $\mathbf{s}_{i}=(s_{i,x},s_{i,y},s_{i,z})\in\mathbb{R}^{3}$ contains scalar translations along the different dimensions and $\boldsymbol{\theta}_{i}=(\theta_{i,xy},\theta_{i,yz},\theta_{i,zx})\in[-\pi;\pi)$ contains the rotation angles to be applied in different coordinate planes, see figure 2.

The projection matrix $W_{i}$, that models the parallel-beam tomographic projection of the object under the rigid geometrical transform $\mathbf{a}_{i}$, can be written as

[TABLE]

where $R(\mathbf{a}_{i})$ is a matrix implementing rotations and shifts and $S_{z}$ sums the voxel values of three-dimensional arrays along the $z$-dimension. Note that only $S_{z}$ and not $R$ is specific to the considered parallel-beam geometry. Hence, it would be easy to implement alignment problems also for cone-beam setups, for example, by simply exchanging $S_{z}$ with the cone-beam projector for a fixed source position.

As shifts and rotations do not commute, we have to follow some convention on the order in which the different transforms are applied. We choose the following:

Rotation in the $z$-$x$-plane (tomographic rotation) 2. 2.

Rotation in the $y$-$z$-plane (pitch rotation) 3. 3.

Shifts in along the $x,y,z$-dimensions 4. 4.

Rotation in the $x$-$y$-plane (in-plane rotation)

As the naming suggests it, we take the $y$-axis as the tomographic axis and apply the desired rotation about this axis first before modeling additional (undesired) motions. The benefit of this ordering is that trivial components of the alignment transforms $\mathbf{a}$ can be easily identified and removed. These correspond to constant transforms of the entire object over the whole tomographic acquisition (see e.g. [3]):

• •

$\{\theta_{i,zx}\}$: The mean corresponds to a constant rotation about the tomographic axis

• •

$\{\theta_{i,xy}\},\{\theta_{i,yz}\}$: Joint sinusoidal variations $(\theta_{i,xy},\theta_{i,yz})=a(\sin(\theta_{i,zx}+b),\cos(\theta_{i,zx}+b))$ for $a>0$, $b\in\mathbb{R}$ correspond to a uniform tilt of the tomographic axis

• •

$\{s_{i,x}\}$: Sinusoidal components $s_{i,x}=c\sin(\theta_{i,zx}+d)$ correspond to a constant translation perpendicular to the tomographic axis

• •

$\{s_{i,y}\}$: The mean corresponds to a constant translation along the tomographic axis

• •

$\{s_{i,z}\}$: Shifts along the projection dimension; redundant in parallel-beam geometry

In order to ensure unambiguous results, we remove the above components from any simulated and reconstructed alignment transforms $\mathbf{a}$ by substracting the orthogonal projections onto the redundant modes.

A.3 Computation of rigid transforms and interpolation

We consider discrete objects $\mathbf{u}\in\mathbb{R}^{N_{x}\times N_{y}\times N_{z}}$, that sample the continuous density $u:\mathbb{R}^{3}\to\mathbb{R}$ on equidistant grid points $\mathbf{x}_{\textup{s}}=(x_{\textup{s}},y_{\textup{s}},z_{\textup{s}})\in\{1,\ldots,N_{x}\}\times\{1,\ldots,N_{y}\}\times\{1,\ldots,N_{z}\}=:D_{\mathbf{u}}$. The action of the geometrical transform on $\mathbf{u}$ can be written as

[TABLE]

for all $\mathbf{x}_{\textup{s}}\in D_{\mathbf{u}}$. Here, $A_{R}(\mathbf{a}_{i})^{-1}:\mathbb{R}^{3}\to\mathbb{R}^{3}$ denotes the inverse of the affine-linear map that implements the rotation and shift parametrized by $\boldsymbol{\theta}_{i},\mathbf{s}_{i}$ and $\tilde{u}$ is some interpolating function for $\mathbf{u}$ that allows its approximation on arbitrary points in $\mathbb{R}^{3}$.

We use trilinear and tricubic interpolators of the form

[TABLE]

for all $\mathbf{x}\in\mathbb{R}^{3}$. The interpolation kernels $k\in\{k_{\textup{lin}},k_{\textup{cubic}}\}$ have the properties $k(0)=1$ and $k(x)=0$ for all $x\in\mathbb{Z}\setminus\{0\}$. Hence, $\tilde{u}$ indeed interpolates $\mathbf{u}$, i.e. $\tilde{u}(\mathbf{x}_{\textup{s}})=u_{\mathbf{x}_{\textup{s}}}$ holds on all sampling points $\mathbf{x}_{\textup{s}}\in D_{\mathbf{u}}$. Furthermore, $k_{\textup{lin}}$ is Lipschitz-continuous and $k_{\textup{cubic}}$ is continuously differentiable with Lipschitz-continuous derivative. Accordingly, the same smoothness properties hold for interpolator defined in (11) and its gradient may be evaluated via

[TABLE]

Differentiability of the interpolation method for $\mathbf{u}$ allows to compute the derivative of $R(\mathbf{a}_{i})(\mathbf{u})_{\mathbf{x}_{\textup{s}}}$ with respect to the alignment parameters $\mathbf{a}_{i}$:

[TABLE]

Evaluating these derivatives for all sampling points $\mathbf{x}_{\textup{s}}\in D_{\mathbf{u}}$ then enables the computation of the sought derivatives $\frac{\partial W_{i}}{\partial\mathbf{a}_{i}}$ via (9).

We emphasize that the tricubic interpolator has just the right regularity to match the lesser differentiability conditions presented in section 2.3. According to in (11), however, this interpolation method requires evaluations of the kernel $k_{\textup{3d}}$ for the 64 nearest- and next-to-nearest neighbors of the evaluation point.

Faster evaluations of the rotation- and shift map $R(\mathbf{a}_{i})(\mathbf{u})$ can be achieved by decomposing it into three subsequent transforms along 2D-coordinate planes, i.e. we exploit that - up to interpolation errors - the following relation holds:

[TABLE]

The transforms are chosen such that each of them requires interpolation only within a specific coordinate plane of the input 3D data. In the case of $R(\mathbf{a}_{i,xy})$ for instance, all interpolation points are of the form $\mathbf{x}=(x,y,z)$ with $z\in\{1,\ldots,N_{z}\}$. (11) becomes

[TABLE]

i.e. the 3D-interpolation is reduced to 2D. Moreover, note that the interpolation weights $k(x-x_{\textup{s}})k(y-y_{\textup{s}})$ are independent of $z$. This enables efficient numerical implementation of (15) by assembling the weights into a sparse matrix and then applying it to a suitably reshaped version of $\mathbf{u}$ by a sparse-matrix-dense-matrix-product.

We refer to the implementation scheme given by (14) and (15) as bilinear (for $k=k_{\textup{lin}}$) or bicubic (for $k=k_{\textup{cubic}}$) due to the two-dimensional nature of the involved interpolations. Note that these methods result in the same smoothness of $\mathbf{a}_{i}\mapsto W_{i}(\mathbf{a}_{i})$ as trilinear and tricubic interpolation, respectively. However, interpolation errors are probably more pronounced for the former ones due to the concatenation of interpolators.

A.4 Tomographic reconstruction method

In the numerical examples presented in section 4, we make the simple choice of quadratic Tikhonov regularization with a gradient penalty term as the cost function in (2):

[TABLE]

Here, $\nabla_{\textup{dis}}$ denotes the discrete gradient operator. With this choice, the first step in Algorithms 1 – 3 amounts to solving the quadratic least-squares problem

[TABLE]

We approximate minimizers by applying the conjugate gradient method to the normal equation corresponding to (17).

As has been extensively discussed in the literature on Tikhonov regularization (see e.g. [14] and references therein), the penalty term in (17) establishes robustness of $\mathbf{u}^{(k+1)}$ against errors in the data $\mathbf{p}^{\textup{obs}}$ and thus in particular to inconsistencies due to misalignment. Accordingly, it prevents overfitting in the tomographic reconstruction step. Our motivation for penalizing the discrete gradient of $\mathbf{u}$ instead of choosing the standard penalty term $\|\mathbf{u}\|_{2}^{2}$ is two-fold: for once, misalignment often manifests in sharp stripe artifacts in the reconstruction, which can be suppressed more efficiently by smoothing penalties. On the other hand, it can be seen from (13) that the evaluation of derivatives $\frac{\partial W_{i}}{\partial\mathbf{a}_{i}}$ requires approximations of the gradient of $\mathbf{u}$. Penalizing large values of $\nabla_{\textup{dis}}\mathbf{u}$ in the tomographic reconstruction thus stabilizes subsequent alignment steps.

Minimizing the Tikhonov functional in (17) can be quite costly if the system is ill-conditioned. Moreover, the approach does not allow to impose a non-negativity constraint for the reconstructed object as is often reasonable to suppress artifacts. We therefore propose an alternative Kaczmarz-type scheme for the tomographic reconstruction in Algorithms 1 – 3:

[TABLE]

We always compute one symmetric Kaczmarz cycle, i.e. $N=2N_{\textup{proj}}$ is twice the number of tomographic projections and $m_{1}=m_{N},m_{2}=m_{N-1},\ldots m_{N/2}=m_{N/2+1}$ holds true. The processing order of the projections $\mathbf{p}^{\textup{obs}}_{m_{1}},\mathbf{p}^{\textup{obs}}_{m_{2}},\ldots$ is chosen according to the multi-level-scheme proposed in [20]. It is shown in [26] that the final Kaczmarz-iterate $\mathbf{u}^{(k)}_{N}$ minimizes of a preconditioned form of the Tikhonov functional in (17). Accordingly, the proposed scheme matches our optimization setting (2). Non-negativity of $\mathbf{u}^{(k+1)}$ can be imposed by setting negative values to zero after each Kaczmarz-iteration in (18).

For the choice of the regularization parameter $\alpha$ in (17) and (18), we propose a heuristic rule: Let $\delta a$ denote an estimate for the perturbations of the alignment parameters to be corrected, i.e. $\delta a\approx\|\mathbf{a}^{\text{(exact)}}-\mathbf{a}^{(0)}\|_{\infty}$. Then tomographic reconstruction from the misaligned data can be expected to result in artifacts roughly on the lenghtscale $\delta a$ and below. In order to suppress these, $\alpha$ is chosen such that the data fidelity- and regularization term in (17) are balanced exactly at this scale, by ensuring $\alpha\approx\|W(\mathbf{a})\mathbf{u}_{\delta a}\|_{2}^{2}/\|\nabla_{\textup{dis}}\mathbf{u}_{\delta a}\|_{2}^{2}$ for Fourier modes $\mathbf{u}_{\delta a}=\exp(\textup{i}\boldsymbol{\xi}\cdot\mathbf{x})$ of frequency $|\boldsymbol{\xi}|=\pi/\delta a$. Using the Fourier-slice-theorem, a simple formula for the choice of $\alpha$ can be obtained from this requirement. Features of size $>\delta a$ are then resolved, whereas misalignment artifacts on scales $<\delta a$ are damped out in the reconstructions. There is no guarantee that this will lead to an optimal choice for $\alpha$, especially when dealing with limited-angle or truncated data. However, the proposed heuristic has proven to be useful in practice, including the latter settings.

A.5 Alignment method

As shown in section A.1, the derivative of the projector $W=(W_{1},\ldots,W_{N_{\textup{proj}}})$ with respect to the alignment parameters $\mathbf{a}=(\mathbf{a}_{1},\ldots,\mathbf{a}_{N_{\textup{proj}}})$ exhibits a block-diagonal structure. Since we do not impose a joint regularizer for $\mathbf{a}$, the alignment steps in Algorithms 1 – 3 decouple into independent iterations for all of the $N_{\textup{proj}}$ projections:

[TABLE]

with \mathbf{s}_{i}^{(k)}:={\big{(}\nabla_{\mathbf{a}_{i}}W_{i}(\mathbf{a}_{i}^{(k)})\mathbf{u}^{(k+1)}\big{)}}^{\textup{T}}\cdot\big{(}W_{i}(\mathbf{a}_{i}^{(k)})\mathbf{u}^{(k+1)}-\mathbf{p}^{\textup{obs}}_{i}\big{)}. The step size parameters $\gamma_{i}^{(k)}$ are chosen according to exact line search for quadratic functionals:

[TABLE]

In order to avoid ill-conditioning due to different magnitudes of translational- and rotational alignment parameters, we measure the increment $s_{i}^{(k)}$ in a weighted 2-norm

[TABLE]

Here, $\odot$ denotes pointwise multiplication and the weights $w_{xy},w_{yz},w_{zx}\in\mathbb{R}_{+}$ for the rotation angles $(\theta_{i,xy},\theta_{i,yz},\theta_{i,zx})$ within the different planes, respectively, are chosen such that the rotation $\theta_{i,\ast}$ corresponds to a mean motion by $w_{xy}\cdot\theta_{i,\ast}$ object voxels.

As the cost function $f$ is not quadratic w.r.t. $\mathbf{a}_{i}$, the exact linesearch method may fail to decrease $f$, i.e. one might have an increase of the residual

[TABLE]

In this case, we sequentially reduce the proposed stepsize $\gamma_{i}^{(k)}$ by a factors of two until the residual is decreased. Note however, that this residual reduction failure never occured for the test cases presented in this work.