Simultaneous Sensing Error Recovery and Tomographic Inversion Using an Optimization-based Approach

TL;DR

This paper introduces an optimization-based method to simultaneously recover sensing errors and perform tomographic inversion, improving image accuracy despite experimental perturbations.

Contribution

It develops new physical models and optimization techniques specifically for correcting center of rotation errors in tomography.

Findings

Effective error correction demonstrated on synthetic data

Improved image reconstruction accuracy

Robustness to known rotational perturbations

Abstract

Tomography can be used to reveal internal properties of a 3D object using any penetrating wave. Advanced tomographic imaging techniques, however, are vulnerable to both systematic and random errors associated with the experimental conditions, which are often beyond the capabilities of the state-of-the-art reconstruction techniques such as regularizations. Because they can lead to reduced spatial resolution and even misinterpretation of the underlying sample structures, these errors present a fundamental obstacle to full realization of the capabilities of next-generation physical imaging. In this work, we develop efficient and explicit recovery schemes of the most common experimental error: movement of the center of rotation during the experiment. We formulate new physical models to capture the experimental setup, and we devise new mathematical optimization formulations for reliable inversion of complex samples. We demonstrate and validate the efficacy of our approach on synthetic data under known perturbations of the center of rotation.

Equations

$$Rf(\tau ,\theta )={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\delta (\tau -x\cos \theta -y\sin \theta )dxdy,$$

$$Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})=\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta\bigl{(}\tau-x_{\theta}^{*}-(x-x_{\theta}^{*})\cos\theta-(y-y_{\theta}^{*})\sin\theta\bigr{)}\>dx\>dy.$$

$$Rf(\tau ,\theta ,0,0)=Rf(\tau -P_{\theta },\theta ,x_{\theta }^{\ast },y_{\theta }^{\ast }),$$

$$P_{\theta }=x_{\theta }^{\ast }(1-\cos \theta )+y_{\theta }^{\ast }\sin \theta .$$

$$R_{\tau ,\theta }(\mathcal{W})=\sum _{v\in \mathcal{V}}{L}_v^{\tau ,\theta }{\mathcal{W}}_v,$$

$${\delta }_{P_{\theta },\sigma }(\tau )=\frac{1}{\sigma \sqrt{2\pi }}{e}^{-\frac{(\tau -P_{\theta }{)}^2}{2{\sigma }^2}},$$

$$\widetilde{Rf}(\tau ,\theta ,0,0)=Rf(\tau ,\theta ,0,0)+O({\sigma }^2)+O\left(\frac{1}{\sigma N_{\tau }}\right).$$

$$\underset{\mathcal{W}\ge 0,{\mathbf{x}}^{\ast },{\mathbf{y}}^{\ast }}{\min }\varphi (\mathcal{W},{\mathbf{x}}^{\ast },{\mathbf{y}}^{\ast })=\frac{1}{2}{\left|\left|\mathbf{L}\mathcal{W}-\mathrm{vec}\left(\mathbf{g}(\mathbf{D},{\mathbf{x}}^{\ast },{\mathbf{y}}^{\ast })\right)\right|\right|}_2^2,$$

$$\underset{\mathcal{W}\ge 0}{\min }\frac{1}{2}{\left|\left|\mathbf{L}\mathcal{W}-\mathrm{vec}\left(\mathbf{D}\right)\right|\right|}_2^2,$$

$$\begin{array}[]{rcl}\nabla\phi(\boldsymbol{\mathbf{{\mathcal{W}}}},{\bf{x}^{*}},{\bf{y}^{*}})&=&\left[\begin{array}[]{l}\displaystyle\nabla_{\boldsymbol{\mathbf{{\mathcal{W}}}}}\phi(\boldsymbol{\mathbf{{\mathcal{W}}}},{\bf{x}^{*}},{\bf{y}^{*}})\\ \displaystyle\nabla_{{\bf{x}^{*}}}\phi(\boldsymbol{\mathbf{{\mathcal{W}}}},{\bf{x}^{*}},{\bf{y}^{*}})\\ \displaystyle\nabla_{{\bf{y}^{*}}}\phi(\boldsymbol{\mathbf{{\mathcal{W}}}},{\bf{x}^{*}},{\bf{y}^{*}})\end{array}\right]\\ &=&\left[\begin{array}[]{ll}\displaystyle{\bf L}^{T}\\ \displaystyle\left(\nabla_{{\bf{x}^{*}}}g({\bf D},{\bf{x}^{*}},{\bf{y}^{*}})\right)^{T}\\ \displaystyle\left(\nabla_{{\bf{y}^{*}}}g({\bf D},{\bf{x}^{*}},{\bf{y}^{*}})\right)^{T}\\ \end{array}\right]\left({\bf L}\boldsymbol{\mathbf{{\mathcal{W}}}}-{\bf g}({\bf D},{\bf{x}^{*}},{\bf{y}^{*}})\right).\end{array}$$

$$\mathrm{vec}\left({\left[{\mathbf{D}}_{\theta ,:}\ast \left(\frac{1}{\sigma \sqrt{2\pi }}\exp \left(\frac{-(\mathcal{T}-P_{\theta }{)}^2}{2{\sigma }^2}\right)\right)\circ \frac{\mathcal{T}-P_{\theta }}{{\sigma }^2}\left(\cos \theta -1\right)\right]}_{\theta =1,\dots ,N_{\theta }}\right),$$

$$\mathrm{vec}\left({\left[{\mathbf{D}}_{\theta ,:}\ast \left(\frac{1}{\sigma \sqrt{2\pi }}\exp \left(\frac{-(\mathcal{T}-P_{\theta }{)}^2}{2{\sigma }^2}\right)\right)\circ \frac{\mathcal{T}-P_{\theta }}{{\sigma }^2}\sin \theta \right]}_{\theta =1,\dots ,N_{\theta }}\right),$$

$$\underset{\mathcal{W}\ge 0,\mathbf{P}}{\min }\varphi (\mathcal{W},\mathbf{P})=\frac{1}{2}{\left|\left|\mathbf{L}\mathcal{W}-\mathbf{g}(\mathbf{D},\mathbf{P})\right|\right|}_2^2,$$

$$\text{SSIM}(a,b)=\frac{(2{\mu }_a{\mu }_b+c_1)(2{\sigma }_{ab}+c_2)}{({\mu }_a^2+{\mu }_b^2+c_1)({\sigma }_a^2+{\sigma }_b^2+c_2)},$$

$$\underset{\mathcal{W}\ge 0}{\min }\frac{1}{2}{\left|\left|\mathbf{L}\mathcal{W}-\mathrm{vec}\left(\mathbf{D}\right)\right|\right|}_2^2+\lambda {\left|\left|\mathcal{L}\mathcal{W}\right|\right|}_2^2$$

$$\underset{\mathcal{W}\ge 0}{\min }\frac{1}{2}{\left|\left|\mathbf{L}\mathcal{W}-\mathrm{vec}\left(\mathbf{D}\right)\right|\right|}_2^2+\lambda {\left|\left|\nabla \mathcal{W}\right|\right|}_1,$$

$$Rf(\tau ,\theta )=Rg({\tau }^{\prime },\theta ),$$

$$\begin{array}[]{rl}Rf(\tau,\theta)&=\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(x,y)\delta(\tau-x\cos\theta-y\sin\theta)dxdy\\ &=\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}g(x-\Delta x,y-\Delta y)\delta(\tau-x\cos\theta-y\sin\theta)dxdy\\ &=\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}g(x^{\prime},y^{\prime})\delta(\tau-(x^{\prime}+\Delta x)\cos\theta-(y^{\prime}+\Delta y)\sin\theta)dx^{\prime}dy^{\prime}\\ &=\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}g(x^{\prime},y^{\prime})\delta(\tau^{\prime}-x^{\prime}\cos\theta-y^{\prime}\sin\theta)dx^{\prime}dy^{\prime}\\ \vspace{0.1cm}\hfil&=\displaystyle Rg(\tau^{\prime},\theta).\end{array}$$

$$Rf(\tau ,\theta )=Rg(\tau ,\theta -\Delta \varphi ).$$

$$\begin{array}[]{rl}Rf(\tau,\theta)&=\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(r,\phi)\delta(\tau-r\cos\phi\cos\theta-r\sin\phi\sin\theta)|r|drd\phi\\ &=\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}f(r,\phi)\delta(\tau-r\cos(\phi-\theta))|r|drd\phi\\ &=\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}g(r,\phi-\Delta\phi)\delta(\tau-r\cos(\phi-\theta))|r|drd\phi\\ &=\displaystyle\int_{-\infty}^{\infty}\int_{-\infty}^{\infty}g(r,\phi^{\prime})\delta(\tau-r\cos(\phi^{\prime}+\Delta\phi-\theta))|r|drd\phi^{\prime}\\ &=\displaystyle Rg(\tau,\theta-\Delta\phi).\end{array}$$

$$\ell (x_{\theta }^{\ast },y_{\theta }^{\ast })=\left\{(\tau -x_{\theta }^{\ast })\left[\begin{array}{c}\cos \theta \\ \sin \theta \end{array}\right]+\left[\begin{array}{c}{x}_{\theta }^{\ast }\\ y_{\theta }^{\ast }\end{array}\right]+s\left[\begin{array}{c}-\sin \theta \\ \cos \theta \end{array}\right]:s\in \mathbb{R}\right\}.$$

$${\tau }^{\prime }=\tau -x_{\theta }^{\ast }+x_{\theta }^{\ast }\cos \theta +y_{\theta }^{\ast }\sin \theta$$

$$Rf(\tau ,\theta ,x_{\theta }^{\ast },y_{\theta }^{\ast })={\int }_{\ell (x_{\theta }^{\ast },y_{\theta }^{\ast })}f|ds|={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }f(x,y)\delta ({\tau }^{\prime }-x\cos \theta -y\sin \theta )dxdy,$$

$$f(0)\approx {\int }_{-\infty }^{\infty }f(x){\delta }_{\sigma }(x)dx,$$

$${\delta }_{\sigma }^{\mathrm{T}}(x;A)=\{\begin{array}{ll}{\delta }_{\sigma }(x)& |x|\le A\\ 0& |x|>A,\end{array}$$

$$E(\sigma ,A,K)=f(0)-\frac{2M}{K}\sum _{k=-N}^{N}f(x_k){\delta }_{\sigma }^{\text{T}}(x_k;A),$$

$$E_{\text{R}}(\sigma )=f(0)-{\int }_{-\infty }^{\infty }f(x){\delta }_{\sigma }(x)dx.$$

$$E_{\text{T}}(\sigma ,A)={\int }_{-\infty }^{\infty }f(x){\delta }_{\sigma }(x)dx-{\int }_{-\infty }^{\infty }f(x){\delta }_{\sigma }^{\text{T}}(x;A)dx.$$

$$E_{\text{D}}(\sigma ,A,K)={\int }_{-\infty }^{\infty }f(x){\delta }_{\sigma }^{\text{T}}(x;A)dx-\frac{2M}{K}\sum _{k=-N}^{N}f(x_k){\delta }_{\sigma }^{\text{T}}(x_k;A).$$

$$E(\sigma ,A,K)=E_{\text{R}}(\sigma )+E_{\text{T}}(\sigma ,A)+E_{\text{D}}(\sigma ,A,K).$$

Full text

{AMS}

68Q25, 68R10, 68U05

1 Introduction

Tomographic imaging has had a revolutionary impact on medicine, physics, and chemistry. Even so, the problem of reconstructing an image from tomographic data remains challenging in many interesting cases, such as when the amount of available data is limited and/or the problem is ill posed (in the sense that canonical metrics used to assess the discrepancy between a reconstruction and the measured data generally possess many local minima). This ill posedness makes reconstructions susceptible to experimental errors, in particular, to errors stemming from mismatches between the experimental configuration and the assumptions of the measurement process. Recovering such errors is crucial for realizing the gains from improvements in measurement and experimental hardware, such as the improved resolution promised by brighter, more coherent next-generation light sources. Sample drift [28] and beam drift [36] are two fundamental sources of error. These errors can often arise from the drift of the center of rotation (CoR) of the imaging stage during data acquisition [2], which is the focus of this work.

In computerized tomography, the object or sample being imaged is placed on a stage and irradiated with parallel beams of x-rays. As the rays pass through the object, they are partially absorbed according to the object’s composition. Radiation that passes through the object unabsorbed is collected by a detector, producing an x-ray “shadow” or “projection” of the object. By rotating the object and repeating the process, one obtains projections of the object as illuminated from several different directions. One then uses these projections to reconstruct an image of the object and its interior [22].

Imperfections in the experimental apparatus can cause the CoR of the imaging stage to vary slightly over the course of this process. When this happens, the projections from different angles will not be properly aligned relative to one another, and failing to account for this misalignment during reconstruction may result in a smeared image (loss of resolution) or, worse, an image that is completely incorrect. Moreover, even a small error in the CoR can yield a large error in the reconstruction, as is illustrated by the experiment of Fig. 1. The image on the left shows the object—a small, solid circle—and the image in the middle shows the result of a standard tomographic reconstruction without addressing the CoR shift from a sequence of measurements in which the CoR is perturbed by a 2% displacement relative to the size of the object domain. The reconstructed image obtained by not correcting for this perturbation looks very different from the true object. We also demonstrate the reconstruction of the approach we propose in this paper on the right side, which almost perfectly resembles the true object.

In the past, it has been possible largely to ignore errors like these, as the drifts have been small compared with the widths of the beams used to illuminate the object. With demand for increasingly finer resolutions leading from higher-quality light sources with narrower beams, however, these errors can no longer be ignored. Potential approaches to this problem include the incorporation of prior knowledge into the reconstruction process using Bayesian frameworks [17, 33, 34] and the use of regularizers to promote sparsity or smoothness of the reconstruction [14, 10, 23, 30, 20]. One can also collect additional data (beyond x-ray projections) that is less sensitive to experimental errors and use this data to help with the reconstruction [7]. Unfortunately, the improvement in tomographic hardware has been so great that these generic approaches, which apply more widely to problems other than CoR drift recovery, are not likely to be sufficient. In the case of CoR drift, we show (in Sec. 5) that standard regularization approaches fail even when the regularization parameter value is chosen in an ideal way.

Much research effort has been devoted to tackling the CoR recovery problem. For example, a typical approach to recovering a single, unknown CoR is to use a pair of projections that are reflections of one another to estimate a detector offset that can be used to shift the projections into the correct positions [40, 27]. This approach is highly sensitive to the accuracy of the mirrored alignment, and its use of such a small number of projections to perform the correction makes it susceptible to noise and other effects of limited data, potentially resulting in a low-contrast recovery. Moreover, mirrored pairs of projections are not always available in practice. Azvedo et al. [2] proposed a method to estimate the CoR based on the preservation of photon counts passing across the sample, which requires almost perfect measurement. The most common technique is to compute cross correlations between projections acquired from successive rotations [11, 1, 18]; however, this approach is limited to simple and relatively homogeneous samples, since there is no good way to rank one feature as “more important” than another when two projections contain very distinguishable features. A similar approach involves manual alignment of the projections using a known “hotspot” in the object as a reference. This process can be labor-intensive and cannot be used if a hotspot cannot be identified. Recently, a new technique known as iterative reprojection has been introduced in which one recovers projection alignments and reconstructs the image simultaneously using an iterative procedure [35, 19, 39, 13]. The basic idea is to alternate between a few iterations of projection alignment and a few iterations of reconstruction until a “forward model” of the experimental setup and an “inverse model” of the reconstruction process are consistent with one another. This approach can yield better accuracy, especially for noisy, limited data.

One deficiency of the iterative reprojection methods that have emerged thus far is that the update of the projection alignment and the update of the reconstruction are not fully coupled. As a result, some mismatch remains between the alignment of the projections and the experimental configuration assumed by the reconstruction process. In this work, we propose to address this problem using a novel joint inversion framework based on optimization in which we explicitly model and recover for CoR drift. Our proposed model is flexible, making no assumptions about when and where the drift happens. Our new approach is easier to automate than existing approaches and can be more robust in the context of poor data quality and limited prior knowledge of the object being imaged.

In Sec. 2, we describe our mathematical forward models of the CoR drift error in the experiment and show how to embed these models into the reconstruction scheme. In Sec. 3, we describe our simultaneous reconstruction approach for recovering the object and the experimental error, including the formulation of the objective function. In Sec. 4, we describe the algorithm for solving the resulting optimization problem and its complexity. In Sec. 5, we present some numerical illustrations comparing the performance of our simultaneous inversion method with that of existing approaches using a few synthetic examples. In Sec. 6, we summarize the proposed method and discuss a few directions for future research.

For simplicity, we confine ourselves to reconstructing 2D images, although our methods work just as well in 3D.

2 Mathematical Model

In this section, we describe our model for the tomographic imaging process. For further details on tomography, we refer the reader to [21].

2.1 Radon Transform

The fundamental mathematical tool in tomography is the Radon transform [32], defined for a compactly supported function $f:\mathbb{R}^{2}\mapsto\mathbb{R}$ by

[TABLE]

where $\delta$ is the Dirac delta function. Throughout this article, we consider the restricted domain given by $\tau\in[0,\infty)$ and $\theta\in[0,2\pi)$, and we assume, as will always be the case in practice, that $f$ is well-enough behaved that the integral makes sense. The Radon transform of a function is frequently called its sinogram.

The projections that are measured in the tomographic x-ray imaging process are values of the Radon transform of the object’s attenuation coefficient $f(x,y)$, a function that describes the propensity of the object to absorb x-rays at each point in the object’s support. The angle $\theta$ encodes the direction from which the beams approach the object; and given that direction, $\tau$ encodes the position of the beam. If the CoR shifts from the origin to $(x_{\theta}^{*},y_{\theta}^{*})\in\mathbb{R}^{2}$, as the imaging stage rotates through the angle $\theta$, the projections we measure come not from (2.1) but from the CoR-shifted Radon transform

[TABLE]

For a derivation of this equation, see Appendix B. As illustrated in Fig. 2, $Rf(\tau,\theta,0,0)$ is the integral of $f$ along the purple line perpendicular to the direction determined by the angle $\theta$ and at a distance $\tau$ from the origin. Alternatively, the purple line is obtained by rotating from its previous position, which overlays with the $y$-axis, $\theta$ degrees with respect to CoR $(0,0)$. However, if the CoR is shifted from $(0,0)$, denoted by the purple dot, to the green dot, we obtain a different projection as $Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})$.

To recover for the difference in CoR, we need a way to convert (2.2) back into (2.1). This is easy: we have

[TABLE]

where

[TABLE]

The simple relationship (2.3)–(2.4) forms the basis for our algorithm. We rotate the imaging stage through an angle $\theta$ with the intent of measuring $Rf(\tau,\theta,0,0)$, but we instead measure $Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})$ due to drift of the CoR. The identity (2.3) shows that all we need to do is translate (in $\tau$) the Radon transform that we measured by an amount $P_{\theta}$, which is related to the true CoR by (2.4). As we later show, given sufficient data, we can estimate the drift-induced $P_{\theta}$ using an optimization procedure.

2.2 Discrete Tomography

In practice, we cannot recover the desired object property (e.g., attenuation coefficient) at all points in space. Instead, we discretize the space (containing the compact object) being imaged into $N\times N$ pixels and vectorize it to an array ${\mathcal{V}}$. Let ${\mathcal{W}}_{v}$ be the value of the object property we intend to recover on the pixel $v\in{\mathcal{V}}$, and let $\boldsymbol{\mathbf{{\mathcal{W}}}}=\{{\mathcal{W}}_{v}\>:\>v\in{\mathcal{V}}\}$ denote the discretized image. Given parameters $\tau$ and $\theta$, we calculate a discrete Radon transform of $\boldsymbol{\mathbf{{\mathcal{W}}}}$ via

[TABLE]

where $L_{v}^{\tau,\theta}$ is the length of the intersection of the beam described by $\tau$ and $\theta$ with the pixel $v$, see Fig. 2. It is this relationship that we must invert to reconstruct our image: knowing $R_{\tau,\theta}(\boldsymbol{\mathbf{{\mathcal{W}}}})$, we wish to find ${\mathcal{W}}_{v}$ for each pixel $v$.

In addition to needing to consider the image discretely, we must also “discretize the beams”: we have access only to $R_{\tau,\theta}$ for a limited number of values of $\tau$ and $\theta$. The values of $\tau$ we use are fixed and equally spaced and are identical for each $\theta$. We denote the collection of values of $\tau$ (sometimes called “beamlets”) by ${\cal T}$ and the collection of values of $\theta$ by $\Theta$; accordingly, $N_{\tau}=|{\cal T}|$ is the number of beamlets, and $N_{\theta}=|\Theta|$ is the number of angles. In reality, the resolution of ${\cal T}$ is decided by the energy level of the radiation source and by the detector resolution [3].

Thus, for each angle $\theta\in\Theta$, our measurement apparatus ideally produces a set of samples of $Rf(\tau,\theta)$, equally spaced in $\tau$, which we take as values of $R_{\tau,\theta}(\boldsymbol{\mathbf{{\mathcal{W}}}})$. Because of the drift in the CoR, however, our samples are actually from $Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})$ instead. We can account for this drift by recomputing the $L_{v}^{\tau,\theta}$ to be consistent with the change in CoR, but we wish to avoid doing so because calculating these values (potentially many times) is expensive. Instead, we calculate them once, assuming that the CoR is at the origin, and recover for the change in CoR by translating the projections.

2.3 Aligning the Projections

In principle, translating $Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})$ in the $\tau$ variable by $P_{\theta}$ is an easy task—the answer is just $Rf(\tau-P_{\theta},\theta,x_{\theta}^{*},y_{\theta}^{*})$—but since we have access only to samples of $Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})$ at $\tau\in{\cal T}$, we cannot do this. Since the samples come from equally spaced points, a natural idea is to effect the translation by using the discrete Fourier transform, which is mathematically equivalent to forming a trigonometric interpolant to $Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})$ through the points $\tau\in{\cal T}$ and translating the interpolant. This approach works, but the images recovered in this way in practice are contaminated with ringing artifacts [4, pp. 209]. To fix this problem, we apply a low-pass filter—a normalized Gaussian filter with standard deviation $\sigma$—to damp the high-order Fourier coefficients.

An alternative way to understand the translation process is as follows. Translating $Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})$, viewed as a function of $\tau$ only, by $P_{\theta}$ is equivalent to convolving $Rf$ with a Dirac delta function centered at $P_{\theta}$: $Rf(\tau-P_{\theta})=(Rf\ast\delta_{P_{\theta}})(x)$, where $\delta_{P_{\theta}}(\tau)=\delta(\tau-P_{\theta})$. We regularize $\delta_{P_{\theta}}$ by replacing it with a Gaussian,

[TABLE]

and instead compute the convolution $(Rf\ast\delta_{P_{\theta},\sigma})(\tau)$, which we do by sampling $Rf$ and $\delta_{P_{\theta},\sigma}$ on the same equally spaced grid and using the discrete Fourier transform.

How should one choose the hyperparameter $\sigma>0$? Let $\widetilde{Rf}(\tau,\theta,0,0)$ be our approximation to $(Rf\ast\delta_{P_{\theta},\sigma})(\tau)$111As in the preceding paragraph, the $Rf$ in $Rf\ast\delta_{P_{\theta},\sigma}$ in this definition is understood to be the CoR-shifted Radon transform. Our notation $\widetilde{Rf}(\tau,\theta,0,0)$ emphasizes our hope that $\widetilde{Rf}(\tau,\theta,0,0)\approx Rf(\tau,\theta,0,0)$. obtained from the discrete Fourier transform. By the error analysis presented in Appendix C, if the (attenuation coefficient of the) object being imaged is twice continuously differentiable, then

[TABLE]

The second term on the right-hand side of (2.6) represents the error incurred in our Gaussian regularization of the delta function, while the third term represents the discretization error due to the limited number of measurements. A trade-off between the two exists: a narrower Gaussian (smaller $\sigma$) implies a more faithful translation, but it requires more measurements (greater $N_{\tau}$) to approximate accurately.

We choose to have the full width at half maximum (FWHM) of the Gaussian regularizer to cover one unit of beamlet width. Since the FWHM for a Gaussian with variance $\sigma^{2}$ is given by ${\mathrm{FWHM}}\approx 2.355\sigma$, this means that we choose $\sigma=1/2.355\approx 0.42$. Observe that this choice is supported by Fig. 3, which displays the approximation error in Eqn. (2.6) in the simulated beam data from Fig. 4. For the values of $N_{\tau}$ considered, the choice ${\mathrm{FWHM}}=1$ is roughly the point at which the second error term in (2.6) takes over from the first.

2.4 Numerical Illustration

To illustrate the effects of the approximation (2.6) in the context of a complete image, we perform a simulation using a standard MRI test image. We simulate beam data from a variety of angles, assuming that as the object rotates, it does so about a CoR $(x^{*}_{\theta_{i}},y^{*}_{\theta_{i}})=(x^{*}_{\theta_{j}},y^{*}_{\theta_{j}})\neq(0,0)$ where $\theta_{i}\in 1,\dots,N_{\theta}$ and $\theta_{j}\in 1,\dots,N_{\theta}$ represent two different angles; in other words, there is only one CoR drift in the beginning of the experiment. The results are displayed in Fig. 4a, where $N=128$, $N_{\tau}=\lfloor\sqrt{2}N\rfloor$, and $N_{\theta}=30$ equally spaced in $[0,2\pi)$. The leftmost panel shows the test image; the red dot marks the CoR $(x_{\theta}^{*},y_{\theta}^{*})$. The second panel shows the sinogram obtained from this data. The third panel shows the sinogram that we would have obtained in the absence of the drift in CoR. The fourth panel shows the sinogram that we obtain after performing the recovering translations using the approach outlined in the preceding subsection. The drift-free and drift-recovered sinograms are in excellent agreement, even in the presence of the errors described by (2.6).

In Fig. 4b, we repeat the experiment of Fig. 4a except that the CoR changes each time the object is rotated instead of just once at the beginning of the experiment; in other words, $(x^{*}_{\theta_{i}},y^{*}_{\theta_{i}})\neq(x^{*}_{\theta_{j}},y^{*}_{\theta_{j}})\neq(0,0)$. The trail of red dots in the leftmost panel shows how the CoR drifts across the image as the data is acquired. The drift-free and drift-recovered sinograms are again in excellent agreement. From a practical perspective, the additional “jitter” due to the multiple shifts of CoR—readily apparent in the measured sinogram—complicates the alignment process significantly. Nevertheless, the algorithm we propose next is able to handle multiple shifts.

3 Optimization-Based Reconstruction Algorithm

In the experiments of Sec. 2.4, we assumed exact knowledge of the CoR $(x_{\theta}^{*},y_{\theta}^{*})$ and the corresponding translation parameter $P_{\theta}$. In practice, we must estimate these parameters from the data. As described in Sec. 1, the best techniques used to date generally do this iteratively by alternating between using the current parameters $P_{\theta}$ to reconstruct the image and using the reconstruction to update the parameters. Instead, we propose to estimate the parameters and reconstruct the image concurrently.

The most straightforward approach would be to try to estimate $x_{\theta}^{*}$ and $y_{\theta}^{*}$ and obtain $P_{\theta}$ from these estimates using (2.4); we refer to this as the explicit approach, since we try to find the CoR explicitly. Alternatively, we can skip the estimation of $x_{\theta}^{*}$ and $y_{\theta}^{*}$ and try to solve for $P_{\theta}$ directly; we call this the implicit approach. One major difference between the explicit and implicit approaches is the different dimension of unknown variables. In general, each $P_{\theta}$ corresponds to a pair of $(x_{\theta}^{*},y_{\theta}^{*})$, which results in $N_{\theta}$ number of CoRs to be recovered. If one knows which angles have the CoR drifts, however, the number of $(x_{\theta}^{*},y_{\theta}^{*})$ pairs can be reduced to the exact number of drifts.

To exploit the correlations between $(x_{\theta}^{*},y_{\theta}^{*})$ and $\boldsymbol{\mathbf{{\mathcal{W}}}}$, we formulate the reconstruction problem as simultaneously recovering the CoR and recovering the object. Therefore, the final explicit reconstruction problem is

[TABLE]

where ${\bf{x}^{*}}=[x^{*}_{\theta}]_{\theta=1}^{N_{\theta}}$, ${\bf{y}^{*}}=[y^{*}_{\theta}]_{\theta=1}^{N_{\theta}}$, $\boldsymbol{\mathbf{{\mathcal{W}}}}\geq 0$ is due to the physical nature of mass, ${\bf D}\in R^{N_{\theta}\times N_{\tau}}$ is the measurement data, ${\bf g}({\bf D},{\bf{x}^{*}},{\bf{y}^{*}})=\left[D_{\theta,\tau}*\delta_{P_{\theta}(x_{\theta}^{*},y_{\theta}^{*}),\sigma}\right]_{\theta,\tau}\in\mathbb{R}^{N_{\theta}\times N_{\tau}}$ is the translated sinogram by (2.4), and ${\bf L}=[L_{v}^{\tau,\theta}]\in\mathbb{R}^{N_{\theta}N_{\tau}\times N^{2}}$ is determined based on the standard 2D Radon mapping with an implied CoR of $(0,0)$.

As a side note, tomographic reconstruction without CoR error recovery is typically formulated as

[TABLE]

which is equivalent to Eq. (3.7) for $({\bf{x}^{*}},{\bf{y}^{*}})=0.$

The first-order derivative of the objective function (3.7) is

[TABLE]

In (3.9), $\nabla_{{\bf{x}^{*}}}g({\bf D},{\bf{x}^{*}},{\bf{y}^{*}})$ and $\nabla_{{\bf{y}^{*}}}g({\bf D},{\bf{x}^{*}},{\bf{y}^{*}})$ are

[TABLE]

and

[TABLE]

respectively, where $\circ$ is the Hadamard product, and ${\bf D}_{\theta,:}$ is the $\theta\,th$ row of ${\bf D}$. A major difference compared with the single CoR-shift case is the increased computational cost. Instead of having only two extra parameters, this case has $2N_{\theta}$ extra parameters.

For now, we have been focusing on explicitly recovering the coordinates of the CoRs $(x_{\theta}^{*},y_{\theta}^{*})$. Alternatively, we can reformulate the optimization problem as finding the optimal shifts $P_{\theta}$ for each angle, which results in the implicit problem

[TABLE]

where ${\bf P}=[P_{\theta}]_{\theta=1,\dots,N_{\theta}}$ and ${\bf g}({\bf D},{\bf P})=\left[D_{\theta,\tau}*\delta_{P_{\theta},\sigma}\right]_{\theta,\tau}\in\mathbb{R}^{N_{\theta}\times N_{\tau}}$. The derivative of the objective function in (3.10) is similar to (3.9) but simpler, so we will not describe it here. In the case of multiple CoR shifts, the advantage of this formulation is the reduced number of extra parameters by a factor of 2. Notice that given a fixed CoR, the Radon transform has the property that different transformations (e.g., translation or rotation) of an object will result in different sinograms; see Appendix A. These properties suggest that any pair of transformed object and its corresponding sinogram is an optimal solution of the reconstruction. Therefore, the proposed optimization problem, either explicit or implicit, will not have a unique solution. In other words, if two objects are related by equal (up to translations and rotations), then the two objects, together with their corresponding sinograms, are equivalent from the perspective of the optimization problems in (3.7) and (3.10).

4 Optimization Complexity and Computational Expense

Since the complexity of Eq. (2.4) is negligible, we analyze only the computational complexity of the implicit approach (3.10). The calculation of Eq. (3.10) requires about $N^{2}N_{\theta}$ flops given the relationship $N_{\tau}=\lfloor\sqrt{2}N\rfloor$, and it includes 1 misfit calculation and 3 FFTs required in operator ${\bf g}$. Figure 5 shows a log-log plot of the computational time of one (function, gradient) evaluation for increased number of $N$, given $N_{\theta}=1$. With a model fit to the time result, the time complexity is on the order of $N^{2}$ which is consistent with our analytical approximation considering $N_{\theta}\ll N$ in general.

The optimization solver we use in the numerical experiments is an inexact truncated-Newton (TN) method [29] with a preconditioned conjugate gradient (PPCG) [31, pp. 460] subproblem solver to compute the search direction. The result, summarized in Alg. 1, is a large-scale optimization solver with memory efficiency well suited for the high computational complexity required by the resulting nonlinear, nonconvex optimization problem.

In our implementation, we do not form the Hessian $\nabla^{2}\phi$; instead, we approximate the Hessian-vector product $(\nabla^{2}\phi)^{T}d$ required in PPCG by taking finite differences with $\nabla\phi$ values. An estimation of the complexity of TN is provided as follows. Each outer iteration needs the following computations with respect to the number of unknown parameters:

• •

1 infinity-norm calculation, 1 vector addition, and 2 (function, gradient) evaluations

• •

a number of PPCG iterations, with cost per inner iteration given by

• –

1 (function, gradient) evaluation,

• –

4 inner products, and

• –

5 “vector+constant$\cdot$vector” operations

In our experience, on average, 5 PPCG iterations are required per outer TN iteration. The overall cost of one TN iteration for solving (3.10) amounts to $(1+5+20+25+1)(|{\mathcal{V}}|+N_{\theta})=52(|{\mathcal{V}}|+N_{\theta})$ floating-point operations, plus 7 (function, gradient) evaluations, whose complexity is $O(N^{2})$.

5 Numerical Results

In this section, we examine the performance of the algorithm for both the explicit and implicit cases. The primary goal of our tests is to measure the performance of how well we recover the CoR shift with respect to different initializations, different objects, and different levels of noise in the data. All numerical experiments are performed on a platform with 32 GB of RAM and two Intel E5430 Xeon CPUs. Throughout all the numerical tests, we fix the experimental setup for all the tests as $N=128$, $N_{\tau}=\lfloor\sqrt{2}N\rfloor$, and $N_{\theta}=30$; that is, the object size is $128\times 128$ pixels, and the tomographic data is measured by collecting 30 projections over a full rotation of $2\pi$ radians. We also choose the stopping criteria of TN to be $\left|\left|\nabla\phi\right|\right|<10^{-5}$.

As indicated in Sec. 3, our proposed optimization problem does not have a unique solution in the sense that the reconstructed object can be a translated or rotated version of the ground truth. Therefore, the error metric we use to measure the reconstruction quality cannot be simply the mean squared error. To resemble the human visual system, we utilize the structural similarity (SSIM) metric from [37] to quantify the reconstruction quality. Given two images $a$ and $b$ of the same dimension, the SSIM index is a measure of the similarity between $a$ and $b$ and is defined by

[TABLE]

where $c_{1}$ and $c_{2}$ are small, positive constants; $\mu_{a}$ and $\sigma_{a}^{2}$ are the mean and variance, respectively, of $a$; $\mu_{b}$ and $\sigma_{b}^{2}$ are the mean and variance, respectively, of $b$; and $\sigma_{ab}$ is the covariance of $a$ and $b$. Notice that for nonnegative $a$ and $b$, $\text{SSIM}(a,b)\geq 0$. The closer the value of $\text{SSIM}(a,b)$ is to $1$, the greater is the similarity between $a$ and $b$.

We emphasize that in this work our focus is on the development of a joint framework to account for general CoR shifts. To isolate the contributions of our proposed approach to the reconstruction performance, we do not consider any additional regularization technique in this work. Additionally, the nonnegative constraint on $\boldsymbol{\mathbf{{\mathcal{W}}}}$ serves as a soft regularizer, which means we do not need to add a regularizer and can instead show the benefit of our approach by itself [8, 5]. As suggested by Fig. 6, where we test the performance of TN on various settings for the traditional tomography problem (3.8), we observe that as long as the dimension of Null(${\bf L}$) is close to the number of zeros of the object, the reconstruction quality is stable and satisfactory based on SSIM and the objective value of problem (3.8). Since the object domain always has a zero support to account for the field of view during the rotation, the dimension of Null(${\bf L}$) typically is proportional to the number of zeros in $\boldsymbol{\mathbf{{\mathcal{W}}}}^{*}$, and this explains the satisfactory reconstruction result from seemingly underdetermined system.

5.1 Analysis of Initialization

In this section, we test the sensitivity of the proposed algorithm to the initialization of $\boldsymbol{\mathbf{{\mathcal{W}}}}$.

We first consider the case of a single CoR shift that happens at the beginning of the experiment on a standard Shepp-Logan phantom object. The CoR is denoted by the red dot on top of the ground truth in Fig. 7 (far left). The upper left panel of Fig. 7 shows the standard reconstruction without considering the shift of the CoR, and we observe that it is far from the ground truth. We also include reconstructions obtained by adding standard regularizers. In particular, we consider $L2$- and $TV$-regularizers [25, 15] on the variable $\boldsymbol{\mathbf{{\mathcal{W}}}}$, respectively given by

[TABLE]

and

[TABLE]

where $\mathcal{L}\in\mathbb{R}^{N^{2}\times N^{2}}$ is the Laplacian operator, $\nabla$ denotes the spatial gradient, and $\lambda>0$ is the regularizer parameter that balances the misfit term and regularization term. One approach to choose $\lambda$ would be the L-curve method [16]. Instead, however, we choose $\lambda\in[10^{-10},10^{2}]$ to be the value that gives the best reconstruction as measured by the SSIM index (see Fig. 8). This is an idealistic choice of $\lambda$ that is impractical. The results given in the upper middle and upper right panels demonstrate this assertion and show that standard regularization techniques are not able to improve the reconstruction quality in the presence of even a single CoR shift. Since even an optimal $\lambda$ performs poorly for CoR recovery, for the rest of our experiments we do not report results from regularization techniques.

We also compare our reconstruction result with the most popular approach for the case of a single CoR drift [27]. The basic idea is to utilize the fact that when the CoR is the origin, the projection at angle $0^{\circ}$ should be a reflection of the projection at angle $180^{\circ}$. Therefore, the CoR shift can be estimated by aligning these two projections. In our case, as is often the case in practice, exactly reflecting directions are not available; instead, we choose the pair of angles $1^{\circ}$ and $175^{\circ}$ from the available 30 projections. We align these two projections by cross-correlation [40] and use the estimated CoR shift to translate the collected sinogram. The reconstruction from the shifted sinogram is reported in the lower left panel. We also report the reconstruction using the true CoR in the lower right panel. The lower middle panel shows the reconstruction obtained by the proposed explicit approach, which clearly resembles the ground truth the best.

The next test is to optimize the explicit problem (3.7) and the implicit problem (3.10) in order to recover $(x_{\theta}^{*},y_{\theta}^{*})$ and ${\bf P}$, respectively. Given an initial $(x_{\theta}^{*},y_{\theta}^{*})$ as the center, left top corner (NW), left bottom corner (SW), right top corner (NE), and right bottom corner (SE) of the object domain, respectively, and initialization $\boldsymbol{\mathbf{{\mathcal{W}}}}^{*}=0$, the performance is shown in Fig. 9. In this test, the explicit approach consistently recovers objects that are better compared with the implicit approach, which struggles for certain initializations. This observation is also reflected in the objective value and the reconstruction quality. In addition, the explicit approach takes fewer iterations to converge to smaller objective values compared with the implicit approach. The explicit approach also provides relatively better reconstruction quality as measured by SSIM index.

In Fig. 10, we also show the iterative progress of the optimization problems (3.8), (3.7), and (3.10) respectively for the single CoR-shift case. The left graph shows the progress of the objective value for the three approaches. Standard tomographic reconstruction, which refers to reconstruction without error recovery, converges much more slowly than the other two with a much higher objective value, while the explicit approach reduces the objective value even further than the implicit approach because of the reduced number of variables. The right graph compares the progress of the SSIM index for the three approaches, which is somewhat consistent with the function value reduction. We observe that the explicit approach provides the best reconstruction quality in terms of SSIM compared with the other two approaches, while the standard approach performs the worst because it does not attempt to recover for the error in CoR.

The next test concerns solving the reconstruction problem with multiple CoR shifts in the sense that every rotation has its own CoR. This problem setting is illustrated in Fig. 12, where the ground truth is again the Shepp-Logan phantom and the trajectory of the dynamic CoR is denoted by connected red dots. The second panel shows the reconstruction without CoR recovery, and the third panel shows the reconstruction with error recovery. Again, we tested the performance of the explicit and the implicit algorithm; the results are shown in Fig. 13. This time we observe that the reconstruction quality of the implicit approach is comparable to the explicit approach but is more sensitive to the initialization. In the multiple CoR-shifts case, however, the convergence of the implicit approach is faster because of its reduced number of variables compared with explicit approach. Again, the performance difference between the explicit and implicit approaches is reflected in the objective value and the corresponding SSIM index. More importantly, as shown in Fig. 11, both approaches are able to remove the jitters introduced by multiple CoR shifts and return a clean sinogram. The consistency between the objective value $\phi$ and SSIM index suggests that, in the case of limited knowledge about the object, the objective value can be reliable to judge the reconstruction quality of different local minima.

We note that as long as the shift parameters ${\bf P}$ are initialized to zero, the performance is relatively stable and less susceptible to the initialization of the object. The reason is that choosing the initial shift as zeros is equivalent to providing the initialization of the object as the reconstruction without error recovery, which is closer to the true solution than the reconstruction from a randomly perturbed sinogram. Therefore, choosing zeros for the initial shifting parameters performs better, and this initialization will be our choice for the rest of the tests.

Similar to the single CoR case, we compare the iterative progress provided by the standard reconstruction, the explicit approach, and the implicit approach for the multiple CoR-shifts case; the results are shown in Fig. 14. As expected, the standard approach converges rapidly to a local minimum that has a high objective value. More importantly, the implicit approach provides a better objective value than does the explicit approach. The reconstruction quality (right panel of Fig. 14) shows comparable performance between the explicit and the implicit approaches.

5.2 Noisy Data

In practice, the experimental data is known to be corrupted by noise. In this section, we test the performance of our proposed approach for contaminated data with respect to different levels of (mean zero) Gaussian noise. We start the test for the single CoR shift case as shown in Fig. 15, where the ground truth in this case is the standard MRI image. Similarly, the second and the third panels show the reconstructions without and with error recovery, respectively.

In this case, we compare the performance of the explicit and implicit approaches with the state-of-the-art alignment approach [13] provided by TomoPy 1.1.3 [12], a widely used tomographic data processing and image reconstruction library. The algorithm proposed in [13] is an iterative method that alternates between the alignment of the sinogram and the reconstruction of the object. The alignment step minimizes the misfit between the current sinogram and the newly simulated sinogram based on the current reconstruction. In Fig. 16, we illustrate the reconstruction results returned by the three methods as we increase the noise level in the experimental data; the columns represent an increase in the noise level from 4% to 22%. We observe that the explicit and implicit approaches are comparable with each other. However, the SSIM index for both approaches suggests that the implicit approach provides better reconstruction quality. This observation can be explained by the fact that although the explicit approach can provide a higher contrast image, its ability to remove noise outside of the sample region is not as good as that of the implicit approach. Overall, the performance of our proposed approaches outperforms the alternating algorithm [13] provided by TomoPy.

The last study is to test the performance of our proposed approach for recovering the multiple CoR shifts, as illustrated in Fig. 17. As measured by SSIM, the implicit approach clearly outperforms the other two approaches while the reconstructions of the alternating alignment do not approach the ground truth at all. In Fig. 18, we illustrate the reconstruction results returned by the three methods respectively as we increase the noise level in the same way as the single CoR shift case illustrated earlier. Again, the explicit and implicit approaches are able to provide a reasonable reconstruction, and the implicit approach in this case is able to provide a better reconstruction for the noisier case. Similarly to the single CoR shift case, the performance of our proposed approaches outperforms that of the alternating algorithm dramatically. Therefore, in the case of limited knowledge of the CoR shifts with noisy measurement data, the implicit approach works the best.

6 Conclusions

In this work, we propose a simultaneous inversion framework to address a common yet challenging experimental error in tomographic reconstruction, namely, a shift of the center of rotation. We also derive an analytical model to describe the general center-of-rotation shifts and how the sinograms with and without such shifts relate to each other. We derive an analytical model for the reconstruction framework to recover the center-of-rotation shifts together with the reconstruction of the object. The resulting optimization framework can be tuned to explicitly or implicitly reveal the true CoR locations. It is solved by a derivative-based optimization solver such as a truncated Newton algorithm.

Our numerical results show that when limited and noisy experimental data is available, significant improvements of the reconstruction quality can be achieved by performing the proposed simultaneous inversion. We compare our new approach with a state-of-the-art error reconstruction technique implemented in TomoPy and standard reconstruction techniques with a regularizer approach. For both comparisons, we observe superior results obtained by our approach. Given prior knowledge of how center-of-rotation shifts happen, explicit and implicit approaches can be chosen to maximize the performance, based on their different computational complexities. In particular, the implicit approach shows a stronger capability in more general scenarios such as noisy data.

Discussion: The simultaneous inversion we propose in this work is flexible and does not require extensive knowledge about the experimental system or the object. Even for the single CoR shift, the proposed approach does not require that the sample be completely contained in the field of view or that no nonlinear effects such as diffraction or refraction occur. It can be naturally extended to partially address more complicated experimental error such as per-scan drift. In parallel-beam tomography, the trajectories of the penetrating x-rays result in tomographic projections that form a stack of parallel 2D slices. Traditionally, 3D objects are recovered by reconstructing each slice independently [38]. Therefore, the extension of the proposed approach to a full 3D reconstruction is trivial in parallel-beam tomography. Alternatively, since each 2D slice shares the same CoR, the proposed approach can also be extended to a direct 3D reconstruction without adding an extra degree of freedom. One aspect we can further explore is to incorporate prior knowledge about the experimental system, such as the continuous dynamic of the center-of-rotation movement, into the joint framework as a regularizer on the error parameters. Overall, the challenge we are facing is the increased computational complexity. One way to mitigate this difficulty is to consider a parallel implementation that could be invoked inside a truncated-Newton solver or at the level of processing blocks of experimental data (e.g., per angle). A multigrid approach [24] could also benefit the computational performance.

Acknowledgments

We are grateful to Stefan Vogt, Si Chen, Doga Gursoy, and Francesco De Carlo for introducing us to this problem and for valuable discussions in the preparation of this paper. This material was based upon work supported by the U.S. Department of Energy, Office of Science, under contract DE-AC02-06CH11357.

Appendix A Properties of Radon Transform

We briefly describe a few properties of the Radon transform that explain the equivalence of reconstruction via translation and rotation.

Symmetry: $Rf(\tau,\theta)=Rf(-\tau,\theta+\pi)$.

Translation: Let $g(x-\Delta x,y-\Delta y)=f(x,y)$. Then

[TABLE]

where $\tau^{\prime}=\tau-\Delta x\cos\theta-\Delta y\sin\theta$.

[TABLE]

Rotation:

Let $g(r,\phi-\Delta\phi)=f(r,\phi)$ where $(r,\phi)$ is the polar coordinates of $(x,y)$. Then,

[TABLE]

[TABLE]

Appendix B Derivation of the CoR-Shifted Radon Transform

Here, we present a brief derivation of (2.2). Assume, without loss of generality, that the beam apparatus is initially positioned to measure the Radon transform along the direction corresponding to $\theta=0$. By Fig. 2, this amounts to integrating the function $f$ representing our object along the vertical lines $x=\tau$, $\tau\in\mathbb{R}$. To measure along the direction corresponding to a general angle $\theta$, we rotate the beam apparatus counterclockwise by $\theta$.222As described earlier, in the actual experimental setup we rotate the object, not the beam; from a mathematical point of view, it does not matter which we rotate. If the center of rotation is $(x_{\theta}^{*},y_{\theta}^{*})$, the line (beam) $x=\tau$ rotates to the line

[TABLE]

By definition, $Rf(\tau,\theta,x_{\theta}^{*},y_{\theta}^{*})$ is the integral of $f$ along $\ell(x_{\theta}^{*},y_{\theta}^{*})$. An easy calculation shows that, ignoring sign, this line is a distance

[TABLE]

from the origin. From the geometry of the rotation process we see that the line perpendicular to $\ell(x_{\theta}^{*},y_{\theta}^{*})$ is inclined at an angle $\theta$ above the (positive) horizontal axis. Hence,

[TABLE]

from which (2.2) follows immediately upon substituting in for $\tau^{\prime}$.

Appendix C Error Analysis for the Projection Alignment Process

In this section, we consider the error in the approximation $Rf(\tau,\theta,0,0)\approx\widetilde{Rf}(\tau,\theta,0,0)$ described in Sec. 2.3. Recall that $\widetilde{Rf}(\tau,\theta,0,0)$ incorporates both the Gaussian regularization of the Dirac delta and the discretization of the convolution integral using the discrete Fourier transform.

Moving to a more abstract setting, the question is fundamentally one of assessing the error in the approximation $f(x-P)\approx(f\ast\delta_{P,\sigma})(x)$, where $f$ is some compactly supported function on $\mathbb{R}$ and where we evaluate the right-hand side by discretizing the convolution integral in a particular way. We lose no generality in assuming that $P=0$, so we can actually analyze $f(x)\approx(f\ast\delta_{\sigma})(x)$, where we use the shorthand $\delta_{\sigma}=\delta_{0,\sigma}$. Since we wish to measure the error with respect to the uniform norm, we need only to consider the error at a single point $x$. We lose no generality in assuming that point is $x=0$, so we are led to the problem of understanding the error in

[TABLE]

again, under a specific discretization of the integral on the right-hand side, which in particular includes truncating the domain of the Gaussian $\delta_{\sigma}$.

Thus, we assume that $\operatorname*{supp}f\subset[-L,L]$ for some $L>0$, and we define

[TABLE]

which is $\delta_{\sigma}$ with its support truncated to $[-A,A]$. Our goal is to understand the error $E(\sigma,A,K)$ in the approximation to $f(0)$ obtained by using the $K$-point (composite) trapezoid rule to discretize the integral in (C.13) on the interval $[-M,M]$, where $M=L+A$:

[TABLE]

where $K=2N+1$ is odd333We make this choice solely for notational convenience. Everything we do applies for even $K$ as well. and $x_{k}=k(2M/K)$, $-N\leq k\leq N$. (This is mathematically equivalent to evaluating the convolution at one of the points $x_{k}$ using the discrete Fourier transform.) In particular, we want to know how $E(\sigma,A,K)$ depends on $\sigma$ and $K$.

$E(\sigma,A,K)$ comprises three sources of error. The first is the regularization error that we incur because $\sigma\neq 0$:

[TABLE]

The second is the truncation error that comes from our limiting the support of the Gaussian:

[TABLE]

The third is the discretization error due to our quadrature rule:

[TABLE]

We have

[TABLE]

We will bound $E(\sigma,A,K)$ by bounding each of the terms on the right-hand side of (C.14) in turn.

Assume that $f$ is twice continuously differentiable. Then, Taylor’s theorem (with Lagrange’s form of the remainder) enables us to write

[TABLE]

for $x\in\mathbb{R}$, where $\xi(x)$ is some function that satisfies $0\leq\xi(x)\leq x$ for $x>0$ and $x\leq\xi(x)\leq 0$ for $x<0$. From symmetry of $x\delta_{\sigma}(x)$, it thus follows that

[TABLE]

and hence we have

[TABLE]

for the regularization error.

The truncation error is even easier to handle. We have

[TABLE]

and so

[TABLE]

by one of the many available bounds for the complementary error function [6].

The discretization error is trickier. Because $\delta_{\sigma}^{\textrm{T}}$ is discontinuous on $[-M,M]$, we cannot apply most standard results, since these results assume continuity of the integrand. Instead, we appeal to the following result from [9], which bounds the error in the trapezoid rule approximation to the integral of a function of bounded variation.

Proposition C.1** ([9, Theorem 2.1]).**

If $f:[a,b]\to\mathbb{R}$ is of bounded variation, then

[TABLE]

*where $V(f)$ is the total variation of $f$ over $[a,b]$. *

This result immediately yields a bound of $(h/2)V(f)$ on the error in the composite trapezoid rule, where $h$ is the spacing between the grid points.

Because $f$ is twice continuously differentiable, it is of bounded variation, and so is $\delta_{\sigma}^{\textrm{T}}$. Hence, their product is of bounded variation as well. Moreover, we have

[TABLE]

(The norms and total variations are taken over $[-M,M]$.) Since

[TABLE]

and

[TABLE]

it follows that

[TABLE]

The grid spacing in our trapezoid rule approximation is $h=2M/K$. It follows that

[TABLE]

Combining (C.14)–(C.17), we arrive at the bound

[TABLE]

This shows that

[TABLE]

the implied limits in the big-O symbols being $\sigma\to 0$ and $K\to\infty$. This is sufficient to yield the desired qualitative understanding of the results presented in Sec. 2.3.