Alternative Fan-beam Backprojection and Adjoint Operators

TL;DR

This paper introduces new analytic formulas for fan-beam backprojection and its adjoint in tomography, derived from parallel tomography theorems, offering more robust performance in noisy data scenarios.

Contribution

It provides novel Bessel-Neumann series-based formulations for fan-beam backprojection, improving computational efficiency and noise robustness over traditional methods.

Findings

Achieves $O(N^{2.3729})$ computational complexity.

Demonstrates increased robustness with noisy data.

Validates formulations through numerical simulations.

Abstract

We present in this work alternative analytic formulations for the fan-beam tomographic backprojection operation and its associated adjoint transform in standard (equiangular) and linear (equidistant) detector geometries. The proposed formulations are obtained from a recent backprojection theorem in parallel tomography. Such formulations are written as a Bessel-Neumann series in the frequency domain that can be implemented as an $O(N^{2.3729})$ matrix multiplication. Proofs are provided together with numerical simulations compared with conventional fan-beam $O(N^{3})$ backprojection representations showing more robustness when dealing with highly noisy data.

Figures

Equations

$$\mathcal{R}f(t,\theta )={\int }_{{\mathbb{R}}^2}f(\mathbf{x})\delta (t-\mathbf{x}\cdot {\mathbf{\xi }}_{\theta })\mathrm{d}\mathbf{x},$$

$$\mathcal{R}{f}_{\mathsf{P}}(t,\theta )={\int }_0^{2\pi }{\int }_{{\mathbb{R}}_{+}}{f}_{\mathsf{P}}(\rho ,\varphi )\delta (t-\rho \cos (\theta -\varphi ))\rho \mathrm{d}\rho \mathrm{d}\varphi .$$

$$\mathcal{R}f(t,\theta )=\mathcal{R}f(-t,\theta +\pi ).$$

$$\forall u\in X,\forall v\in Y,⟨Tu,v{⟩}_Y=⟨u,T^{\ast }v{⟩}_X.$$

$${\mathcal{R}}^{\ast }p(\mathbf{x})={\int }_0^{2\pi }p(\mathbf{x}\cdot {\mathbf{\xi }}_{\theta },\theta )\mathrm{d}\theta .$$

$$\widehat{\mathcal{B}p}(\sigma {\mathbf{\xi }}_{\theta })=\frac{\hat{p}(\sigma ,\theta )}{\sigma },$$

$$w(\gamma ,\beta )=\mathcal{R}f(D\sin \gamma ,\beta +\gamma ),$$

$$g(s,\beta )=\mathcal{R}f(\frac{sD}{\sqrt{s^2+D^2}},\beta +\arctan \frac{s}{D}),$$

$$\begin{array}{rl}w(\gamma ,\beta )& =w(-\gamma ,\beta +2\gamma +\pi ),\\ g(s,\beta )& =g(-s,\beta +2\arctan \frac{s}{D}+\pi ).\end{array}$$

$${\mathcal{R}}_{\mathsf{s}}f(\gamma ,\beta )={\mathcal{M}}_{\mathsf{s}}\mathcal{R}f(\gamma ,\beta ),{\mathcal{R}}_{\ell }f(s,\beta )={\mathcal{M}}_{\ell }\mathcal{R}f(s,\beta ),$$

$$\begin{array}[]{ll}\mathcal{M}_{\mathsf{s}}p(\gamma,\beta)=p(D\sin\gamma,\beta+\gamma),\\ \mathcal{M}_{\ell}p(s,\beta)=p(\dfrac{sD}{\sqrt{s^{2}+D^{2}}},\beta+\arctan\dfrac{s}{D}).\end{array}$$

$$\begin{array}{rl}⟨{\mathcal{R}}_{\mathsf{s}}f,w{⟩}_{V_{\mathsf{s}}}& =⟨{\mathcal{M}}_{\mathsf{s}}\mathcal{R}f,w{⟩}_{V_{\mathsf{s}}}\\ & =⟨\mathcal{R}f,{\mathcal{M}}_{\mathsf{s}}^{\ast }w{⟩}_V=⟨f,{\mathcal{R}}^{\ast }{\mathcal{M}}_{\mathsf{s}}^{\ast }w{⟩}_U.\end{array}$$

$${\mathcal{R}}_{\mathsf{s}}^{\ast }={\mathcal{R}}^{\ast }{\mathcal{M}}_{\mathsf{s}}^{\ast }\text{and}{\mathcal{R}}_{\ell }^{\ast }={\mathcal{R}}^{\ast }{\mathcal{M}}_{\ell }^{\ast }.$$

$${\mathcal{R}}_{\mathsf{s}}^{\ast }w(\mathbf{x})={\int }_0^{2\pi }\frac{1}{L_{\beta }}w({\gamma }_{\beta },\beta )\mathrm{d}\beta ,$$

$$\begin{array}{rl}⟨{\mathcal{R}}_{\mathsf{s}}{f}_{\mathsf{P}},w{⟩}_{V_{\mathsf{s}}}& ={\int }_0^{2\pi }{\int }_0^{\pi }w(\gamma ,\beta ){\mathcal{R}}_{\mathsf{s}}{f}_{\mathsf{P}}(\gamma ,\beta )\mathrm{d}\gamma \mathrm{d}\beta \\ & ={\int }_0^{2\pi }{\int }_0^{\pi }w(\gamma ,\beta ){\int }_0^{2\pi }{\int }_{{\mathbb{R}}_{+}}{f}_{\mathsf{P}}(\rho ,\varphi )\delta (L_{\beta }\sin ({\gamma }_{\beta }-\gamma ))\rho \mathrm{d}\rho \mathrm{d}\varphi \mathrm{d}\gamma \mathrm{d}\beta ,\end{array}$$

$$⟨{\mathcal{R}}_{\mathsf{s}}{f}_{\mathsf{P}},w{⟩}_{V_{\mathsf{s}}}={\int }_0^{2\pi }{\int }_{{\mathbb{R}}_{+}}f(\rho ,\varphi ){\int }_0^{2\pi }{\int }_0^{\pi }\frac{1}{L_{\beta }}w(\gamma ,\beta )\delta (\sin ({\gamma }_{\beta }-\gamma ))\mathrm{d}\gamma \mathrm{d}\beta \rho \mathrm{d}\rho \mathrm{d}\varphi .$$

$${\mathcal{R}}_{\ell }^{\ast }g(\mathbf{x})=\frac{1}{D}{\int }_0^{2\pi }\frac{1}{U_{\beta }}\sqrt{s_{\beta }^2+D^2}g(s_{\beta },\beta )\mathrm{d}\beta ,$$

$$U_{\beta }=\frac{D^2-{\mathbf{r}}_{\beta }\cdot \mathbf{x}}{D^2},s_{\beta }=\frac{\mathbf{x}\cdot {\mathbf{\xi }}_{\beta }}{U_{\beta }}.$$

$${\mathcal{B}}_{\mathsf{s}}w(\mathbf{x})={\int }_0^{2\pi }\frac{1}{L_{\beta }^2}w({\gamma }_{\beta },\beta )\mathrm{d}\beta ,{\mathcal{B}}_{\ell }g(\mathbf{x})={\int }_0^{2\pi }\frac{1}{U_{\beta }^2}g(s_{\beta },\beta )\mathrm{d}\beta ,$$

$${\mathcal{M}}_{\mathsf{s}}^{\ast }w(t,\theta )=w(\gamma (t),\theta -\gamma (t)){\mathsf{J}}_{\mathsf{s}}(t),$$

$${\mathcal{M}}_{\ell }^{\ast }g(t,\theta )=g(s(t),\theta -\gamma (t)){\mathsf{J}}_{\ell }(t),$$

$$\gamma (t)=\arcsin \frac{t}{D},s(t)=\frac{tD}{\sqrt{D^2-t^2}},{\mathsf{J}}_{\mathsf{s}}(t)=\frac{1}{\sqrt{D^2-t^2}},{\mathsf{J}}_{\ell }(t)=\frac{D^3}{(D^2-t^2{)}^{3/2}}.$$

$$\widehat{{\mathcal{R}}_{\mathsf{s}}^{\ast }w}(\sigma {\mathbf{\xi }}_{\theta })=\widehat{{\mathcal{R}}^{\ast }{\mathcal{M}}_{\mathsf{s}}^{\ast }w}(\sigma ,\theta )=\frac{2}{\sigma }\widehat{{\mathcal{M}}_{\mathsf{s}}^{\ast }w}(\sigma ,\theta ).$$

$$c_n(\theta )=\frac{1}{2\pi }{\int }_0^{2\pi }z(\gamma ,\theta )e^{-in\gamma }\mathrm{d}\gamma ,$$

$$b_n=2\pi [c_n+(-1{)}^n\overline{c_n}]\text{for}n\ge 1,\text{and}{b}_0=2\pi c_0.$$

$$\widehat{{\mathcal{M}}_{\mathsf{s}}^{\ast }w}(\sigma ,\theta )=\sum _{n=0}^{\infty }{b}_n(\theta )J_n(D\sigma ),$$

$$\begin{array}{rl}\widehat{{\mathcal{M}}_{\mathsf{s}}^{\ast }w}(\sigma ,\theta )& ={\int }_{\mathbb{R}}w(\gamma (t),\theta -\gamma (t))e^{-it\sigma }{\mathsf{J}}_{\mathsf{s}}(t)\mathrm{d}t\\ & ={\int }_0^{2\pi }w(\gamma ,\theta -\gamma )e^{-iD\sigma \sin \gamma }\mathrm{d}\gamma ={\int }_0^{2\pi }z(\gamma ,\theta )e^{-iD\sigma \sin \gamma }\mathrm{d}\gamma .\end{array}$$

$$\begin{array}{rl}\widehat{{\mathcal{M}}_{\mathsf{s}}^{\ast }w}(\sigma ,\theta )& =\sum _n{c}_n(\theta ){\int }_0^{2\pi }{e}^{i[n\gamma -D\sigma \sin \gamma ]}\mathrm{d}\gamma \\ & =2\pi \sum _n{c}_n(\theta )J_n(D\sigma ),\end{array}$$

$$\widehat{{\mathcal{R}}_{\mathsf{s}}^{\ast }w}(\sigma {\mathbf{\xi }}_{\theta })=\frac{2}{\sigma }\sum _{n=0}^{\infty }{b}_n(\theta )J_n(D\sigma ),$$

$$\widehat{{\mathcal{B}}_{\mathsf{s}}w}(\sigma {\mathbf{\xi }}_{\theta })=\widehat{\mathcal{B}{\mathcal{M}}_{\mathsf{s}}^{-1}w}(\sigma ,\theta )=\frac{1}{\sigma }\widehat{{\mathcal{M}}_{\mathsf{s}}^{-1}w}(\sigma ,\theta ).$$

Full text

Keywords: Fan-beam, backprojection, tomographic reconstruction, Radon transform.

AMS subject classifications: 44A12, 41A58.

Published in: J. of Inverse and Ill-Posed Problems, 2023.

Corresponding author: patricio.guerrero@kuleuven.be

1. Introduction

Fan-beam tomographic measurements are used in different modalities of non-destructive imaging, as those obtained using X-rays. A typical tomographic device using fan-beam geometry is shown in Figure 1. This is a widely used and known technique, and there are many reconstruction algorithms for this configuration. After being generated with a given aperture angle and a fixed source-detector distance, the wavefront hits the sample originating a signal on the detector. Different propagation regimes can be considered with a varying distance [13, 15], although we will consider a pure mathematical signal idealized as the fan-beam Radon transform of the given object.

For the tomographic parallel case by means of the classical backprojection formulation for image reconstruction, we can use the recent backprojection slice theorem formulation [16]. It is a formula that reduces a complete backprojection of a 2D sinogram from a computational cost of $O(N^{3})$ to $O(N^{2}\log N)$ with $N\times N$ pixels in the final reconstructed image and $N$ tomographic rotations. In this work, we take advantage of that formulation for two other popular fan-beam geometries, that are a) equispaced angles within the fan with uniform size $\Delta\gamma$, and b) equispaced points in a linear detector with a mesh step $\Delta s$. As indicated in [19], we refer to each acquisition as standard fan and linear fan, respectively, and these are illustrated in Figure 1.(a) and 1.(b). The linear case is easier to be implemented at a synchrotron beamline with an appropriate optical setup, whereas the former is more used with medical or industrial X-rays scans. For completeness, we discuss both standard and linear cases.

Related work, contributions and overview

Further work has been done after the conventional $O(N^{3})$ backprojection strategies for fan-beam sinograms were established in [9], where both the fan-beam filtered backprojection (FBP) algorithm and a fan-to-parallel rebinning algorithm is formalized. In [28] a generalized Fourier slice theorem (FST) for fan-beam geometries is introduced presenting the same $O(N^{3})$ computational complexity as a conventional backprojection in the parallel case. In addition it requires a rebinning operation on the measured data from fan to parallel geometry and then it generalize the FST. Further relations on the frequency domain are obtained in [10] where a fan-to-parallel rebinning is also necessary and the Gerchberg-Papoulis iterative algorithm is used to obtain the final reconstruction; the method is mentioned to be computationally expensive although the computational complexity is not explicitly described. A different approach is established in [3] without direct rebinning in the image domain but in the frequency domain of the reconstruction; however, a weighted backprojection of the data is needed and no algorithm with lower cost than $O(N^{3})$ is provided to this end. In [25], a series formulation for the fan-beam FBP algorithm is presented where the backprojection is written as the first order approximation of the FBP series. Using a hierarchical approach, [7] presents an $O(N^{2}\log N)$ algorithm based on efficient sampling of the fan-beam geometry and without rebinning to parallel projections, although a series of intermediate interpolations are needed on every level of the hierarchical structure.

Fast algorithms are either based on deriving a mathematical expression relying on low complexity auxiliary transformations, e.g., the fast Fourier transform, or hierarchically, taking a divide-and-conquer approach [6]. While the method in [7] relies on the latter approach, this work relies on the former. We present here an analytic formulation for the fan-beam backprojection operation as well as for the associated adjoint transform, considering both standard (Theorem 2) and linear (Theorem 4) geometries. Such formulations can be implemented numerically as a multiplication of two matrices, allowing the use of fast matrix multiplication algorithms [11, 26] resulting in a computational complexity of $O(N^{2.3729})$. It is worth to recall then that our method is computationally more expensive than the $O(N^{2}\log N)$ hierarchical method of [7] and the purpose of this work is an alternative analytic and simple formulation of the fan-beam backprojection and adjoint transforms still less expensive than the usual $O(N^{3})$ formulation.

This work is organized as follows. Section 2 presents a review of conventional fan-beam backprojection operators for both standard and linear geometries and the relationship with their associated adjoint transforms, which are defined as bounded operators on Hilbert spaces. In section 3 we rewrite the rebinning strategy in terms of the inverse rebinning operator as a two-step backprojection formula. In this section, the adjoint transforms of both geometries are also formalized as a composition of two operators and therefore, we show that this formulation is equivalent to the one in section 2 from the uniqueness of the adjoint of bounded operators. Section 4 presents our fan-beam adjoint and backprojection theorems as a Bessel-Neumann series representation (referred here as BN series method) and their proofs, while section 5 numerically supports our method with simulated data where the backprojection is computed as a straightforward matrix multiplication based on the BN series method and its performance is compared with the conventional rebinning strategy.

2. Fan-beam backprojections

2.1. The 2D Radon transform

The two-dimensional parallel Radon transform is defined [5] as the linear operator $\mathcal{R}\colon U\to V$ where $U$ is the Schwartz space $\mathcal{S}(\mathbb{R}^{2})$ of rapidly decreasing functions defined on $\mathbb{R}^{2}$, so-called feature space; and $V$, the sinogram space, is the Schwartz space defined on the unit cylinder $\mathcal{S}(\mathbb{R}\times\mathbb{R}/2\pi\mathbb{Z})$ where $\mathbb{R}/2\pi\mathbb{Z}$ is the quotient group of reals modulo $2\pi$.

Let $(t,\theta)\in\mathbb{R}\times\mathbb{R}/2\pi\mathbb{Z}$ and $f\in U$. $\mathcal{R}f$ is defined as the integral operator

[TABLE]

where $\delta$ is the Delta distribution and $\bm{\xi}_{\theta}=(\cos\theta,\sin\theta)^{T}$. We will denote by $f_{\mathsf{P}}\in U$ a feature function acting on $\bm{x}=(\rho,\phi)\in\mathbb{R}_{+}\times[0,2\pi[$ expressed in polar coordinates. $\mathcal{R}f_{\mathsf{P}}$ is then written as

[TABLE]

Note that $\mathcal{R}f$ is an even function on the unit cylinder, that is,

[TABLE]

Adjoint and backprojection.

Let $T\colon X\to Y$ be a bounded linear operator between Hilbert spaces $X$ and $Y$. Its adjoint operator is defined as the operator $T^{*}\colon Y\to X$ that verifies

[TABLE]

The adjoint and backprojection operators of the Radon transform are closely related as follows [5]. Expanding the definitions of $U$ and $V$ to $L^{2}$ spaces [12], the adjoint of $\mathcal{R}$ is the operator $\mathcal{R}^{*}\colon V\to U$ defined for $p\in V$ and $\bm{x}\in\mathbb{R}^{2}$ by

[TABLE]

With standard $L^{2}$ inner products on $U$ and $V$, it is clear that $\mathcal{R}^{*}$ verifies (4).

In the FBP algorithm [19], the backprojection operation is applied after filtering the data with an adequate convolutional filter. Such backprojection operator, denoted $\mathcal{B}$ is defined by limiting the last integral to $[0,\pi]$, and from the symmetry property (3) we have $\mathcal{R}^{*}\equiv 2\mathcal{B}$.

Computing $b=\mathcal{B}p$ could be extremely expensive for discrete versions of the sinogram $p$, when $\bm{x}$ covers a domain with a large number of points (pixels in practice) and also $(t,\theta)$ covers a large number of pixels and a variety of angles (according to Crowther’s criterion [4]). This is the case for synchrotron tomographic projections using high-resolution detectors with more than $2048\times 2048$ pixels and more than $2048$ angles [17]. Recently [16], a low-complexity formulation for computing $b$ was obtained in the frequency domain using polar coordinates, the Backprojection Slice Theorem (BST). Such theorem reads for $p\in V$ and $\sigma>0$,

[TABLE]

where $\widehat{\mathcal{B}p}$ is the 2D Fourier transform of $\mathcal{B}p$ in polar frequency coordinates and $\hat{p}$ is the 1D Fourier transform of $p$ with respect to $t$. The action of $\{\mathcal{R},\mathcal{B}\}$ is presented in Figure 2. It is a well known fact that $\{p,f\}$ are related through the FST [19, 5], while $\{p,b\}$ through the BST [16].

2.2. Fan-beam geometries

There exist two main fan-beam parameterizations for the Radon transform [9], the first $\mathcal{R}_{\mathsf{s}}\colon U\to V_{\mathsf{s}}$ is referred as the standard fan-beam transform, where the domain of sinograms $V_{\mathsf{s}}$ lies within the set $\mathbb{R}/\pi\mathbb{Z}\times\mathbb{R}/2\pi\mathbb{Z}$. The second case, $\mathcal{R}_{\ell}\colon U\to V_{\ell}$ is referred as the linear fan-beam transform, with sinograms domain $V_{\ell}$ varying within the set $\mathbb{R}\times\mathbb{R}/2\pi\mathbb{Z}$.

Let us consider $(\gamma,\beta)\in\mathbb{R}/\pi\mathbb{Z}\times\mathbb{R}/2\pi\mathbb{Z}$ a pair of angles and $s\in\mathbb{R}$. Acting on a feature function $f\in U$, the linear operators $\mathcal{R}_{\mathsf{s}}$ and $\mathcal{R}_{\ell}$ are defined respectively as

[TABLE]

and

[TABLE]

where $D>0$ is the distance of the focal point source to the center of rotation located at the origin, see Figure 1. Hence $w\in V_{\mathsf{s}}$ and $g\in V_{\ell}$. The following symmetry relationships hold:

[TABLE]

Since both operations represent a change of variables in the classical parallel sinogram, we can use the following notation:

[TABLE]

where $\mathcal{M}_{\mathsf{s}}\colon V\to V_{\mathsf{s}}$ and $\mathcal{M}_{\ell}\colon V\to V_{\ell}$ are rebinning operators acting on a parallel sinogram $p$, respectively as

[TABLE]

From the experimental point of view, we assume that the origin is positioned at the center of rotation, where the sample is centered. In both geometries, for simplicity and without any loss of generality, a virtual detector will be considered centered at the origin. In the standard geometry, the angle $\gamma$ from the central ray indicates a position in the circle centered at the beam source with radius $D$, this circle acts as the detector in this geometry and then it’s clearly $2\pi$-periodic on $\gamma$. Whereas in the linear geometry, $s$ is the signed height on the linear detector perpendicular to the central ray and centered at the origin, see Figure 1.

The change of variables operation $\mathcal{M}_{\mathsf{s}}$ is depicted in Figure 2. $\mathcal{M}_{\mathsf{s}}$ and $\mathcal{M}_{\ell}$ are well defined operators because $\beta+\gamma$ and $\beta+\arctan\dfrac{s}{D}$ both belong to $\mathbb{R}/2\pi\mathbb{Z}$. Since it is true that $\mathcal{R}_{\mathsf{s}}=\mathcal{M}_{\mathsf{s}}\mathcal{R}$, a conventional functional relation for space $U$ and $V_{\mathsf{s}}$ provides us with the following statement,

[TABLE]

Therefore, it follows that $\mathcal{R}_{\mathsf{s}}^{*}\ =\mathcal{R}^{*}\mathcal{M}_{\mathsf{s}}^{*}$. In this work we propose a Fourier approach for equations

[TABLE]

2.3. Integral representation of the adjoint transforms

Considering both fan-beam geometries, we provide two main results for both the adjoint and backprojection operators of $\mathcal{R}_{\mathsf{s}}$ and $\mathcal{R}_{\ell}$. In the following we denote by $\bm{r}_{\beta}$ the Cartesian coordinates of the source.

Lemma 1**.**

The operator $\mathcal{R}_{\mathsf{s}}^{*}\colon V_{\mathsf{s}}\to U$, defined for $w\in V_{\mathsf{s}}$ and $\bm{x}\in\mathbb{R}^{2}$ by

[TABLE]

is the the adjoint of $\mathcal{R}_{\mathsf{s}}$ in the sense of (4). Here, $L_{\beta}=\left\|\bm{r}_{\beta}-\bm{x}\right\|_{2}$ is the Euclidean distance of the source with the backprojected position $\bm{x}$ and $\gamma_{\beta}$ is the angle of such point within the fan, i.e., $\cos\gamma_{\beta}=\dfrac{1}{DL_{\beta}}\bm{r}_{\beta}\cdot(\bm{r}_{\beta}-\bm{x})$.

Proof.

Let $f_{\mathsf{P}}\in U$ and $w\in V_{\mathsf{s}}$. Using the polar representation (2) of $\mathcal{R}$ and the relationship ${\rho\cos(\beta+\gamma-\phi)=L_{\beta}\sin(\gamma_{\beta}-\gamma)}$ (see [9]), then

[TABLE]

which, from Fubini’s theorem, becomes

[TABLE]

Therefore, since $\delta(\sin(\gamma_{\beta}-\gamma))=\delta(\gamma_{\beta}-\gamma)$ for $\gamma\in[0,\pi[$ and a fixed $\gamma_{\beta}$, after performing the above $\gamma$-integration, $\mathcal{R}_{\mathsf{s}}^{*}$ defined in (10) verifies (4). ∎

Lemma 2**.**

The operator $\mathcal{R}_{\ell}^{*}\colon V_{\ell}\to U$, defined for $g\in V_{\ell}$ and $\bm{x}\in\mathbb{R}^{2}$ by

[TABLE]

is the adjoint of $\mathcal{R}_{\ell}$ in the sense of (4). $U_{\beta}$ is the ratio of the scalar projection of $\bm{r}_{\beta}-\bm{x}$ on the central ray to the source-origin distance, and $s_{\beta}$ is the corresponding height $s$ for $\bm{x}$ at a source angle $\beta$, i.e.,

[TABLE]

Proof.

As in the proof of Lemma 1, this follows directly from (4) and the polar representation (2). ∎

The relationship between the adjoint transforms in (10) and (11) and the backprojection operators, denoted here by $\mathcal{B}_{\mathsf{s}}$ and $\mathcal{B}_{\ell}$ respectively for the standard and linear geometry follows. The backprojection operators are applied e.g., on the FBP algorithm as in the parallel case after the filtering process [9]. Such operators are respectively defined as the weighted integrals

[TABLE]

therefore, the difference with the adjoint transforms lies in the weight factors.

3. Two-step backprojection formulas

The two-step adjoint formulations for each geometry, namely $\mathcal{R}_{\mathsf{s}}^{*}=\mathcal{R}^{*}\mathcal{M}_{\mathsf{s}}^{*}$ and $\mathcal{R}_{\ell}^{*}=\mathcal{R}^{*}\mathcal{M}_{\ell}^{*}$ derived from (8) will be detailed here. It is clear that the adjoint of the rebinning operators $\{\mathcal{M}_{\mathsf{s}},\mathcal{M}_{\ell}\}$ plays an important role in such formulations.

Lemma 3** (Adjoint of a rebinning operator).**

Let $\mathcal{M}\mathrel{\mathop{\mathchar 58\relax}}X\to Y$ be a bijective and differentiable with continuous inverse rebinning operator between two Hilbert spaces $X$ and $Y$. Its adjoint $\mathcal{M}^{*}\colon Y\to X$ is then given by $\mathcal{M}^{*}=\mathsf{J}\,\mathcal{M}^{-1}$, where $\mathsf{J}$ is the Jacobian determinant of the rebinning operation and $\mathcal{M}^{-1}$ its inverse rebinning.

Proof.

If $u\in X$ and $v\in Y$, a direct verification of $\langle\mathcal{M}u,v\rangle_{Y}=\langle u,\mathcal{M}^{*}v\rangle_{X}$ through an integral representation and a change of variables will prove the result. Bijectivity of $\mathcal{M}$ is needed to preserve the range of integration in both domains of $X$ and $Y$ and for the inverse to exist. ∎

Lemma 4**.**

The adjoint operators of $\mathcal{M}_{\mathsf{s}}$ and $\mathcal{M}_{\ell}$ are respectively defined for $w\in V_{\mathsf{s}}$ and $g\in V_{\ell}$ by

[TABLE]

[TABLE]

where

[TABLE]

Proof.

As both operators are bijective [19] and continuously differentiable, this is an immediate application of Lemma 3. ∎

The formal adjoints of $\mathcal{R}_{\mathsf{s}}$ and $\mathcal{R}_{\ell}$ are presented in Lemmas 1 and 2. The problem with these formulations is the difficulty for a low-cost implementation algorithm. To circumvent this problem, we use the fact that $\{\mathcal{R},\mathcal{R}^{*}\}$ are bounded operators between Hilbert spaces $U$ and $V$, here considered to be $L^{2}$ spaces. Hence, we obtain the following result.

Theorem 1**.**

Considering the two fan-beam geometries, the adjoint operators of $\mathcal{R}_{\mathsf{s}}$ and $\mathcal{R}_{\ell}$, given explicitly by expressions (10) and (11), are also given respectively by $\mathcal{R}_{\mathsf{s}}^{*}=\mathcal{R}^{*}\mathcal{M}_{\mathsf{s}}^{*}$ and $\mathcal{R}_{\ell}^{*}=\mathcal{R}^{*}\mathcal{M}_{\ell}^{*}$.

Proof.

This is an immediate consequence of the uniqueness of the adjoint for bounded operators on Hilbert spaces [21]. In fact, since $\mathcal{R}^{*}$ and $\mathcal{M}_{\mathsf{s}}^{*}$ are bounded, the composition is also bounded. The same applies for $\mathcal{M}_{\ell}^{*}$. ∎

To conclude this section, a straightforward two-step rebinning formulation for the backprojection operators is provided here as $\mathcal{B}_{\mathsf{s}}=\mathcal{B}\mathcal{M}_{\mathsf{s}}^{-1}$ and $\mathcal{B}_{\ell}=\mathcal{B}\mathcal{M}_{\ell}^{-1}$ having an evident justification. The difference mainly lies on the Jacobian determinants $\mathsf{J}_{\mathsf{s}}$ and $\mathsf{J}_{\ell}$ as we can write $\mathcal{M}_{\mathsf{s}}^{*}=\mathcal{M}_{\mathsf{s}}^{-1}\mathsf{J}_{\mathsf{s}}$ and $\mathcal{M}_{\ell}^{*}=\mathcal{M}_{\ell}^{-1}\mathsf{J}_{\ell}$. From the adjoint transform expressed as the two-step operation in Theorem 1, we derive our main result in the following section.

4. Main backprojection theorems

4.1. Standard geometry

The main result of this work is based on a Fourier approach for equations (9) and then we make use of the BST formula (5) and the identity $\mathcal{R}^{*}\equiv 2\mathcal{B}$ for a fixed $\theta$ and $\sigma>0$ in the form

[TABLE]

In [20, eq. 2.5.6], a formula relating a parallel sinogram and its backprojection is obtained, similarly to the BST, with a kernel $1/\sigma^{n-1}$ for an arbitrary dimension $n$, also in the frequency domain. This formulation requests however, contrary to the BST, the sinogram $p$ to be the in the range of the Radon operator $\mathcal{R}$, i.e., it should exist some $f\in U$ such that $p=\mathcal{R}f$. This is not necessarily true in the derivation of (15), but on the other hand we can safely apply the BST here.

We start announcing our Fourier-based fan-beam adjoint theorem for standard geometry, and then we show how it is easily adapted to the linear case. From (15), we observe the need of computing $\widehat{\mathcal{M}_{\mathsf{s}}^{*}w}$ that is written as a series representation in the following Lemma.

Lemma 5**.**

Given a standard fan-beam sinogram $w\in V_{\mathsf{s}}$, first define $z\in V_{\mathsf{s}}$ by $z(\gamma,\theta)=w(\gamma,\theta-\gamma)$ for all $(\gamma,\theta)\in\mathbb{R}/\pi\mathbb{Z}\times\mathbb{R}/2\pi\mathbb{Z}$. As $z$ is $2\pi$-periodic with respect to $\gamma$ we can write its Fourier series expansion with the Fourier coefficients

[TABLE]

from where we define the coefficients

[TABLE]

Then we have a Bessel-Neumann (BN) series description of $\widehat{\mathcal{M}_{\mathsf{s}}^{*}w}$ as

[TABLE]

where $(J_{n})$ is a sequence of Bessel functions of the first kind.

Proof.

Given $w\in V_{\mathsf{s}}$, from (12) and (14) the 1D Fourier transform of $\mathcal{M}_{\mathsf{s}}^{*}w(t,\theta)$ with respect to $t$ is

[TABLE]

After expanding $Z$ with $c_{n}$ we have

[TABLE]

from where the result claims using (16) and the property $J_{-n}(x)=(-1)^{n}J_{n}(x)$. ∎

We can now state our fan-beam adjoint theorem for standard geometries.

Theorem 2** (Standard fan-beam adjoint Theorem).**

The adjoint operator ${\mathcal{R}_{\mathsf{s}}^{*}\colon V_{\mathsf{s}}\to U}$ acting on $w\in V_{\mathsf{s}}$ can be expressed in the Fourier domain as the series

[TABLE]

where $(\sigma,\theta)\in\mathbb{R}_{+}\times\mathbb{R}/2\pi\mathbb{Z}$ and the coefficients $b_{n}$ are computed as in Lemma 5.

Proof.

This is an immediate consequence of (15) and Lemma 5. ∎

Remark 1**.**

The backprojection operator $\mathcal{B}_{\mathsf{s}}$ can also be written as a BN series. In fact, from $\mathcal{B}_{\mathsf{s}}=\mathcal{B}\mathcal{M}_{\mathsf{s}}^{-1}$ we have

[TABLE]

Then, the Jacobian determinant $\mathsf{J}_{\mathsf{s}}$ in the derivation (17) needs to be compensated. Due to $\mathrm{d}t=D\cos\gamma\,\mathrm{d}\gamma$ in (17), we define $\dot{z}=\mathcal{A}w$ by

[TABLE]

and then we compute the Fourier coefficients $(\dot{b}_{n})$ from $\dot{z}$ in the same way as in (16). Following the same derivation of Theorem 2, we have the representation

[TABLE]

Moreover, from the operation $\beta=\theta-\gamma$ we note that to perform the backprojection with the BN series, we only need a short-scan sinogram where $\beta\in[0,\pi+2\bar{\gamma}[$, with $\bar{\gamma}=\arcsin{\bar{x}/D}$ the upper bound of $\gamma$ such that the beam intersects the domain of the phantom bounded by $\mathinner{\!\left\lVert\bm{x}\right\rVert}_{2}\leq\bar{x}$.

Remark 2**.**

The backprojected image obtained with BN series (20) assumes a null DC component due to the kernel $1/\sigma$. However, (20) can still be used to reconstruct an arbitrary object $f$ with the correct DC component $\kappa$ by estimating it using the volume conservation property of the Radon transform [19] that states

[TABLE]

Then, using the rebinning operator $\mathcal{M}_{\mathsf{s}}$, $\kappa$ can be estimated from fan-beam projections $\mathcal{R}_{\mathsf{s}}f$ by

[TABLE]

The last integral is computed with a fixed $\theta$ but an average over $\theta$ is also possible giving a more robust $\kappa$ with e.g., noisy projections.

Remark 3**.**

Convergence of the BN series (20) is guaranteed due to the fact that $J_{n}(x)$ has a pointwise convergence to [math] since $|J_{n}(x)|\leq|\frac{1}{2}x|^{n}/n!$ for all $x\in\mathbb{R}_{+}$ [1, eq. 9.1.62].

4.2. Linear geometry

The adjoint operator for the linear fan-beam transform follows. Switching between fan-beam geometries requires only a one-dimensional interpolation on the first variable. In fact, taking $g\in V_{\ell}$, the rebinning operator on the first variable $\mathcal{L}\colon V_{\ell}\to V_{\mathsf{s}}$ acting as

[TABLE]

provides a sinogram $w$ on the space $V_{\mathsf{s}}$. Such bijective operator, according to Lemma 3, has an adjoint operator given by

[TABLE]

where

[TABLE]

The following Theorem enables us to provide a adjoint algorithm for the linear fan-beam geometry.

Theorem 3**.**

The adjoint operators $\mathcal{R}_{\mathsf{s}}^{*}$ and $\mathcal{R}_{\ell}^{*}$ are related, for all $g\in V_{\ell}$ by

[TABLE]

where $\tau(\gamma)=D\sec^{2}\gamma$.

Proof.

From (21) we obtain $\mathcal{R}_{\mathsf{s}}^{*}=(\mathcal{L}\mathcal{R}_{\ell})^{*}=\mathcal{R}_{\ell}^{*}\mathcal{L}^{*},$ from where $\mathcal{R}_{\ell}^{*}=\mathcal{R}_{\mathsf{s}}^{*}(\mathcal{L}^{*})^{-1}$ follows. Using (22) and the fact that $(\mathcal{L}^{*})^{-1}=\mathsf{J}(s(\gamma))^{-1}\mathcal{L}$ with $s(\gamma)=D\tan\gamma$, equation (23) is obtained. ∎

Therefore, our resulting formulation for the linear fan-beam adjoint transform is obtained by the following construction.

Theorem 4** (Linear fan-beam adjoint Theorem).**

The operator $\mathcal{R}_{\ell}^{*}\colon V_{\mathsf{s}}\to U$ acts on $g\in V_{\ell}$ as the following two-steps operation.

1. (i)

*Let $\dot{w}(\gamma,\beta)=\tau(\gamma)w(\gamma,\beta)$ with $w=\mathcal{L}g$ be in $V_{\mathsf{s}}$. * 2. (ii)

The adjoint operator is obtained as the BN series (18) where the coefficients $(b_{n})$ in (16) are obtained from $\dot{z}(\gamma,\theta)=\dot{w}(\gamma,\theta-\gamma)$.

Proof.

Each step is justified respectively by Theorems 3 and 2. ∎

Finally, the backprojection operator $\mathcal{B}_{\ell}$ can also be computed in the same way as in Remark 1, by first writing $\mathcal{B}_{\ell}g=\mathcal{B}_{\mathsf{s}}\mathcal{L}g$ and then computing $\mathcal{B}_{\mathsf{s}}$ using (20).

Figure 3 illustrates the action of our method, where a backprojection $b\in U$ is obtained in two-steps. The unknown backprojection $\mathcal{B}_{\ell}g$ is obtained in the frequency domain from (18), using polar coordinates by two main operations to obtain sinogram $\dot{z}\in V_{\mathsf{s}}$ through $\mathcal{L}$ defined in (21) (interpolation on the first variable) and $\mathcal{A}$ defined in (19) (interpolating on the second variable).

4.3. Mapping to Cartesian coordinates

A backprojected fan-beam sinogram using Theorems 2 or 4 respectively for standard and linear geometry is obtained in polar coordinates $(\sigma,\theta)\in\mathbb{R}_{+}\times\mathbb{R}/2\pi\mathbb{Z}$ in the Fourier domain that needs to be mapped to Cartesian coordinates before applying the inverse Fourier transform. This problem is also present in parallel geometry if using the Fourier reconstruction method based on the FST [19] or using the BST approach [16]. The evident strategy is a straightforward (linear) interpolation with the expense of important errors at higher frequencies.

An error estimate related to this polar-to-Cartesian mapping for the parallel case using the FST is provided in [18, sec. V.2]. If $\tilde{f}$ is the reconstruction of $f$ using the Fourier method, the error $\|f-\tilde{f}\|_{L_{2}}$ related to the mentioned mapping can be bounded by $c\|f\|_{\mathcal{H}}$ with some constant $c>0$ and a convenient norm on a Sobolev space $\mathcal{H}$. If a parallel sinogram $g(t,\theta)$ is backprojected within this error bound, then (not conversely, as explained below) it verifies the optimal sampling conditions $\{\Delta t\leq\pi/\Omega,\Delta\theta\leq\pi/(\rho\Omega)\}$ for an object with support included in the disk of radius $\rho$ and $\Omega$ the highest frequency following the Nyquist sampling criterion. Equivalent sampling conditions are also available for fan-beam sinograms in [19, sec. IV.3] as $\{\Delta\gamma\leq\pi/(D\Omega),\Delta\beta\leq\pi(D+\rho)/(\rho D\Omega)\}$ in standard geometry and the obvious equivalent in the linear case. It is worth to mention that the highest frequency $\Omega$ is obtained in such references with the Fourier Transform defined with the normalization as in (5). Conversely, it is also shown in [18] that the optimal angular sampling is enough to guarantee the error bound but not the radial sampling in polar frequencies; this can be solved by oversampling by $2$ the Fourier domain with zero-padding. In our case, this is largely verified as the oversampling is performed anyway by a factor at least $2$ until convergence of the BN series. Under these conditions, linear interpolation in the Fourier domain is enough within classical interpolations methods as it will be shown in Figure 7(a) where linear interpolation is compared with quadratic splines where results presented negligible differences. In the parallel case, the BST approach is shown to reach competitive results to other $O(N^{2}\log N)$ methods and overcome them in some scenarios with only linear interpolation from polar to Cartesian coordinates and an adequate zero-padding [16].

More importantly, advanced polar-to-Cartesian mapping strategies in the context of the Fourier method with FST have shown to compete or even overcome FBP strategies while still presenting the lower $O(N^{2}\log N)$ complexity. For example, a gridding strategy [19, sec. V.2.2] can be used by weighting the polar data with an adequate function and then the mapping to Cartesian coordinates is performed by convolution with such weight function. This method handles better interpolations errors than polynomial-based interpolation and reach equivalent results as FBP [14]. Although no error bounds are presented in such work, the claim is proved by numerical results. The gridding algorithm is also implemented in the academic reconstruction library tomopy and applied to reconstruct synchrotron tomographic data in different synchrotron facilities [8, 14] instead of FBP. Finally, another strategy is to regularize the passage to Cartesian coordinates using total variation regularization based on the approach developed in [27] for parallel tomography with FST, still with an $O(N^{2}\log N)$ complexity. This technique is motivated by compressed sensing and sparsity priors successfully developed for magnetic resonance imaging (MRI) reconstructions where there is also Fourier data on polar coordinates and it is shown that even with strongly undersampled Fourier data, significantly less polar samples than the optimal conditions given above, high quality reconstructions can be obtained [2]. These two mentioned strategies are available to be used directly on our reconstruction method for fan-beam tomography since the BN series provides polar data on the Fourier domain of the feature object as FST/BST for parallel tomography or MRI does.

5. Numerical results

Simulations are concentrated on the linear case since the standard case is one rebinning operation less than the linear, so results would be slightly better. The Shepp-Logan phantom in Figure 5(a) is employed. We first computed the parallel Radon transform (1) to then obtain a linear fan-beam sinogram through operation $\mathcal{M}_{\ell}$ defined in (7).

Discrete configuration.

The domain of the feature phantom is the unit disk $\mathinner{\!\left\lVert\bm{x}\right\rVert}_{2}\leq 1$ where we have a $N\times N$ uniformly sampled quadrilateral mesh domain. For the domain of $U$ to compute (1), we have $N\times N_{\theta}$ points uniformly spaced with $(t,\theta)\in[-1,1]\times[0,\pi]$.

For the linear fan-beam geometry, we consider the distance source-origin $D=8$, then $\bar{s}=D(D^{2}-1)^{-1/2}$ is the highest detector position that measures the sample on the unit disk subtending an angle $\bar{\gamma}=\arcsin{1/D}$. We recall that a virtual detector is supposed to be centered, or scaled, at the origin. Due to the symmetry property (6) and to Remark 1, $(s,\beta)$ needs to cover the sets $[-\bar{s},\bar{s}]\times[0,\pi+2\bar{\gamma}[$. For simplicity, we sampled these sets as in the parallel geometry uniformly with $N$ and $N_{\theta}$ points. Finally, here we will use $N=N_{\theta}=512$.

Implementation and complexity.

Implementation details of the two steps of Theorem 4 follow. Step (i) is a sequence of two 1D interpolations, $\mathcal{L}$ followed by $\mathcal{A}$ to obtain $\dot{z}$. $\mathcal{L}$ is a rebinning on the spatial variable from linear to standard geometry, so only needed in the linear case, and $\mathcal{A}$ is a linear operation on the angular variable then not causing strong rebinning errors with linear interpolation. More importantly, both operations are not expensive having complexity $O(NN_{\theta})$ and they were easily parallelized as they are both one-dimensional.

Step (ii) is more challenging, we need to compute a truncated BN series (20) until convergence. Of course, the number of terms in the sequence $\{J_{n}\}$ is the same as that of $\{b_{n}\}$, which itself is the number of samples of $s$ in the sinogram $g$ following the Nyquist sampling theorem [9]. Convergence of the series can thus be reached by first zero-padding the linear variable $s$ of $g$, or equivalently, $\gamma$ of $\dot{z}$ after performing operation $\mathcal{A}$.

Computing coefficients $\{b_{n}\}$ requires $N_{\theta}$ independent 1D Fourier transforms to obtain (5) followed by the easy operation (16). They were also parallelized in $\theta$ having a complexity $O(N\log N)$ each with Fast Fourier techniques. Importantly, the sequence $\{J_{n}(D\sigma)\}$ in (20) not depending on the sample is computed only once for a given setup and can be used as a lookup table for the computation of (20), not being costly at computing time.

After truncating sum (20), computing it results simply in a matrix multiplication of the matrices composed with elements $\{\dot{b}_{n}(\theta)\}$ and $\{J_{n}(D\sigma)\}$ for sampled values of $\{\theta,\sigma\}$. The frequency variable $\sigma$ was sampled following the Nyquist sampling theorem. Such matrix multiplication has a complexity of $O(N^{3})$ but can be reduced to $O(N^{2.3729})$ with fast matrix multiplication algorithms [11, 26], although we used the direct approach in our simulations. We recall that conventional approaches present a cost of $O(N^{3})$. Even using a conventional $O(N^{3})$ matrix multiplication algorithm, the simplicity and straightforward parallel implementation of matrix multiplication related to a standard backprojection algorithm makes this approach still interesting. Finally, the Fourier domain of the backprojected image is obtained at a polar grid followed by linear interpolation to Cartesian coordinates (see section 4.3 for comments on this crucial step) and the 2D inverse (fast) Fourier Transform; operations with respective complexity of $O(N^{2})$ and $O(N^{2}\log N)$.

Reconstruction results.

To obtain image reconstructions, we filtered the parallel sinogram with the conventional ramp filter [9] before applying the rebinning operation $\mathcal{M}_{\ell}$ to obtain a filtered fan-beam sinogram. We proceeded in this way for sake of simplicity because the filtering processes for fan-beam sinograms, detailed in [9, 19], is not in the scope of this work but only the backprojection operation.

Two backprojections methods were implemented, the conventional rebinning through parallel geometry by operation $\mathcal{M}_{\ell}^{-1}$ in (13) followed by the standard $O(N^{3})$ backprojection $\mathcal{B}$ as presented in theorem 1 and our BN series method. The rebinning strategy, used here to evaluate the accuracy of our method, is often employed in available libraries that include fan-beam tomography as in [22]. The normalized mean square error (MSE), computed as $100\mathinner{\!\left\lVert u-\tilde{u}\right\rVert}_{2}^{2}/\mathinner{\!\left\lVert\tilde{u}\right\rVert}_{2}^{2}$, and the normalized mean absolute error (MAE), as $100\mathinner{\!\left\lVert u-\tilde{u}\right\rVert}_{1}/\mathinner{\!\left\lVert\tilde{u}\right\rVert}_{1}$, for a ground truth image $\tilde{u}$ and corrupted image $u$ are used as error metrics.

With noiseless data almost perfect reconstructions were obtained with the rebinning method (image not shown) mainly due to fact that the sinogram was obtained through $\mathcal{M}_{\ell}$. With the BN series the results are also almost exact, shown in Figure 5(b), the associated central $(y=0)$ profile in Figure 5(a) and a particular profile detail in Figure 5(b) where both methods show to reconstruct the same level of target features. Spatial resolution is estimated with the Fourier Ring Correlation (FRC) method with the half-bit criterion as suggested in [23]. Three FRC curves are calculated for the phantom and reconstructions with both methods in Figure 8. Negligible sub-pixel difference is obtained between the three curves especially in the region of the intersection with the threshold curve. Namely, $(1.151,1.153,1.139)$ pixels are respectively obtained as spatial resolution for the phantom, BN series and rebinning methods with a difference of around $0.01$ pixel between phantom and both methods. This shows that no significant resolution is lost with the BN method related to the conventional $O(N^{3})$ rebinning method. Reconstruction MSE and MAE are still better with the rebinning method with noiseless data, as indicated in Figure 7 for the null sinogram MSE case.

To evaluate the noise response of the algorithm, Poisson noise was added to sinograms with increasing sinogram MSE; reconstructed images are shown in Figure 7 with no significant visual difference between both methods. However, as the noise increases, the BN series method behaves better as shown in Figure 7(a) where both reconstruction MSE appears to grow linearly with noise but the BN series with a weaker slope. This is supported with reconstruction MAE behavior in figure 7(b), where the BN series overcomes the rebinning method even earlier (with less noise). This could be due to the effect of 2D interpolations of operator $\mathcal{M}_{\ell}^{-1}$ on noisy sinograms. The polar-to-Cartesian mapping on the Fourier domain is performed both with linear and quadratic splines interpolation on Figure 7(a) as commented in Section 4.3 to show that with the adequate sampling, linear interpolation is enough within polynomial interpolation but can be further improved with advanced interpolation as gridding or regularized strategies.

Both algorithms were implemented on a multithreading architecture for each step. We used the Python library scipy [24] to compute Fast Fourier transforms, interpolations and Bessel coefficients. Besides the theoretical interest of our formula (20), computing time can be reduced with optimized code related to conventional backprojections thanks to the lower computational complexity, when fully exploiting fast matrix multiplication and fast Fourier transform algorithms.

6. Conclusions

We have proposed in this work a low-complexity analytic backprojection formulation for fan-beam tomography with standard and linear detectors to be easily implemented as a matrix multiplication algorithm as well as their associated formulations of the adjoint transform on Hilbert spaces. The linear case being more interesting in a synchrotron context was developed as a two-step process. The first step consists in a sequence of simple one-dimensional operators where a standard fan-beam sinogram is obtained. The second step, the backprojection operation, is performed by writing the Fourier transform of the backprojected image as a Bessel-Neumann series on the frequency variable $\sigma$ weighted by $1/\sigma$. The coefficients of the expansion are in fact the Fourier coefficients of the sinogram obtained in the first step. This approach, based on the low cost backprojection BST formula [16] for parallel sinograms, presents a reduction on the computational cost related to conventional fan-beam backprojections when fast matrix multiplication and fast Fourier transforms algorithms are used, and a simpler formulation as a two-matrix multiplication. Iterative reconstruction algorithms where the backprojection process is computed on each iteration, as the EM algorithm [19] can also take advantage of this formulation to obtain robust reconstructions. A second interesting feature is the absence of a strong rebinning process from fan to parallel beam projections as is mostly done in other fan-backprojection algorithms. This results in better handling of noisy data where interpolations can have strong negative effects as shown in our numerical simulations when linearly increasing the noise level and our method outperforms the conventional rebinning strategy.

Acknowledgements.

PG was supported by the FWO-SBO MetroFleX project (grant agreement S004217N).