New iterative reconstruction methods for fan-beam tomography

Abstract

In this paper, we present a novel class of iterative reconstruction methods for severely angular undersampled or/and limited-view tomographic problems with fan-beam scanning geometry. The proposed algorithms are based on a new analytical transform which generalizes Fourier-slice theorem to divergent-beam scanning geometries. Using a non-rigid coordinate transform, divergent rays can be reorganized into parallel ones. Therefore, one can employ a simpler parallel-beam projection model instead of more complicated divergent-beam geometries. Various existing iterative reconstruction techniques for divergent-beam geometries can be easily adapted to the proposed framework. The significant advantage of this formulation is the possibility of exploiting efficient Fourier-based recovery methods without rebinning of the projections. In case of highly sparse measurements (few-view data), rebinning methods are not suitable due to error-prone angular interpolation involved. In this work, three new methods based on the novel analytical framework for fan-beam geometry are presented: the Gerchberg-Papoulis algorithm, the Neumann decomposition method and its total variation regularized version. Presented numerical experiments demonstrate that the methods can be competitive in reconstructing from few-view noisy tomographic measurements.

Keywords:

• Radon transform
• Fourier methods
• few-view
• tomography
• inverse problems
• integral equations
• numerical solutions
• computational methods

AMS Subject Classifications:

• 65N21
• 44A12

1. Introduction

Direct reconstruction methods in X-ray computed tomography (CT) are proven to be fast and simple to implement [1,2]. They are based on a one-step approximation of a solution to inverse problem by backprojecting all measured tomographic projections. In order to overcome the blurring effect of the back-summation procedure, projections are convolved with a high-pass filter (e.g. a ramp filter). In practice, commonly used direct methods include Filtered Back Projection (FBP) algorithm for parallel and fan-beam geometries and Feldkamp–Davis–Kress (FDK) algorithm for the cone-beam geometry [3]. The filtered backprojection procedure can be considered as a discrete approximation to inverse Radon transform formula which maps the space of measurements to the space of objects.

The FBP algorithm for parallel-beam geometry can also be introduced using Fourier-slice theorem (FST) [1,3,4]. In two-dimensional (2D) case, FST states that the one-dimensional (1D) Fourier Transform (FT) of a single projection gives the corresponding orthogonal ’slice’ in 2D spectral domain. Therefore by taking 1D FT of an each subsequent projection, one can fill the 2D frequency spectrum with the corresponding Fourier coefficients in polar coordinates. In order to reconstruct an image, interpolation from polar to Cartesian coordinates (a regridding procedure) and the 2D inverse FT are required. When the measured projection space is sampled richly and in accordance with the Nyquist–Shannon theorem, the reconstruction with direct Fourier method tends to be similar to the convolution-backprojection approach. However, when the projection data are sparse and undersampled, the Fourier-based reconstruction can result in aliasing artefacts. Normally, the error originates from the regridding step which can lead to frequencies overlapping due to angular undersampled measurements. This phenomenon is inherent for higher frequencies, which are responsible for the edges in images. In practice, the convolution-based reconstruction methods are used due to straightforward implementation [3,4]. However, the Fourier-based reconstruction approaches are computationally efficient due to Fast Fourier Transform (FFT) [3].

In general, however, all direct methods fail to reconstruct from noisy sparse measurements due to the ill-posed and ill-conditioned nature of the inverse problem [5]. When the sampling rate is lower than it is required by the Nyquist–Shannon theorem (twice the frequency maximum present in data), a high level of noise and the loss of spatial resolution are expected in the reconstructed images. In clinical CT, a limited sampling rate is associated with the requirements to minimize the ionizing radiation dose. In material science, one would like to improve temporal resolution of dynamic (4D) experiments by performing faster acquisition with shorter exposures [6].

Iterative Image Reconstruction (IIR) techniques are better adapted to deal with ill-posed inversion than direct methods as they can account for more accurate modelling of the acquisition process and allow the use of a priori information [5]. Due to inaccurate and noisy measurements, IIR methods often require an additional regularity on the solution to ensure well posedness of inversion. For instance, the spatial regularization can be imposed through smoothness (Tikhonov ${\ell }_2$-norm or Laplacian) [5], gradient sparsity (Total Variation (TV) semi-norm) [7], piecewise-smooth ${\ell }_{1-2}$ [8], various time-frequency-transformations [8], patch-based priors [6] and others.

The majority of IIR techniques are based on the projection-image space alternating error-correcting refinements where a priori information is applied in the image space only. The class of IIR methods based on the FT has an additional flexibility of embedding prior information in both spatial and frequency domains. Some examples of those methods are non-uniform FFT (NUFFT) [9,10], equally sloped tomography method [11,12], and Gerchberg-Papoulis (GP) algorithm [2,13–17]. In this paper, we will focus on the GP method – a powerful and efficient FT-based IIR technique which is well suited for the reconstruction from sparse tomographic measurements. The GP method was initially introduced to solve 1D problem of out-of-band signal extrapolation [18,19]. Gerchberg proved convergence of the method in the absence of noise [18] and the regularized version was proposed later in [14]. The significant advantage of the GP algorithm over the main class of projection-image space IIR methods is the ability of suppressing noise and image artefacts by applying various constraints in both spectral and image domains. In the spatial domain, conventional regularization techniques can be used: non-negativity, smoothness, defined spatial support, bounded intensity values and others. In the spectral domain, frequency filtering, band-limiting (e.g. zeroing coefficients outside known frequency support), adaptive to sparse measurements k-space interpolation and other techniques can be applied. The use of the GP method for super-resolution problems [20,21] provides an insight to similar challenges in tomographic reconstruction, such as the limited data and missing wedge reconstruction [8]. The GP method enables computationally efficient and high-quality reconstructions from few-view tomographic data. The disadvantage of the GP method is its inapplicability to divergent-beam measurements due to the limitations of the classical FST.

Previously, there were different attempts to generalize the classical FST to divergent-beam geometries [22–24]. These techniques were developed to avoid rebinning of divergent rays into parallel rays [4]. The rebinning step can be expensive memory-wise for larger data-sets since it requires storing the result of interpolation. Notably, when only a few projections available and the reconstruction problem is severely angular undersampled, the rebinning step will result in significant errors due to interpolation. Therefore, the rebinning of projections should be avoided in this case.

In [22], a novel approach to generalization of the FST was proposed which is different from the methods in [23,24]. The core of a method is in the simple non-rigid coordinate transform in image space which allows the further use of the classical FST for divergent-beam geometries [22]. This formulation relaxes the difficulty which is associated with the use of the FST for divergent-beam geometries and therefore any FT-based IIR method can be implemented using the properties of the Generalized FST (GFST).

In this paper, along with the Fan-Beam GP (FBGP) method [16], we present two new iterative algorithms which are also based on the GFST: Neumann Decomposition Shepp-Logan (NDSL) method and its TV-regularized version (NDSL-TV). The NDSL method was successfully used previously for the parallel-beam geometry [25,26] and it is based on the decomposition of inverse projection operator into the Neumann series of approximate operators [27,28]. Here we reformulated the NDSL method to fan-beam geometry using the GFST and presented the TV-regularized version of the method (NDSL-TV). To make our presentation more self-consistent, we introduce the GFST in Section 2.1. The FBGP method is presented in Section 2.2 and details to improve its performance given in Section 2.3. The NDSL and NDSL-TV methods are presented in Section 2.4. The numerical experiments for comparison of the FBGP, NDSL and NDSL-TV methods are given in Section 3. Discussion is given in Section 4 and conclusions are given in Section 5.

2. Theory

2.1. The GFST

Here we introduce the GFST [22] in application to fan-beam geometry with a flat detector; however, the theorem can be extended to other divergent-beam geometries (e.g. the cone-beam geometry). The core of the GFST lies in a non-linear coordinate transform which helps to translate the fan-beam problem into the parallel-beam one. For each angular scanning position, the object undergoes a unique forward coordinate deformation where divergent rays become parallel and therefore the classical FST can be used directly. In order to return to the original coordinate system, the inverse deformation must be applied.

Here we consider the standard acquisition scheme for the 2D fan-beam geometry in CT (see Figure 1). Let the domain of unknown function g(x, y) is bounded by the unit circle $x^2+y^2\le 1$, where fan-beam source E rotates along the circular path and (x, y) is the coordinate system rotated by the angle $\mathit{\beta }$, such as$$x(\mathit{\beta })=-Dsin\mathit{\beta },y(\mathit{\beta })=Dcos\mathit{\beta };$$

where $D=OE$ is a distance from E to the origin of the rotated coordinate system (s, p). The detectors are located on the straight line $D_1{D}_2$. For the fixed angle $\mathit{\beta }$, each ray EA inside the fan-beam can be characterized by its angle $\mathit{\gamma }$ from the central ray, and its distance q from O. The angle between the x-axis and the beam is denoted by $\mathit{\xi }$, while $\mathit{\theta }$ denotes an angle between the x-axis and the line perpendicular to the beam. The measured signals on the detectors line $D_1{D}_2$ can also be considered on the line of virtual detectors $D_1^{{}^{\prime }}{D}_2^{{}^{\prime }}$.

Figure 1. Fan-beam scanning geometry with a flat detector.

We propose the following projective (’deforming’) non-linear transformation of the point (s, p) in the rotated coordinate system:(1) $$\begin{array}{c}\hfill \left\{\begin{array}{c}u=s/(1-\frac{p}{D})=s/Q;\hfill \\ v=p.\hfill \end{array}\right.\end{array}$$(1)

The transform (1(1) $$\begin{array}{c}\hfill \left\{\begin{array}{c}u=s/(1-\frac{p}{D})=s/Q;\hfill \\ v=p.\hfill \end{array}\right.\end{array}$$(1) ) changes rays in the fan-beam in (x, y) coordinate system into straight lines parallel to the v-axis. Therefore, for the corresponding angle $\mathit{\beta }$, the initial image g(x, y) will be deformed into a new image g(X(u, v), Y(u, v)):(2) $$\begin{array}{c}\hfill \left\{\begin{array}{c}x=X(u,v)=u(1-\frac{v}{D})cos\mathit{\beta }-vsin\mathit{\beta };\hfill \\ y=Y(u,v)=u(1-\frac{v}{D})sin\mathit{\beta }+vcos\mathit{\beta }.\hfill \end{array}\right.\end{array}$$(2)

The inverse substitution in (2(2) $$\begin{array}{c}\hfill \left\{\begin{array}{c}x=X(u,v)=u(1-\frac{v}{D})cos\mathit{\beta }-vsin\mathit{\beta };\hfill \\ y=Y(u,v)=u(1-\frac{v}{D})sin\mathit{\beta }+vcos\mathit{\beta }.\hfill \end{array}\right.\end{array}$$(2) ): $u=U(x,y),v=V(x,y)$, where $(r,\mathit{\varphi })$ are polar coordinates of (x, y) and $(q,\mathit{\theta })$ are the coordinates of a fan-beam passing through this point. Notably, the measured signal becomes zero if $|s|\ge D/\sqrt{D^2-1}$. Then fan-beam projection can be expressed as follows:(3) $$\begin{array}{cc}\hfill f_{\mathit{\beta }}(u)& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathit{\delta }(q-U(x,y))g(x,y)\mathrm{d}x\mathrm{d}y\hfill \\ \hfill & ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\left|J\right|\mathit{\delta }(q-U)g(u,v)\mathrm{d}u\mathrm{d}v.\hfill \end{array}$$(3)

The transition from the coordinates (x, y) to (u, v) is carried out with the help of the Jacobian matrix J, which is defined using (2(2) $$\begin{array}{c}\hfill \left\{\begin{array}{c}x=X(u,v)=u(1-\frac{v}{D})cos\mathit{\beta }-vsin\mathit{\beta };\hfill \\ y=Y(u,v)=u(1-\frac{v}{D})sin\mathit{\beta }+vcos\mathit{\beta }.\hfill \end{array}\right.\end{array}$$(2) ) as:(4) $$\begin{array}{cc}\hfill J& =\left[\begin{array}{cc}\mathit{\partial }x/\mathit{\partial }u& \mathit{\partial }x/\mathit{\partial }v\\ \mathit{\partial }y/\mathit{\partial }u& \mathit{\partial }y/\mathit{\partial }v\end{array}\right]\hfill \\ \hfill & =\left(1-\frac{v}{D}\right)=Q.\hfill \end{array}$$(4)

Finally, the Equation (3(3) $$\begin{array}{cc}\hfill f_{\mathit{\beta }}(u)& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathit{\delta }(q-U(x,y))g(x,y)\mathrm{d}x\mathrm{d}y\hfill \\ \hfill & ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\left|J\right|\mathit{\delta }(q-U)g(u,v)\mathrm{d}u\mathrm{d}v.\hfill \end{array}$$(3) ) in a new system of coordinates can be written as:(5) $$\begin{array}{cc}\hfill f_{\mathit{\beta }}(u)& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }|Q|\mathit{\delta }(q-U)g(u,v)dudv\hfill \\ \hfill & =\frac{1}{cos\mathit{\gamma }(u)}{\int }_{-\infty }^{\infty }g(X(u,v),Y(u,v))\mathrm{d}v.\hfill \end{array}$$(5)

In the derivation of Equation (5(5) $$\begin{array}{cc}\hfill f_{\mathit{\beta }}(u)& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }|Q|\mathit{\delta }(q-U)g(u,v)dudv\hfill \\ \hfill & =\frac{1}{cos\mathit{\gamma }(u)}{\int }_{-\infty }^{\infty }g(X(u,v),Y(u,v))\mathrm{d}v.\hfill \end{array}$$(5) ), we used the formulae:(6) $$\begin{array}{cc}\hfill \mathit{\delta }(q-U)& =\mathit{\delta }(q-rcos(\mathit{\theta }-\mathit{\varphi }))\hfill \\ \hfill & =\mathit{\delta }(Dsin\mathit{\gamma }-cos\mathit{\gamma }[rcos(\mathit{\beta }-\mathit{\varphi })-tan\mathit{\gamma }rsin(\mathit{\beta }-\mathit{\varphi })])\hfill \\ \hfill & =\frac{1}{|cos\mathit{\gamma }|}\mathit{\delta }(u(1+\frac{rsin(\mathit{\beta }-\mathit{\varphi })}{D})-rcos(\mathit{\beta }-\mathit{\varphi }))\hfill \\ \hfill & =\frac{1}{|cos\mathit{\gamma }|}\mathit{\delta }(uQ-rcos(\mathit{\beta }-\mathit{\varphi }))=\frac{1}{\left|Qcos\mathit{\gamma }(u)\right|}\mathit{\delta }(u-u^{\prime }),\hfill \end{array}$$(6)

where $u^{\prime }=\frac{rcos(\mathit{\beta }-\mathit{\varphi })}{Q}=s/Q$. From the geometrical formulation of the problem (see Figure 1) $\mathit{\gamma }<\mathit{\pi }/2$, which implies that $|cos\mathit{\gamma }|=cos\mathit{\gamma }$ (notice that in fan-beam geometry $Q>0$, $Q=|Q|$). As a result, the Jacobian module $\left|Q\right|$ can be cancelled with the multiplier ${|Q|}^{-1}$ and removed from the delta-function.

Furthermore, we rename the integrand g(u, v) in the integral (5(5) $$\begin{array}{cc}\hfill f_{\mathit{\beta }}(u)& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }|Q|\mathit{\delta }(q-U)g(u,v)dudv\hfill \\ \hfill & =\frac{1}{cos\mathit{\gamma }(u)}{\int }_{-\infty }^{\infty }g(X(u,v),Y(u,v))\mathrm{d}v.\hfill \end{array}$$(5) ) into $h_{\mathit{\beta }}(u,v)$ and define a new function $f_{\mathit{\beta }}^{\prime }(u)=f_{\mathit{\beta }}(u)cos\mathit{\gamma }(u)$. Here the index $\mathit{\beta }$ appears since the (u, v) coordinates depend on the angle $\mathit{\beta }$. The FT of the function $f_{\mathit{\beta }}^{\prime }(u)$ leads to the generalized formulation of the FST for fan-beam geometry [22]:(7) $$\begin{array}{cc}\hfill {\widetilde{f^{\prime }}}_{\mathit{\beta }}({\mathit{\nu }}_u)& ={\int }_{-\infty }^{\infty }\left[{\int }_{-\infty }^{\infty }{h}_{\mathit{\beta }}(u,v)\mathrm{d}v\right]e^{-2\mathit{\pi }i{\mathit{\nu }}_{u}u}\mathrm{d}u\hfill \\ \hfill & ={\tilde{h}}_{\mathit{\beta }}({\mathit{\nu }}_u,{\mathit{\nu }}_v){|}_{{\mathit{\nu }}_v=0}.\hfill \end{array}$$(7)

Thus, the GFST is now expressed in (u, v) coordinates which are obtained from the initial coordinates (s, p) by the non-linear projective variable transformation (1(1) $$\begin{array}{c}\hfill \left\{\begin{array}{c}u=s/(1-\frac{p}{D})=s/Q;\hfill \\ v=p.\hfill \end{array}\right.\end{array}$$(1) ). This is a one-to-one transform between the domain of the unknown function g(x, y) and the deformed function $h_{\mathit{\beta }}(u,v)$. Hence, the FT of the deformed fan-beam projection (multiplied by the cosine of the incidence angle of each ray in the fan-beam) coincides with the central section of the 2D Fourier-image of the deformed object $h_{\mathit{\beta }}(u,v)$.

In Figure 2, we illustrate the deformation procedure (1(1) $$\begin{array}{c}\hfill \left\{\begin{array}{c}u=s/(1-\frac{p}{D})=s/Q;\hfill \\ v=p.\hfill \end{array}\right.\end{array}$$(1) ) applied to a simple 2D phantom which consists of two circles and a rectangle. Spatial transformation from (s, p) coordinates into the deformed (u, v) coordinates is performed using bilinear interpolation. In Figure 2(b), direct deformation of the phantom (a) is given and (c) shows the inverse deformation of (b). Notably, when the inverse transform is applied to the deformed image (b) the object returns to its initial geometrical shape (c). In Figure 2(c), one can see that the parallel rays in the deformed space (b) have become fan-rays in the original coordinate system. Therefore, one can scan an object using fan-beam geometry and after transformation (1(1) $$\begin{array}{c}\hfill \left\{\begin{array}{c}u=s/(1-\frac{p}{D})=s/Q;\hfill \\ v=p.\hfill \end{array}\right.\end{array}$$(1) ) the resulting deformed object (b) is suitable for parallel-beam geometry. Furthermore, in image (c) one can notice that parallel and fan-rays cross each other exactly at the same points along the x-axis, which proves that the deformation has been applied correctly. Therefore in the deformed space, the classical FST can be applied directly.

Figure 2. (a) original phantom, $N_x=N_y=256$; (b) deformed phantom, $\mathit{\beta }={180}^{\circ },D=1.5$; (c) inverse deformation of (b).

Figure 3. (a) original image of Lena (cropped in a circle shape), $N_x=N_y=256$; (b) deformed image, $\mathit{\beta }={180}^{\circ },D=1.4$; (c) inverse deformation of (b). The minor isotropic blurring effect of two bilinear interpolations (forward-backward deformation steps) is noticeable in (c).

For the sake of interesting morphing effects of non-linear deformation (1(1) $$\begin{array}{c}\hfill \left\{\begin{array}{c}u=s/(1-\frac{p}{D})=s/Q;\hfill \\ v=p.\hfill \end{array}\right.\end{array}$$(1) ), we repeat the experiment with an image of Lena (see Figure 3). Although image (c) looks similar to the original image (a), it has been deformed twice by means of bilinear interpolation and it is actually a smoothed version of (a). The repetitive (forward-backward) interpolation steps of GFST-based IIR algorithm can result in over-smoothing of the solution. More accurate interpolation approaches can lead to better results, therefore in our experiments we use bicubic interpolation [29]. In our implementation, we use Matlab function interp2 with ’cubic’ parameter [30].

2.2. Gerchberg-Papoulis algorithm for fan-beam geometry

The classical iterative GP method for parallel-beam scanning geometry is based on the application of the FST and it is presented in Algorithm 1 [2,13–17]. In Algorithm 1, $F_2$ and $F_2^{-1}$ are the operators of direct and inverse 2D Fourier-transform, respectively. The operators ${\mathrm{\Phi }}_s^{(n)}$ and ${\mathrm{\Phi }}_f^{(n)}$ are responsible for applying a priori information (every n-th iteration) in spatial and frequency domains, respectively (see Table 1). Additionally, the operator ${\mathrm{\Phi }}_f$ maps Fourier spectrum of projections on the radial lines to the 2D Cartesian Fourier space. The regridding procedure is an important part of the GP algorithm for reconstruction from undersampled projection data and it requires a special care to avoid interpolation errors [16] (see also Section 2.3). Before giving more technical details on operators ${\mathrm{\Phi }}_{s,f}$, we introduce the GP method for fan-beam geometry (FBGP) [16].

Table 1. The list of actions inside the operators ${\mathrm{\Phi }}_f$ and ${\mathrm{\Phi }}_s$.

${\mathrm{\Phi }}_f$	${\mathrm{\Phi }}_s$
• Strip interpolation on the radial lines (see Section 2.3)	• Non-negativity
• Zeroing of frequencies outside the strip	• Bounded support
• Filtering (11(11) $$\begin{array}{c}\hfill {\tilde{g}}^{(n)}({\mathit{\nu }}_x,{\mathit{\nu }}_y)\to \frac{{\tilde{g}}^{(n)}({\mathit{\nu }}_x,{\mathit{\nu }}_y)}{1+\mathit{\alpha }{\mathit{\nu }}^2},\end{array}$$(11) )

Let us define a deformation matrix T as:(8) $$\begin{array}{c}\hfill T=\left[\begin{array}{cc}{Q}^{-1}& 0\\ 0& 1\end{array}\right],Q=(1-\frac{p}{D}).\end{array}$$(8)

Using matrix T, one can transform the coordinate system (s, p) into a new coordinate system (u, v): $\left(\genfrac{}{}{0pt}{}{u}{v}\right)=T\left(\genfrac{}{}{0pt}{}{s}{p}\right)$. The rotation of the coordinate system (x, y) by $\mathit{\beta }$ (see Figure 1) is defined by the following rotation matrix:(9) $$\begin{array}{c}\hfill R_{\mathit{\beta }}=\left[\begin{array}{cc}cos\mathit{\beta }& sin\mathit{\beta }\\ -sin\mathit{\beta }& cos\mathit{\beta }\end{array}\right].\end{array}$$(9)

Then the space of the deformed object can be represented using new coordinates $(x^{\prime },y^{\prime })$ as:(10) $$\begin{array}{c}\hfill \left(\genfrac{}{}{0pt}{}{x^{\prime }}{y^{\prime }}\right)=R_{\mathit{\beta }}T{R}_{-\mathit{\beta }}\left(\genfrac{}{}{0pt}{}{x}{y}\right),\left(\genfrac{}{}{0pt}{}{x}{y}\right)=R_{-\mathit{\beta }}{T}^{-1}{R}_{\mathit{\beta }}\left(\genfrac{}{}{0pt}{}{x^{\prime }}{y^{\prime }}\right),\end{array}$$(10)

where $R_{-\mathit{\beta }}$ is the inverse of $R_{\mathit{\beta }}$. Note that one can rotate an object first and then deform it (two interpolations involved) or the deformation can be performed directly from the current source position (one interpolation). To reduce computational cost and minimize interpolation errors, we use the later case with bicubic interpolation instead of bilinear. Therefore, only two 2D interpolations are needed for forward and inverse deformations (10(10) $$\begin{array}{c}\hfill \left(\genfrac{}{}{0pt}{}{x^{\prime }}{y^{\prime }}\right)=R_{\mathit{\beta }}T{R}_{-\mathit{\beta }}\left(\genfrac{}{}{0pt}{}{x}{y}\right),\left(\genfrac{}{}{0pt}{}{x}{y}\right)=R_{-\mathit{\beta }}{T}^{-1}{R}_{\mathit{\beta }}\left(\genfrac{}{}{0pt}{}{x^{\prime }}{y^{\prime }}\right),\end{array}$$(10) ). Apart from the forward-inverse deformations for angle $\mathit{\beta }$, the FBGP method (see Algorithm 2) is exactly the same as the classical GP method (see Algorithm 1).

2.3. Some important aspects of the GP and the FBGP algorithms

The 1D FT can be applied to a measured projection which results in a set of discrete points lying along the radial line in the frequency domain.

Figure 4. (a) The set of radial lines for four projections in polar coordinates with Cartesian coordinate system superimposed. The symmetric strip of $2S_w$ width is centred on each radial line determining the trust region where interpolation is performed. All frequency coefficients outside the strip are considered unreliable and therefore ignored during the interpolation (zeroed before the first 2D IFT). The strip interpolation is proven to be much more accurate than the classical isotropic interpolation methods for undersampled projection datasets; (b) Zoomed strip region shows two 1D interpolations employed within the strip to regrid radial data points to Cartesian grid.

Therefore, for every projection angle, the spectral domain is filled with data-contained radial lines (see Figure 4(a)). In order to apply the 2D inverse FT (IFT) to $\tilde{g}({\mathit{\nu }}_r,{\mathit{\nu }}_{\mathit{\varphi }})$, it is necessary to interpolate values lying on the radial lines in polar coordinate system to Cartesian coordinate system $({\mathit{\nu }}_x,{\mathit{\nu }}_y)$. Traditional spatially invariant regridding techniques are extremely sensitive to outliers [3,10,11]. In the case of limited projections, the spectral domain is sampled non-uniformly with bigger gaps (the missing wedges) further from the centre of coordinates (see Figure 4(a)). However, due to noise present in the data, the missing wedges contain some frequencies which lead to image artefacts after application of the 2D IFT. Therefore, one should consider applying interpolation within some trust region only where data are more reliable. To increase accuracy, we employ a strip interpolation method [16] which is proven to be successful to recover images with less artefacts when the data spectrum is sparse.

In our implementation, we choose a symmetrical strip with defined width $2S_w$, which is centred on each radial line in polar coordinate system (see Figure 4). Inside the strip, two 1D interpolations are performed, the nearest neighbour interpolation across the strip and linear interpolation along the strip (see Figure 4(b)). All data points within the strip are considered to have equal weights; however, additional weighting factors can be assigned to more distant from the central radial line data points.

To further improve the accuracy of the strip interpolation, the width of a strip continuously shrinks in iterations with a period $m_0$ and some multiplier $0<{\mathit{\delta }}_s\le 1$. At each n-iteration, the condition $mod(n,m_0)=0$ is checked and if fulfilled the width of a strip is updated as $S_w=S_{n-1}·{\mathit{\delta }}_s$. Therefore if $m_0=1$ we update the strip every iteration. The initial width $2S_w$ and the multiplier ${\mathit{\delta }}_s$ are important parameters which can affect the convergence of the GP algorithm. The choice of $S_w$ is usually an empirical one and also can depend on the sampling rate. The main idea of the iterative change of the width is in the ability of the algorithm to adapt to newer estimates by improving the accuracy of interpolation (this is equivalent to the refinement of the trust region based on new data). Normally, a wide strip leads to a slower convergence and too narrow strip can result in large reconstruction errors (or even divergence of the algorithm). Previously, it was established [16] that it is important to iteratively reduce the strip’s width in iterations until a small width value is reached. Alternatively, if the width of a strip is kept constant in iterations, then the algorithm can stagnate on earlier iterations and even diverge.

Apart from the non-standard interpolation procedure (included into operator ${\mathrm{\Phi }}_f^{(n)}$, see Table 1), there are additional elements which can help to obtain better reconstruction and accelerate convergence of the GP method. For instance, the measured spectrum can be bounded using the Nyquist theorem in order to remove some spectral coefficients, e.g. to zero frequencies above the Nyquist frequency. It is also important to regularize reconstruction and it can be done conveniently in the Fourier domain [14]. High frequencies can be suppressed in accordance with the level of noise in projections:(11) $$\begin{array}{c}\hfill {\tilde{g}}^{(n)}({\mathit{\nu }}_x,{\mathit{\nu }}_y)\to \frac{{\tilde{g}}^{(n)}({\mathit{\nu }}_x,{\mathit{\nu }}_y)}{1+\mathit{\alpha }{\mathit{\nu }}^2},\end{array}$$(11)

where $\mathit{\alpha }$ is a regularization parameter and ${\mathit{\nu }}^2={\mathit{\nu }}_x^2+{\mathit{\nu }}_y^2$. The regularization parameter can be found empirically or using a discrepancy principle [31]. In the image domain (operator ${\mathrm{\Phi }}_s^{(n)}$), one can apply positivity and bounded support constraints. Notably one can also regularize in image space, however we do not implement this in our modification of the GP algorithm.

In order to quantify reconstructed images, we use the following relative error:(12) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_1=\sqrt{\frac{\Vert g-\widehat{g}{\Vert }_2}{\Vert \widehat{g}{\Vert }_2}}·100\%,\end{array}$$(12)

where g is a reconstructed image and $\widehat{g}$ is an exact phantom.

2.4. Iterative method of Neumann decomposition with Shepp-Logan filter for fan-beam geometry

In this section, we present a different reconstruction approach which is also based on the proposed deformation (8(8) $$\begin{array}{c}\hfill T=\left[\begin{array}{cc}{Q}^{-1}& 0\\ 0& 1\end{array}\right],Q=(1-\frac{p}{D}).\end{array}$$(8) ).

The tomographic reconstruction problem can be formulated as a solution to the system of linear algebraic equations [1–5]:(13) $$\begin{array}{c}\hfill \mathit{f}=\mathrm{A}\mathit{g}+\mathit{\delta },\end{array}$$(13)

where $\mathit{g}\in {\mathbb{R}}^N$ is a vectorized representation of unknown object (N is the total number of image elements), $\mathit{f}\in {\mathbb{R}}^M$ is a vectorized sinogram, $\mathit{\delta }\in {\mathbb{R}}^M$ is a noise component present in data and $\mathrm{A}\in {\mathbb{R}}^{M\times N}$ is a sparse system projection matrix which maps the ’space of objects’ into the ’space of observations’. The formal solution to the problem (13(13) $$\begin{array}{c}\hfill \mathit{f}=\mathrm{A}\mathit{g}+\mathit{\delta },\end{array}$$(13) ) can be written as $\mathit{g}={\mathrm{A}}^{-1}\mathit{f}$, where ${\mathrm{A}}^{-1}$ is the inverse or pseudo-inverse of $\mathrm{A}$. However, since $\mathrm{A}$ is strongly ill-conditioned for undersampled data the direct inversion is not feasible and also not unique in presence of noise $\mathit{\delta }$ [5]. The approximate solution can be found through the Moore–Penrose pseudo-inversion or iteratively using the Least-Squares (LS) criterion where the ${\ell }_2$-norm is minimized by residual ${\Vert \mathit{f}-\mathrm{A}\mathit{g}\Vert }_2^2$. Alternatively, one can decompose ${\mathrm{A}}^{-1}$ matrix into a Neumann series [27] and arrive at the iterative scheme which can be adapted to solve the reconstruction problem (13(13) $$\begin{array}{c}\hfill \mathit{f}=\mathrm{A}\mathit{g}+\mathit{\delta },\end{array}$$(13) ) [25,26,28,32].

The projection matrix $\mathrm{A}$ can be decomposed as $\mathrm{A}=\mathrm{B}{A}_0,\mathrm{B}=\mathrm{A}{A}_0^{-1}$, where ${\mathrm{A}}_0^{-1}$ is some known approximation to the unknown operator ${\mathrm{A}}^{-1}$. Using the Neumann series decomposition [27] one can re-write operator ${\mathrm{B}}^{-1}$ as(14) $$\begin{array}{c}\hfill \mathit{g}={\mathrm{A}}_0^{-1}\sum _{k=0}^{\infty }{\mathit{\tau }}^k{(\mathrm{I}-\mathrm{B})}^k\mathit{f}={\mathrm{A}}_0^{-1}{\mathrm{B}}^{-1}\mathit{f}={\mathrm{A}}^{-1}\mathit{f},\end{array}$$(14)

where $\mathrm{I}$ is the identity matrix and $\mathit{\tau }$ is a positive relaxation parameter. Convergence of the series (14(14) $$\begin{array}{c}\hfill \mathit{g}={\mathrm{A}}_0^{-1}\sum _{k=0}^{\infty }{\mathit{\tau }}^k{(\mathrm{I}-\mathrm{B})}^k\mathit{f}={\mathrm{A}}_0^{-1}{\mathrm{B}}^{-1}\mathit{f}={\mathrm{A}}^{-1}\mathit{f},\end{array}$$(14) ) is ensured by $\Vert \mathrm{A}\Vert ={sup}_{\mathit{f}}\frac{\Vert \mathrm{A}\mathit{f}\Vert }{\Vert \mathit{f}\Vert }<1$.

For tomographic reconstruction problem, we can re-write (14(14) $$\begin{array}{c}\hfill \mathit{g}={\mathrm{A}}_0^{-1}\sum _{k=0}^{\infty }{\mathit{\tau }}^k{(\mathrm{I}-\mathrm{B})}^k\mathit{f}={\mathrm{A}}_0^{-1}{\mathrm{B}}^{-1}\mathit{f}={\mathrm{A}}^{-1}\mathit{f},\end{array}$$(14) ) as [25,28]:(15) $$\begin{array}{c}\hfill {\mathbf{g}}^{n+1}={\mathbf{g}}^n+\mathit{\tau }{\mathrm{A}}_0^{-1}\mathrm{\Delta }{\mathbf{r}}^n,\end{array}$$(15)

where $\mathrm{\Delta }{\mathit{r}}^n=f-\mathrm{A}{\mathit{g}}^n$ is a discrepancy vector and ${\mathrm{A}}_0^{-1}$ is replaced with the known approximation of the continuous inverse Radon transform ${\mathrm{A}}_0^{-1}=:{\mathcal{R}}_0^{-1}$ given as [1]:(16) $$\begin{array}{c}\hfill g(x,y)={\mathcal{R}}_0^{-1}f(\mathit{\xi },p)=-\frac{1}{2{\mathit{\pi }}^2}{\int }_0^{\mathit{\pi }}\frac{\mathrm{d}\mathit{\xi }}{Q^2}{\int }_{-a}^{a}f(\mathit{\xi },p)\mathit{\psi }(p-p_0)dp,\end{array}$$(16)

where $p_0=-xsin\mathit{\xi }+ycos\mathit{\xi }$, $a=(O{D}_1^{{}^{\prime }})$ or $2a=(D_1^{{}^{\prime }}{D}_2^{{}^{\prime }})$ (see Figure 1), and $\mathit{\psi }(p)$ is a Shepp-Logan (SL) filter:(17) $$\begin{array}{c}\hfill \mathit{\psi }(k{h}_p)=\frac{4}{h_p^2(4k^2-1)},k=0,\pm 1,\pm 2,\dots ,\end{array}$$(17)

where $h_p$ is a step on p and Q is defined by formula (8(8) $$\begin{array}{c}\hfill T=\left[\begin{array}{cc}{Q}^{-1}& 0\\ 0& 1\end{array}\right],Q=(1-\frac{p}{D}).\end{array}$$(8) ). Note that the element to deal with fan-beam projections using the GFST is embedded into Q operator in (16(16) $$\begin{array}{c}\hfill g(x,y)={\mathcal{R}}_0^{-1}f(\mathit{\xi },p)=-\frac{1}{2{\mathit{\pi }}^2}{\int }_0^{\mathit{\pi }}\frac{\mathrm{d}\mathit{\xi }}{Q^2}{\int }_{-a}^{a}f(\mathit{\xi },p)\mathit{\psi }(p-p_0)dp,\end{array}$$(16) ).

The classical FBP algorithm with the SL filter in the inner integral is expressed in (16(16) $$\begin{array}{c}\hfill g(x,y)={\mathcal{R}}_0^{-1}f(\mathit{\xi },p)=-\frac{1}{2{\mathit{\pi }}^2}{\int }_0^{\mathit{\pi }}\frac{\mathrm{d}\mathit{\xi }}{Q^2}{\int }_{-a}^{a}f(\mathit{\xi },p)\mathit{\psi }(p-p_0)dp,\end{array}$$(16) ) and iterative method (15(15) $$\begin{array}{c}\hfill {\mathbf{g}}^{n+1}={\mathbf{g}}^n+\mathit{\tau }{\mathrm{A}}_0^{-1}\mathrm{\Delta }{\mathbf{r}}^n,\end{array}$$(15) ) is defined as Neumann Decomposition Shepp-Logan (NDSL) algorithm. Notably, the NDSL method shares similarities with the Landweber reconstruction method [5], which is a well-known iterative technique to solve the LS problems. The main difference between the two methods is the nature of the operator ${\mathrm{A}}_0^{-1}$ (15(15) $$\begin{array}{c}\hfill {\mathbf{g}}^{n+1}={\mathbf{g}}^n+\mathit{\tau }{\mathrm{A}}_0^{-1}\mathrm{\Delta }{\mathbf{r}}^n,\end{array}$$(15) ). The Landweber algorithm uses ${\mathrm{A}}_0^{-1}:={\mathrm{A}}^T$, where ${\mathrm{A}}^T$ is a backprojection (adjoint) operator without filtering. However, one can consider the NDSL algorithm (15(15) $$\begin{array}{c}\hfill {\mathbf{g}}^{n+1}={\mathbf{g}}^n+\mathit{\tau }{\mathrm{A}}_0^{-1}\mathrm{\Delta }{\mathbf{r}}^n,\end{array}$$(15) ) as an accelerated version of the Landweber method [33], when high-pass filtering of the sinogram’s residual is performed before the backprojection.

Additionally, to ensure stable convergence of the NDSL method (15(15) $$\begin{array}{c}\hfill {\mathbf{g}}^{n+1}={\mathbf{g}}^n+\mathit{\tau }{\mathrm{A}}_0^{-1}\mathrm{\Delta }{\mathbf{r}}^n,\end{array}$$(15) ) we regularize it using the TV penalty [7]. The reconstruction-regularization alternating process is based on the operator splitting method [12,13] and summarized in Algorithm 3. The TV minimization problem in Algorithm 3 is solved using explicit Rudin–Osher–Fatemi (ROF) scheme [7] with a small time step parameter to ensure convergence. In this experiment, the regularization parameter $\mathit{\lambda }$ has been found empirically by minimizing ${\mathrm{\Delta }}_1$ error. For real data case, one can use Morozov’s principle instead [31].

The deformation procedure (see Section 2.1) is included in the forward (fan-beam geometry) model $\mathrm{A}{\mathit{g}}^n$ (on every n-iteration, the object ${\mathit{g}}^n$ is deformed). The scheme (15(15) $$\begin{array}{c}\hfill {\mathbf{g}}^{n+1}={\mathbf{g}}^n+\mathit{\tau }{\mathrm{A}}_0^{-1}\mathrm{\Delta }{\mathbf{r}}^n,\end{array}$$(15) ) converges fast and also robust to noisy and undersampled data due to added TV term.

3. Numerical results

In this section, we present preliminary numerical results for three iterative GFST-based algorithms for fan-beam geometry (see Section 2.1). The compared algorithms include: FBGP (see Section 2.2), NDSL, NDSL-TV (see Section 2.4) and the classical Simultaneous Iterative Reconstruction Technique (SIRT) [1] for fan-beam geometry. Notably, in SIRT we do not use any of the proposed deformation principles. The reconstruction error ${\mathrm{\Delta }}_1(\%)$ (12(12) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_1=\sqrt{\frac{\Vert g-\widehat{g}{\Vert }_2}{\Vert \widehat{g}{\Vert }_2}}·100\%,\end{array}$$(12) ) is used to quantify the results.

3.1. Synthetic models and algorithms initialization

In order to test our algorithms accurately, we used two analytical phantoms with different properties. Model 1 is a smooth phantom and consists of six symmetric Gaussians (see Figure 5(a–b)). Model 2 is a piecewise-constant phantom and consists of six symmetric circles with uniform intensity (see Figure 5(c–d)). The compared algorithms can demonstrate a different behaviour to these models.

Figure 5. Synthetic phantoms used in numerical experiments: (a) model 1 (six symmetric Gaussians); (b) surface representation of model 1; (c) model 2 (six symmetric circles); (d) surface representation of model 2. The corresponding projection data is obtained analytically using the fan-beam geometry in Figure 1.

Figure 6. Various initializations of $g^{(0)}(x,y)$ for iterative algorithms. Reconstructions of model 2 from 11 projections ($K=11$) with noise level $\mathit{\kappa }=3\%$; (a) backprojection; (b) FBP; (c) The Baranov’s approximation [34].

It is also important to establish a good initialization (the first approximation to the seeking solution) for the proposed IIR schemes (see Algorithms 1(1) $$\begin{array}{c}\hfill \left\{\begin{array}{c}u=s/(1-\frac{p}{D})=s/Q;\hfill \\ v=p.\hfill \end{array}\right.\end{array}$$(1) –3(3) $$\begin{array}{cc}\hfill f_{\mathit{\beta }}(u)& ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\mathit{\delta }(q-U(x,y))g(x,y)\mathrm{d}x\mathrm{d}y\hfill \\ \hfill & ={\int }_{-\infty }^{\infty }{\int }_{-\infty }^{\infty }\left|J\right|\mathit{\delta }(q-U)g(u,v)\mathrm{d}u\mathrm{d}v.\hfill \end{array}$$(3) ). Potentially, if the first approximation to the solution is well chosen, it can have a positive impact on the speed of convergence and the final error. Here we compare four different initializations for our algorithms: zero image $g^{(0)}(x,y)=0$, backprojected image (see Figure 6 (a)), SL filtered backprojected (FBP) image (see Figure 6(b)), and initialization using the Baranov’s criteria of minimal projections [34] (see Figure 6(c)). The reconstructions in Figure 6 are performed using only 11 fan-beam projections (K) in $(0,2\mathit{\pi })$ angular range in radians with $3\%$ of randomly distributed noise ($\mathit{\kappa }$) added to the projections.

Our experiments have shown that the initial approximation for $g^{(0)}(x,y)$ using the Baranov’s criteria [34] (see Figure 6(c)) produces the best results (convergence is faster and the final error is lower) for all compared algorithms. Therefore, we will use this initialization in our experiments.

3.2. Numerical comparisons

We start our numerical simulations with a thorough search for all parameters to each compared algorithm. For the FBGP algorithm, the following parameters have been empirically established: the width of the strip $S_w=2h_p$, the multiplier ${\mathit{\delta }}_s=0.9h_p$ and the frequency of the width updates $m_0=1$ (see Section 2.3). For the NDSL and NDSL-TV algorithms, the most important parameters are the step size $\mathit{\tau }$ which controls the convergence of iterations and the regularization parameter $\mathit{\lambda }$ for the TV minimization in Algorithm 3. To find $\mathit{\tau }$, we performed several full reconstructions with different values of $\mathit{\tau }$ (see Figure 7). Notably, the parameter $\mathit{\lambda }$ has been found prior to this experiment and kept fixed to establish $\mathit{\tau }$ for the NDSL-TV method. We established that the optimal parameters $\mathit{\tau }=0.1$ and $\mathit{\tau }=0.07$ give the best results for the NDSL and NDSL-TV methods, respectively (see Figure 7). For the SIRT method, we use 150 iterations which is sufficient for the method to converge.

Figure 7. Optimization procedure to find the optimal step parameter $\mathit{\tau }$ for the NDSL (left) and NDSL-TV (right) methods. Model 1, $K=11$, $\mathit{\kappa }=3\%$.

Once all parameters have been identified through empirical optimization procedure for each algorithm, we proceed with numerical comparisons of methods with respect to various levels of noise ($\mathit{\kappa }=3\%,7\%$) and the number of projections ($K=11,22,33$) in $(0,2\mathit{\pi })$ angular range in radians. We apply reconstruction methods to both models 1 and 2. In Figure 8, we demonstrate the best achieved reconstruction errors ${\mathrm{\Delta }}_1(\%)$ for $\mathit{\kappa }=3\%$ and in Figure 9 for $\mathit{\kappa }=7\%$.

One can see that the NDSL-TV algorithm outperforms the unregularized NDSL method, FBGP and the SIRT methods. The NDSL-TV reconstruction results in the smallest errors for the limited data conditions $K\le 22$. In the case when more projections available ($K>22$), the efficiency of the NDSL-TV and the FBGP algorithms reduces (no significant enhancement achieved), while the NDSL algorithm has almost linear correlation of errors with respect to the number of projections.

Figure 8. Reconstructions of model 1 (left) and model 2 (right) using different number of projections ($K=11,22,33$) and the same level of noise $\mathit{\kappa }=3\%$. The NDSL-TV method outperforms (in terms of ${\mathrm{\Delta }}_1(\%)$) all other algorithms significantly.

Figure 9. Reconstructions of the model 1 (left) and model 2 (right) using the different number of projections ($K=11,22,33$) and the same level of noise $\mathit{\kappa }=7\%$. The NDSL-TV method outperforms (in terms of ${\mathrm{\Delta }}_1(\%)$) all other algorithms.

Figure 10. Reconstructions of model 1 (left) and model 2 (right) with respect to the iteration number and ${\mathrm{\Delta }}_1(\%)$ error. For both models the NDSL-TV method gives stable convergence to the point of stagnation.

To demonstrate the rate of convergence for the compared algorithms, we show the curves of residual errors with respect to the iteration number for both models (see Figure 10). One can see that the NDSL method is not stable and diverges quickly after ten iterations. The NDSL-TV method reaches the lowest error and ensures stable convergence to the point of stagnation (we iterated up to 150 iterations). The FBGP method is the quickest in convergence during first iterations; however, it slowly diverges afterwards and the minimal error is still higher than for the NDSL-TV method. The highest error is obtained with SIRT for both models. Notably we used the SIRT algorithm in its classical form with zeroed image as an initialization. Additionally, no a priori information is integrated into the algorithm compare to other methods. Consequently, it is expected that errors can be higher for SIRT.

In Figures 11 and 12, we demonstrate the reconstructed images for models 1 and 2, respectively. The reconstructions are linked to the curves presented in Figure 10, the images shown for the lowest errors achieved.

Figure 11. Reconstruction of the smooth model 1, $K=11$, $\mathit{\kappa }=3\%$ using the methods: (a) SIRT, ${\mathrm{\Delta }}_1=36.5\%$; (b) FBGP, ${\mathrm{\Delta }}_1=16.5\%$; (c) NDSL, ${\mathrm{\Delta }}_1=25.3\%$; (d) NDSL-TV, ${\mathrm{\Delta }}_1=11.2\%$. Note the flatten centres (peaks) of the reconstructed Gaussians with the NDSL-TV method.

Figure 12. Reconstruction of the discontinuous model 2, $K=11$, $\mathit{\kappa }=3\%$ using the methods: (a) SIRT, ${\mathrm{\Delta }}_1=43.2\%$; (b) FBGP, ${\mathrm{\Delta }}_1=29.6\%$; (c) NDSL, ${\mathrm{\Delta }}_1=35.7\%$; (d) NDSL-TV, ${\mathrm{\Delta }}_1=18.6\%$. Note the smooth shapes of the reconstructed circles with the FBGP method.

Since the compared methods are intrinsically different (especially the FBGP and the NDSL methods) they own different reconstruction characteristics and therefore the reconstructed images are dissimilar (see Figures 11 and 12). The smooth model 1 is reconstructed quite well with the FBGP and the NDSL-TV methods (see Figure 11). Notably, the NDSL-TV algorithm flattens the centres of Gaussians. This is mainly due to the piecewise-constant edge preserving properties of the TV penalty. Although, the total reconstruction error is higher with the FBGP method than with the NDSL-TV method the centres of Gaussians are recovered better with the FBGP method. SIRT reconstruction is full of aliasing artefacts; however, the shapes of Gaussians are preserved. The discontinuous model 2, however is reconstructed very poorly with the SIRT algorithm, image contains many artefacts and the uniform shapes of circles are distorted (oversmoothed). It is known that SIRT can implicitly smooth an image in iterations. The model is reconstructed much better with the NDSL-TV method due to piecewise-constant nature of the phantom itself (see Figure 12). The FBGP method fails to reconstruct the uniform structures of circles and smooths the boundaries (a bit similar to SIRT reconstruction, but not so extreme). The unregularized NDSL method, as expected, results in a very noisy reconstruction for both models.

4. Discussion

Although the performance of the FBGP algorithm has not been superior among the algorithms presented here, it can be further improved. The strip interpolation is one of the most important elements of the FBGP algorithm that affects its performance. Overall, the FBGP method remains a powerful and computationally efficient technique for reconstruction of the undersampled data. Arguably, the repetitive use of interpolation (even of higher orders) in the deformation procedure can affect the quality of the reconstructed images. In order to decide which interpolation is more suitable for the task and check if the error accumulates in iterations, one needs to initiate a thorough investigation of the issue. We will consider this in our future work. On this stage of research, the proposed framework can be used for fast prototyping and testing various IIR techniques (including the Fourier-based ones) for divergent-beam geometries.

5. Conclusion

In this work, we presented a novel class of iterative solutions for fan-beam geometry based on the generalized Fourier Slice Theorem. The proposed approach helps to generalize parallel-beam specific state-of-the-art reconstruction techniques to divergent-beam geometries. We have shown that the Fourier-based iterative methods, such as the Gerchberg-Papoulis algorithm can be adapted to reconstruct divergent-beam data using the proposed framework. This approach avoids rebinning step which can be undesirable for few-view tomography. The simple implementation of the generalized projection slice theorem makes it ideal for prototyping various iterative algorithms. The acceleration using the Fast Fourier Transform routines makes this approach more feasible than the classical iterative methods for reconstruction of big datasets.

Our preliminary numerical experiments have shown the potential in applicability of the generalized slice theorem to severely angularly undersampled reconstruction problems in X-ray tomography. Three new iterative methods for fan-beam scanning geometry have been introduced. Methods have been compared with each other and also with the classical algebraic fan-beam algorithm. The TV-regularized Neumann decomposition method has demonstrated the best performance among compared methods. Our future work will include improving these methods, applying to reconstruction of real data and generalization to cone-beam case.

Additional information

Funding

This work has been supported partly by the Engineering and Physical Sciences Research Council [grant number EP/J010553/1], [grant number EP/J010456/1], [grant number EP/I02249X/1].

Notes

No potential conflict of interest was reported by the authors.