Fast gridding projectors for analytical and iterative tomographic reconstruction of differential phase contrast data

Filippo Arcadu; Marco Stampanoni; Federica Marone

doi:10.1364/OE.24.014748

1. Introduction

The gridding method is a technique designed to retrieve a 2D or 3D signal from samples of its Fourier transform located on a non-Cartesian lattice [1]. This method has been applied for decades to image reconstruction problems in radio-astronomy [2, 3], magnetic resonance imaging [1, 4, 5] and absorption computer tomography (CT) [6–16].

A tomographic forward or backprojector based on the gridding method has O(N²log₂N) complexity (N is the number of rows/columns of an image supposed to be square), leading to a significantly higher computational efficiency than a space-based implementation of the same operator (O(N³) complexity) [17]. Several studies have also shown that CT analytical [9, 10, 15] and iterative [12, 13] algorithms, adopting gridding projectors, can provide high quality reconstructions, while substantially reducing the calculation times (up to 3 orders of magnitude).

This work introduces a gridding forward and backprojector designed to perform analytical and iterative reconstruction of tomographic datasets, where the forward model is represented by the derivative of the Radon transform. These operators are, in particular, applied to the reconstruction of differential phase contrast (DPC) tomographic data acquired with a grating interferometer.

A grating interferometer is an imaging setup yielding three complementary types of information [18–20]: attenuation, differential phase and dark field signals. DPC focuses exclusively on the phase information related to the real part of the object refractive index and provides a high sensitivity to electron-density variations, down to 0.18 e/nm³ [21]. For this reason, DPC based on a grating interferometer is particularly suitable for the visualization of soft-tissue specimens.

DPC data can be analytically reconstructed by means of filtered backprojection without need for integration of the projection information, if an imaginary filter function corresponding to a Hilbert transform in the image space (HFBP) is used instead of the usual Ram-Lak filter [19]. This approach is however not suited for the reconstruction of DPC underconstrained data because of insufficient accuracy. Recent interest in the biomedical field for grating interferometry, when investigating dose sensitive specimens, has lead to several studies aimed at extending established CT iterative algorithms to the DPC case. In this regard, ad-hoc forward projectors working with a blob [22] and a cubic B-spline [23] basis and iterative schemes relying on different algorithms (the regularized maximum likelihood [22], the separable paraboloidal surrogate [24], the alternate direction method of multipliers [25] and a combination of statistical model and thresholding [26]) have been proposed.

So far, all analytical and iterative reconstruction methods developed for DPC utilize space-based projectors, characterized by a complexity of O(N³) and, therefore, strongly limiting their computational performance. The projectors proposed here greatly reduce the computational times required for analytical and, in particular, iterative reconstruction, while preserving the accuracy of the results.

1.1. Contributions

The contributions of this manuscript are summarized as follows:

Design of a gridding forward projector acting as derivative of the Radon transform.
Optimization of the parameters by means of the minimal oversampling strategy [5].
Comparison of accuracy and efficiency with a standard space-based implementation of the same projector.
Validation of the proposed operators through analytical and iterative reconstruction of simulated and experimental DPC data.

2. Notation and preliminaries

This work focuses on parallel beam geometry, even though gridding can be extended to fan-and cone-beam geometries [14–16].

The study object, i.e., an axial 2D slice of the whole 3D volume, is represented by a finite integrable real function f = f(x) = f(x, y): Ω ⊆ ℝ² → ℝ with bounded support S = supp(f) ≡ {x ∈ ℝ² |f(x) ≠ 0}.

Its CT forward projection or X-ray transform or Radon transform (in 2D, these transforms coincide [27]) is given by:

P_{θ} (t) : = ℛ {f (x)} (θ, t) : = \int_{Ω} d x f (x) δ (x \cdot e_{θ} - t),

where · represents the inner product, θ ∈ [0, π) and e_θ = (cosθ, sinθ). The variables used in (1) are displayed in Fig. 1(a).

Fig. 1 (a) Coordinate system and variables used. (b) DPC sinogram example.

Download Full Size | PDF

The derivative of the Radon transform is considered with respect to the variable t and is given by:

\frac{d}{d t} P_{θ} (t) : = ℛ^{(1)} {f (x)} (θ, t) : = \frac{d}{d t} \int_{Ω} d x f (x) δ (x \cdot e_{θ} - t) .

The collection of all projections (P_{θ} or dP_{θ}/dt) in parallel beam geometry is also called sinogram. Figure 1(b) shows an example of DPC sinogram; the pixels of each projection are placed along the row or channel direction.

The 2D Fourier transform of a generic function g(a, b) is indicated with ĝ or ℱ_a,b{g}; ℱ_a{g} or ℱ_b{g} represents the 1D Fourier transform along the variable a or b, respectively. (x, y), (t, θ), (u, v) and (θ, ω) are the coordinates of the real-space Cartesian, real-space polar, Fourier Cartesian and Fourier polar reference frame, respectively. (···) refers to entries of a continuous function, [···] entries of a discrete function/array and {···} either entries of an operator or a collection of discrete values. Vectors and matrices are indicated with lower- and upper-case bold letters, respectively; both ^† and * correspond to the Hermitian adjoint operator.

3. Proposed method

3.1. Connection between gridding and CT

The gridding method (GM) retrieves a signal from samples of its Fourier transform located on a non Cartesian lattice: the samples are, first, convolved with a smooth and rapidly decaying window function, the inverse fast Fourier transform (IFFT) is computed and, finally, the contribution of the window function is removed from the signal in real space [1].

The connection between GM and CT reconstruction is provided by the Fourier slice theorem (FST) [27]:

{\hat{P}}_{θ} (ω) = \hat{f} (ω cos θ, ω sin θ) .

The discretized version of (3) states that the fast Fourier transform (FFT) of a projection at angle θ yields equidistant samples of f̂ along the line v = u tan θ. The samples collected from multiple projections form a polar lattice in Fourier space. To retrieve f through IFFT, these polar samples need to be interpolated at Cartesian locations, an operation that can be performed by the GM. CT reconstruction through a filtered gridding backprojector consists in the following steps [6, 8, 10]:

FFT1 of the sinogram along the channel direction ← P̂_{θ};
ramp filtering of the sinogram $\leftarrow {\hat{P}}_{{θ}}^{(f)}$ ;
convolution with window function $\leftarrow {\hat{f}}^{(d)} = \hat{h} * {\hat{P}}_{{θ}}^{(f)}$ ;
IFFT2 of f̂^(d) ← f^(d);
removal of the window function ← f = f^(d)/h.

Accordingly, the gridding-based implementation of ℛ or forward gridding projector (FRP) is given by the same steps listed above (ramp filtering excluded) in the reverse order [9, 13]:

pre-correction ← f^(d) = f/h;
FFT2 of f^(d) ← f̂^(d);
convolution with window function → P̂_{θ} = ĥ * f̂^(d);
IFFT1 of each polar slice of P̂_{θ} ← P_{θ}.

3.2. Differentiated forward gridding projector

The formula of the differentiated forward gridding projector (DFRP) is obtained from the FST (3), the definitions (1), (2) and from the property that derivatives in real space are mapped into multiplications in Fourier space.

\begin{array}{l} ℛ^{(1)} {f} (θ, t) & = \frac{d}{d t} P_{θ} (t) = \frac{d}{d t} \int_{Ω} d ω {\hat{P}}_{θ} (ω) e^{- 2 π i ω t} = \int_{Ω} d ω e^{- 2 π i ω t} (- 2 π i ω) \cdot (\hat{h} * {\hat{f}}^{(d)}) (θ, ω) \\ = ℱ_{ω}^{- 1} {((- 2 π i ω))} \cdot (\hat{h} * ℱ_{x, y} {f / h (x, y)}) (θ, ω), \end{array}

where i is the imaginary unit (i² = −1), h and ĥ are the gridding interpolation kernel and its Fourier transform, f^(d) = f/h is the pre-corrected and oversampled version of f. Choosing the kernel separable, i.e., h(x, y) = ψ(x)ψ(y) ⇒ ĥ(u, v) =ψ̂(u)ψ̂(v), simplifies the implementation of the method. Discretizing (4) results in:

\begin{array}{l} \frac{d}{d t} P_{θ_{m}} & = \frac{d}{d t} ℱ_{{ω}}^{- 1} {{\hat{P}}_{θ_{m}} [ω_{j}]}_{n} = \frac{d}{d t} ℱ_{{ω}}^{- 1} {(\hat{h} * {\hat{f}}^{(d)}) [ω_{j}, θ_{m}]}_{n} \\ = \frac{d}{d t} ℱ_{{ω}}^{- 1} {\sum_{\begin{array}{l} u_{r} \in U_{m j} \\ v_{s} \in V_{m j} \end{array}} \underset{C_{m, j, r, s}}{\underset{︸}{\hat{ψ} [g_{1} (ω_{j, i} - u_{r})] \hat{ψ} [g_{1} (ω_{j, v} - v_{s})]}} {\hat{f}}^{(d)} [u_{r}, v_{s}]} \\ ℱ_{{ω}}^{- 1} {(- 2 π i ω_{j}) \sum_{\begin{array}{l} u_{r} \in U_{m j} \\ v_{s} \in V_{m j} \end{array}} C_{m, j, r, s} ℱ_{{x, y}} {\frac{f^{(o)} [x_{p}, y_{q}]}{ψ [g_{2} x_{p}] ψ [g_{2} y_{q}]}}_{u_{r}, v_{s}}}_{n} \end{array}

where i is the imaginary unit; θ_m = mπ/M ∈ [0, π) for m = 0, 1,...,M −1 and M is number of views; the image is assumed to be square with N pixels (N is even) and to have unit resolution, therefore, x_p = −N/2 + p and y_q = −N/2 + q for p, q = 0, 1,..., N − 1 are the image pixel coordinates; t_n = −N/2+n for n = 0, 1, ..., N −1 (N is also the number of detector cells); the Fourier space is sampled at G evenly spaced points in the range [0.5, 0.5), therefore, u_r = −1/2 + r/G and v_s = −1/2+s/G for r, s = 0, 1, ..., G−1, G = αN for α > 1; α is called oversampling ratio; ω_j,u = ω_j cosθ_m, ω_j,v = ω_j sinθ_m and ω_j = −1/2 + j/G for j = 0, 1, ..., G − 1; U_mj (V_mj) are the sets of pixel coordinates u_r (v_s) inside the interpolation support; g₁ and g₂ are factors mapping distances from the polar and Cartesian grid into look-up table (LUT) indices and depend on the sampling density of ψ and ψ̂.

ℱ_{ω}^{- 1}

and ℱ_{x,y} in (5) are implemented as IFFT1 along ω and FFT2 along x, y, respectively. f^(o) represents a zero-padded or oversampled version of f by a factor α > 1.0. Gridding an oversampled grid prevents the appearance of aliasing artifacts due to tails of the reconstruction wrapping back into the field-of-view (FOV) [7, 10]. A general rule of thumb is to set α = 2.0.

Since the DFRP is a linear operator, the differentiated gridding backprojector (DBRP) is the Hermitian adjoint operator of (5). From a computational point of view, the DBRP performs the same exact operations as the DFRP, but in reverse order.

The accuracy and efficiency of gridding operators depend on the choice of the finite kernel ψ: the smaller the support of ψ, the more accurate the pre-correction in real space; at the same time, the smaller the support of ψ̂, the faster the interpolation stage in Fourier space [7, 10, 13]. The accuracy of a convolution kernel is also strongly determined by the shape (in particular, the rolloff) of the central lobe and the amplitude of the aliasing sidelobes characterizing ψ. The pre-correction removes the rolloff induced by this central lobe, but, in doing so, the aliasing sidelobes are amplified [7, 13].

In previous work [4, 6, 7], different compact kernels have been studied for gridding [7]. The conclusion is that a Kaiser-Bessel (KB) [28] kernel has an optimal compact support in real and Fourier space and provides analytical formulas, that enable shape optimization for ψ and ψ̂ on the basis of a certain criterion, if available. The formulas for the KB kernel are:

ψ_{kb} (x) = \frac{sin \sqrt{{(\frac{π W x}{G})}^{2} - β^{2}}}{\sqrt{{(\frac{π W x}{G})}^{2} - β^{2}}}, {\hat{ψ}}_{kb} (k_{x}) = \frac{G}{W} I_{0} (β \sqrt{1 - {(\frac{2 G k_{x}}{W})}^{2}}),

valid for |k_x| ≤ W/2G and where I₀ is the zero-th order modified Bessel function, W is the size of the convolving kernel and β is the tapering parameter, that determines how fast the KB drops to zero. The accuracy and computational performance of (5) with a KB kernel depends on α, W, β and S, the sampling density of ψ̂. Working with a pre-sampled LUT enables faster calculations and this approach is, therefore, preferable than invoking Eq. (6), every time a point needs to be interpolated in Fourier space.

According to the minimum oversampling criterion, the KB’s shape can be optimized by minimizing the maximum pixel aliasing error [5]:

ε_{i j} = \sqrt{\frac{\sum_{p \neq 0} (h {[i + G p, j + G p]}^{2})}{{(h [i, j])}^{2}}} = \sqrt{\frac{\sum_{p \neq 0} {(ψ [i + G p] ψ [j + G p])}^{2}}{{(ψ [i] ψ [j])}^{2}}} .

The optimal kernel ψ^opt is found as solution of the minimization problem:

ψ^{opt} = \underset{ψ}{arg min} max_{i, j} (ε_{i, j}),

which translates in the following results [5]:

β = π \sqrt{{(\frac{W}{α})}^{2} \cdot {(α - \frac{1}{2})}^{2} - 0.8}, s = {\begin{array}{l} 0.91 / (α \cdot γ) & for NN \\ \sqrt{0.37 / (γ α^{2})} & for LIN \end{array}

where NN and LIN refer, respectively, to the nearest neighbor and linear interpolation scheme to sample the convolution LUT and γ is the maximum allowed sampling error for ψ̂. The minimal oversampling criterion (9) reduces the initial parameter space {W, S, α, β} to {W, γ, α}.

3.3. Algorithm complexity

The cost of the DFRP lies in the convolution and the call of FFT/IFFT. Given an input image of N × N pixels, an oversampling ratio α, a kernel width W and a number of views M, the convolution amounts to approximately WMN operations, whereas the overall call of FFT/IFFT corresponds to α²N² log₂(αN) + αMN log₂(αN) (1-time FFT-2D and M-times IFFT-1D) floating point operations. This second factor usually represents the leading cost term for real datasets.

4. Optimization of the gridding parameters

4.1. Optimization approach

For the pure Radon transform case, an optimal KB kernel, ψ^opt, can be constructed with any setting {W, S, α, β} fulfilling (9). The choice of one specific set among the optimal ones would then be dictated by computational efficiency criteria. However, the minimal oversampling criterion does not account for the influence of the term −2πiω_j on the accuracy of the actual implementation of the DFRP. To assess which set of theoretically optimal gridding parameters provides the most accurate results for the DFRP, the 3D parameter space {W, γ, α} needs to be explored in a brute-force manner, through experiments making use of a family of functions, whose ℛ⁽¹⁾ can be computed analytically. The following family of 2D radially symmetric functions is considered:

f^{(n)} (x) = {\begin{array}{l} {(a^{2} - {‖ x ‖}^{2})}^{n} & if & ‖ x ‖ \leq a \\ 0 & if & ‖ x ‖ > a \end{array} \forall n \in ℕ .

The analytical differentiated Radon transform is given by:

ℛ^{(1)} {f^{(n)}} (θ, t) = {\begin{array}{l} - 2 t \frac{(2 n)!!}{(2 n - 1)!!} {(a^{2} - t^{2})}^{n - 1 / 2} & if & n = 2 k \\ - π t \frac{(2 n)!!}{(2 n - 1)!!} {(a^{2} - t^{2})}^{n - 1 / 2} & if & n = 2 k + 1 \end{array} k \in ℕ, \forall θ \in [0 . π) .

To derive Eq. (11), one has to apply definition (2) to f⁽ⁿ⁾(x) and perform a cosinusoidal substitution to obtain an integral of the form

\int_{0}^{π} d ω {(sin ω)}^{m}

, whose solution is available in [29]. Optimal gridding parameters are evaluated for analytical and iterative reconstruction methods. For the former case, the Hilbert-filtered DBRP is applied to DPC sinograms, created with (11), and the obtained reconstructions are compared with the original phantom (10). For the latter case, an approach disregarding the specific choice of iterative algorithm (and its specific parameters) is followed: the forward projection computed with DFRP is compared to the analytical sinogram and is, then, reconstructed with the Hilbert-filtered DBRP. In this way, both the accuracy of the standalone forward projector and of the coupling forward-adjoint operator can be investigated.

The experiments in Sec. 4.2 make use of f⁽ⁿ⁾(x) with n = {2, 3}, sampled on a grid with 512 × 512 pixels and its analytical ℛ⁽¹⁾ is computed for 805 views, homogeneously distributed in [0, π). The number of views is enough to prevent undersampling artifacts in the reconstructions, since 805 ⋍ 512 * π/2 (according to the FBP Nyquist criterion [27]). Results are not biased by the resolution of the phantoms, because experiments repeated with different grid sizes have yielded the same outcome.

The accuracy of the results is measured by the peak signal-to-noise ratio (PSNR) [30], a score quantifying the overall difference between the computed and reference image. The PSNR is defined as:

PSNR = 10 \cdot {log}_{10} [\frac{max_{i, j} {(O [i, j])}^{2}}{\frac{1}{M N} \cdot \sum_{i = 0}^{M - 1} \sum_{i = 0}^{N - 1} {(O [i, j] - I [i, j])}^{2}}]

where O and I are 2D arrays with M rows and N columns, being, respectively, the oracle and the image to analyze. Since in (12) the difference is located at the denominator, PSNR results are more sensitive to small variations from the oracle than the mean square error. When analyzing sinograms, the PSNR is computed over the entire image, whereas, for reconstructions, it is computed within the resolution circle.

The mean structural similarity index (MSSIM) [31] has also been employed for the image analysis. Although confirming the trends described by the PSNR, the MSSIM results poorly sensitive for this kind of test. For this reason, only PSNR scores are reported in Sec.4.2.

4.2. Optimal gridding parameters

To limit the extent of the brute force search of {W, γ, α}_opt, preliminary experiments, studying how the accuracy of the operators depends on just one parameter, while fixing the other two, have been performed. These experiments focus on the forward projection with DFRP and on the reconstruction of analytical sinograms with Hilbert-filtered DBRP. Results, collected in Figs. 2(a)–2(c), show the following key trends: the choice of α is more critical for the backprojector (Fig. 2(b)); the PSNR decreases only when γ > 10⁻⁴ (Fig. 2(a)); the DFRP is highly influenced by the choice of W (Fig. 2(c)).

Fig. 2 Preliminary experiments to determine the interesting range where to span each parameter in a more general brute-force search. The DFRP projects f⁽ⁿ⁾ with n = 2, 3: the result is compared to the analytical sinogram. The Hilbert-filtered DBRP reconstructs the analytical sinogram: the reconstructed slice is compared to the phantom. (a) Experiments spanning γ ∈ [10⁻¹⁰, 3.9 · 10⁻³] and setting W = 4.45 and α = 3.0; all PSNR curves reach the maximum for γ < 10⁻⁴ (corresponding to S > 50). (b) Experiments spanning α ∈ [1.05, 3] and setting W = 4.45 and γ = 10⁻¹⁰; all PSNR curves reach the maximum for α > 1.8. (c) Experiments spanning W ∈ [1.9, 12.1] and setting γ = 10⁻¹⁰ and α = 3.0.

Download Full Size | PDF

The general brute force search of the 3D parameter space has W ranging in [2.7, 17.7] (20 points), γ in [10⁻⁸, 10⁻³] (10 points) and α in [1.05, 3.0] (20 points). Results can be visualized by means of 2D PNSR maps, displaying the accuracy as a function of two parameters while the third one is fixed, as shown in Figs. 3(a)–3(c). Few relevant trends, mirroring the results of the preliminary experiments, can be recognized. The accuracy of DFRP and DBRP is only weakly influenced by γ, provided that this parameter is < 10⁻⁴ (Figs. 3(a) and 3(b)). DFRP is very sensitive to the choice of W (Fig. 3(b)): the PSNR varies up to 5 dB for a difference of 0.8 in kernel size. The reconstruction accuracy of DFRP sinograms (coupling DFRP-DBRP) is sensitive to the choice of α (Fig. 3(a)) in a wider interval, 1.05 ⩽ α ⩽ 2.30, than for the reconstruction of analytical sinograms (DBRP standalone, Fig. 3(c)), 1.05 ⩽ α ⩽ 1.70.

Fig. 3 Examples of 2D PSNR maps visualizing the results of the brute-force search of the 3D parameter space {W, γ, α}. These experiments made use of the phantom f⁽²⁾. (a) Map for the reconstruction with Hilbert-filtered DBRP of DFRP sinograms. (b) Map for the computation of DFRP sinograms. (c) Map for the reconstruction with Hilbert-filtered DBRP of analytical sinograms.

Download Full Size | PDF

The optimal parameters are found by locating the maximum PSNR values in the 3D volume, created by the brute force search. The PSNR reaches the maximum on a continuous spread area (the intersection of the darkest blue regions in Figs. 3(a)–3(c)), while no local extrema exist. Since the broad PSNR maximum corresponds to multiple parameter configurations, {α, W, γ}_opt is chosen as the peripheral point, requiring the smallest computational time. The optimal parameters for analytical reconstruction, as obtained from the top score for the reconstruction with Hilbert-filtered DBRP of analytical sinograms, are {W, γ, α}_opt₁ = {4.45, 1.7 · 10⁻⁶, 1.75}. The optimal parameters for iterative reconstruction, merging the results for the computation of the forward projection with DFRP and the recostructions of these sinograms with Hilbert-filtered DBRP, are {W, γ, α}_opt₂ = {6.6, 6.0 · 10⁻⁶, 2.38}.

The second set of optimal parameters is more conservative and less computationally efficient, because it involves a bigger kernel size and a higher oversampling ratio compared to the first set. The aliasing error introduced by the DFRP can be amplified by the DBRP at each iteration in iterative schemes, therefore a more conservative choice of parameters is mandatory. The results of the brute-force search were confirmed, when using f⁽ⁿ⁾ with n = 1, 4, 5.

Additional experiments with the Shepp-Logan (SL) [32] phantom have proved that this optimization approach is not biased by the choice of the family function. In the experiment of Tabs. 1(a) and 1(b), the settings {W, γ, α}_nopt₁ and {W, γ, α}_nopt₂, belonging to the PSNR maximum surface for analytical and iterative reconstruction respectively, are tested against the optimal settings. Table 1(a) shows the PSNR score related to the reconstruction with Hilbert-filtered DBRP of the SL analytical sinogram (800 views × 512 pixels) with {W, γ, α}_opt₁ and {W, γ, α}_nopt₁ = {7.6, 1.7 · 10⁻⁷, 2.1}. The PSNR score in Table 1(b) refers to the reconstruction with Hilbert-filtered DBRP of the sinogram (800 views × 512 pixels) computed by the DFRP with {W, γ, α}_opt₂ and {W, γ, α}_nopt₂ = {7.6, 6.0 · 10⁻⁷, 2.7}. The negligible differences in PSNR confirm, that {W, γ, α}_nopt₁ and {W, γ, α}_nopt₂ correspond indeed to the respective PSNR maximum surfaces. The same approach is followed for the experiment in Tabs. 2(a) and 2(b), where, this time, the settings {W, γ, α}_nopt₃ = {11.9, 6.0 · 10⁻⁷, 1.66} and {W, γ, α}_nopt₄ = {8.2, 6.0 · 10⁻⁷, 2.2} are not located on the maximum surfaces for analytical and iterative reconstruction, respectively. The differences of PSNR in Tabs. 2(a) and 2(b), between optimal and suboptimal settings, are certainly not negligible.

Table 1. (a) Reconstruction with Hilbert-filtered DBRP of an analytical SL sinogram (800 views × 512 pixels). (b) Reconstruction with Hilbert-filtered DBRP of a sinogram created by the DFRP (800 views × 512 pixels). The settings {W, γ, α}_nopt₁ and {W, γ, α}_nopt₂ belong to the PSNR maximum surface for analytical and iterative reconstruction, respectively. This is confirmed by the negligible differences in PSNR.

View Table | View all tables in this article

Table 2. (a) Reconstruction with Hilbert-filtered DBRP of an analytical SL sinogram (800 views × 512 pixels). (b) Reconstruction with Hilbert-filtered DBRP of a sinogram created by the DFRP (800 views × 512 pixels). The settings {W, γ, α}_nopt₃ and {W, γ, α}_nopt₄ do not correspond to the PSNR maximum surface for analytical and iterative reconstruction, respectively. The differences in PSNR are not negligible.

View Table | View all tables in this article

5. Comparison with a standard implementation of ℛ⁽¹⁾

5.1. Benchmark procedure

To benchmark the accuracy and efficiency of the proposed operators, a standard space-based implementation of ℛ⁽¹⁾ has been considered: the forward projector introduced by [17] with additional differentiation along the channel direction. Considered a discrete image f[m, n], the original implementation of the Radon Transform of [17], based on linear interpolation, is given by:

P θ (t) = {\begin{array}{l} \frac{1}{sin θ} \sum_{m} [(1 - η) f [m, ⌊ a m + b ⌋] + η f [m, ⌊ a m + b ⌋ + 1]] & if & | sin θ | > \frac{\sqrt{2}}{2} \\ \frac{1}{cos θ} \sum_{n} [(1 - η) f [⌊ \tilde{a} n + \tilde{b} ⌋, n] + η f [⌊ \tilde{a} n + \tilde{b} ⌋ + 1, n]] & if & | sin θ | \leq \frac{\sqrt{2}}{2} \end{array}

where

\begin{array}{l} \begin{array}{l} a = - \frac{cos θ}{sin θ} & b = \frac{t - x_{min} cos θ - y_{min} sin θ}{sin θ} \end{array} \\ \begin{array}{l} \tilde{a} = \frac{1}{a} & \tilde{b} = \frac{t - x_{min} cos θ - y_{min} sin θ}{cos θ} \end{array} \\ η = {\begin{array}{l} a m + b - ⌊ a m + b ⌋ & if & | sin θ | > \frac{\sqrt{2}}{2} \\ \tilde{a} n + \tilde{b} - ⌊ \tilde{a} n + \tilde{b} ⌋ & if & | sin θ | \leq \frac{\sqrt{2}}{2} . \end{array} \end{array}

⌊...⌋ represents the floor operation, x_min and y_min correspond to the smallest abscissa and ordinate of the selected ray within the image grid. ℛ⁽¹⁾ is computed from finite differences of Eq. (13) along the pixel direction, i.e. the t-axis. This operator is named differentiated Radon transform and abbreviated as DRT; its corresponding backprojector is abbreviated as DBP.

Experiments focus on the the impact of undersampling and noise in analytical and iterative reconstructions performed with the proposed operators on one side, DRT and DBP on the other side.

The PSNR and MSSIM are used for the assessment of the reconstruction quality of simulated data. For real data, since a reference image is not available, the signal-to-noise ratio (SNR) and contrast-to-noise ratio (CNR) are averaged for multiple regions-of-interest (ROIs) at different distances from the image center. The SNR is computed within a ROI, supposed to be homogeneous; for the CNR, two different neighbouring ROIs are used, each supposed to be homogeneous. SNR and CNR are less reliable metrics than those relying on a reference image, but they can, nevertheless, give an approximate indication of image quality.

5.2. Validation of the proposed operators for analytical reconstruction

Analytical reconstructions are performed by the Hilbert-filtered DBRP (HF-DBRP) with optimal parameters {W, γ, α}_opt₁ = {4.45, 1.7 · 10⁻⁶, 1.75} and Hilbert-filtered DBP (HF-DBP).

The first experiment deals with a simulated undersampled and noisy dataset of f⁽ⁿ⁾ (n = 1, 2, 3) with 512 × 512 pixels, computed by means of (11). A set of noiseless analytical sinograms is created with a number of views ranging from 20 to 600. Another set is, instead, created with an increasing amount of Poisson noise, with σ ranging from 1% to 50% of the average value, μ, of the sinogram, while the number of views is fixed to 805.

Since reconstructions with HF-DBRP and HF-DBP, for the same condition of either under-sampling or noise, are indistinguishable at visual inspection, the analysis of the results relies on the information provided by PSNR and MSSIM (Figs. 4(a)–4(d)). Figures 4(a) and 4(b) show that the reconstructions with HF-DBRP have generally superior quality for a number of views ≥12% of what is required by the FBP sampling criterion (in this case, 100 views). For very low number of views, the HF-DBRP scores slightly worse than the space-based counterpart. Figures 4(c) and 4(d) show that the gridding algorithm appears quite sensitive to the level of noise affecting the data, since the drop of PSNR and SSIM for increasing noise is steeper compared to the other operator. The PSNR scores are everywhere higher for the gridding reconstruction, whereas, according to the MSSIM, for σ > 20%, the reconstructions with HF-DBP seem to perform better.

Fig. 4 PSNR and SSIM results for the comparison between the Hilbert-filtered DBRP and DBP for the analytical reconstruction of DPC undersampled and noisy sinograms. (a,b) Scores for the “undersampling-test ” (for f⁽¹⁾). (c,d) Scores for the “noise-test ” (for f⁽³⁾); σ on the x–axis is the variance of the Poisson noise corresponding to a percentage of the average value of the noiseless original sinogram.

Download Full Size | PDF

These results indicate that analytical reconstruction with the gridding backprojector is less sensitive to undersampling than to noise. The HF-DBRP may slightly underperform compared to a standard method, when a substantial amount of noise affects the data.

The second experiment tackles the reconstruction of a DPC sinogram of a rat brain, acquired at the TOMCAT beamline of the Swiss Light Source at the Paul Scherrer Institute (Switzerland). The array has 721 views × 1493 pixels. The reconstructions, shown in Figs. 5(a) and 5(b), are indistinguishable at visual inspection; nevertheless, the reconstruction with HF-DBRP shows improved SNR and CNR values (Table 3(a)). When the original sinogram is downsampled to 400 and 300 views, each pair of reconstructions still looks identical and the difference of SNR and CNR becomes rather small (Table 3(b) and 3(c)).

Fig. 5 Analytical reconstructions of the DPC sinogram of a rat brain with a zoomed area. (a) Reconstruction with the Hilbert-filtered DBRP. (b) Reconstruction with the Hilbert-filtered DBP. The light blue and yellow squares in (b) indicate the ROIs used to computed the average SNR and CNR values reported in Table 3.

Download Full Size | PDF

Table 3. SNR and CNR scores for the analytical reconstructions of the rat brain with 721 (a), 400 (b) and 300 (c) views. These values represent averages computed over the ROIs displayed in Fig. 5(b).

View Table | View all tables in this article

The computational efficiency represents the very attractive feature of the proposed operators. Table 4 reports the time elapsed (in seconds) for the reconstruction of DPC sinograms of various sizes. The used hardware is an Intel(R) Core(TM) i7-3520M CPU 2.90GHz. The HF-DBRP can speed-up calculations by two orders of magnitude for small datasets, up to three orders for bigger datasets with respect to a standard space-based operator.

Table 4. Times elapsed to perform the analytical reconstructions of DPC sinograms of different sizes. The used hardware is an Intel(R) Core(TM) i7-3520M CPU 2.90GHz.

View Table | View all tables in this article

5.3. Validation of the proposed operators for iterative reconstruction

Iterative reconstructions are performed by means of the alternate direction method of multipliers (ADMM) [33], solving a standard LASSO (Least Absolute Selection and Shrinkage Operator) problem. This method has been applied to DPC data in [25] with a more sophisticate formulation and a forward projector based on a cubic B-spline basis [23]. Details regarding the ADMM, its usage for tomographic reconstruction and the choice of the parameters are given in the Appendix.

Two implementations of the ADMM are tested in this section and the difference lies in the forward and backprojector used: the DFRP and its Hermitian adjoint operator (ADMM-DFRP) versus the DRT and its Hermitian adjoint operator (ADMM-DRT).

The first experiment focuses on the convergence of the algorithm: an analytical sinogram, created from f⁽²⁾ and having 100 views × 256 pixels with additional Poisson noise of σ = 8%μ (this σ has been chosen to be slightly bigger than that estimated on the real data presented in the following), is reconstructed by means of ADMM and the relative norm of the difference between consecutive iterations, i.e. ${‖ x^{(k + 1)} - x^{(k)} ‖}_{2}^{2} / {‖ x^{(k)} ‖}_{2}^{2}$ , is measured. The ADMM runs for 100 iterations. Figure 6 shows that ADMM-DFRP and ADMM-DRT have almost exactly the same convergence speed and that the solution remains practically unaltered after 20 iterations. This fact proves that the gridding operators do not alterate the algorithm convergence with respect to standard projectors. For the next reconstructions, the ADMM is stopped when ${‖ x^{(k + 1)} - x^{(k)} ‖}_{2}^{2} / {‖ x^{(k)} ‖}_{2}^{2} < ε = 5 \cdot 10^{- 6}$ .

Fig. 6 Convergence test showing the decrease of log₁₀ $({‖ x^{(k + 1)} - x^{(k)} ‖}_{2}^{2} / {‖ x^{(k)} ‖}_{2}^{2})$ as a function of the number of iterations.

Download Full Size | PDF

In a second step, the iterative reconstructions of the noisy and undersampled simulated sinogram, employed for the convergence test, are analyzed. Images in Figs. 7(b) and 7(c) show no visible difference; the PSNR and SSIM analysis (Table 5(a)) reveals a slightly better accuracy for the ADMM-DFRP.

Fig. 7 Iterative reconstructions of a noisy undersampled DPC sinogram of f⁽²⁾. (a) Original phantom. (b) Reconstruction with ADMM-DFRP. (c) Reconstruction with ADMM-DRT.

Download Full Size | PDF

Table 5. (a) PSNR and MSSIM scores for the iterative reconstructions shown in Figs. 7(b) and 7(c). (b) Average SNR and CNR computed for multiple ROIs at different distances from the image center for the iterative reconstructions shown in Figs. 8(a) and 8(b). The SNR and CNR scores were averaged over the ROIs displayed in Fig. 5(b).

View Table | View all tables in this article

Finally the sinogram of the rat brain, also used in the previous section, is reconstructed with the ADMM approach. In this case, the number of projections has been reduced by a factor of 2.5 leading to an undersampled dataset with an array of 300 views x 1493 pixels. Also for real data, the reconstructed slices (Figs. 8(a) and 8(b)) look almost identical at visual inspection; the analysis, reported in Tabs. 5(a) and 5(b), highlights a slight better performance for the ADMM-DFRP.

Fig. 8 Iterative reconstructions of a DPC undersampled sinogram of a rat brain with a zoomed area. (a) Reconstruction with ADMM-DFRP. (b) Reconstruction with ADMM-DRT.

Download Full Size | PDF

From a computational point of view, the proposed operators are ideal for iterative algorithms, since these methods generally call the forward projector and its Hermitian adjoint few times per iteration, often representing the computational bottleneck of the entire reconstruction process. For this last dataset, on an Intel(R) Core(TM) i7-3520M CPU 2.90GHz, a single conjugate-gradient sub-iteration within the ADMM (see Appendix) takes around 40s with the proposed operators (ADMM-DFRP), and around 830s for ADMM-DRT. Since the stopping criterion terminates the procedure after approximately 20 iteration, each involving 15 sub-iterations, the time elapsed for the reconstructions, if run on a single core, is around 3.5h for the ADMMDFRP and 69h for the ADMM-DRT.

6. Summary

Novel projectors have been introduced to address analytical and iterative reconstruction of tomographic differential phase contrast datasets, when the forward model is represented by the derivative of the Radon transform, ℛ⁽¹⁾. These operators are based on the gridding method and feature a complexity of O(N² log₂ N). The gridding parameters, offering the best tradeoff between accuracy and computational performance, have been experimentally optimized on the basis of the minimal oversampling strategy, that requires the utilization of a Kaiser-Bessel kernel for the interpolation in Fourier space.

The performance of the proposed differentiated gridding forward projector (DFRP) and back-projector (DBRP) have been assessed with respect to a standard space-based implementation of ℛ⁽¹⁾ and ℛ^{(1)^*} . The comparison has shown that the differentiated gridding operators can yield analytical and iterative reconstructions of the same quality as space-based operators. The great advantage of the proposed operators lies in the computational efficiency, since calculation times can be reduced by 3 orders of magnitude.

Appendix LASSO solved by the alternate direction method of multipliers

The LASSO problem corresponds to a L₁ penalized linear regression of the form:

\hat{x} = \underset{x}{arg min} [{‖ Ax - b ‖}_{2}^{2} + λ {‖ x ‖}_{1}] .

Considering the reconstruction problem presented in Sec.5.3, A, i.e., A = ℛ⁽¹⁾, is the matrix representation of the forward projector (therefore, A^† is the Hermitian adjoint operator or back-projector), i.e., b is the DPC sinogram, x is the unknown slice of the object and λ ∈ ℝ is the regularization weight.

The LASSO problem can be numerically tackled by the alternate direction method of multipliers (ADMM) [33], that maps (15) into the minimization of the following augmented Lagrangian:

ℒ_{μ} (x, u, m) = \frac{1}{2} {‖ Ax - b ‖}_{2}^{2} + λ \sum_{k} {‖ u_{k} ‖}_{1} + m^{T} (x - u) + \frac{μ}{2} {‖ x - u ‖}_{2}^{2}

where m are the Lagrangian multipliers and u is an auxiliary variable or the dual image.

The ADMM iteratively minimizes ℒ_μ by sequentially solving smaller problems; each iteration involves two sub-optimizations with respect to x and to u, followed by the update of m:

\begin{array}{l} 1 . x^{(k + 1)} \leftarrow \underset{x}{arg min} ℒ_{μ} (x, u^{(k)}, m^{(k)}) \\ 2 . u^{(k + 1)} \leftarrow \underset{u}{arg min} ℒ_{μ} (x^{(k + 1)}, u, m^{(k)}) \\ 3 . m^{(k + 1)} \leftarrow m^{(k)} + μ (x^{(k + 1)} - u^{(k + 1)}) \end{array}

In step (1), the conjugate gradient method (CG) is applied to the following linear system:

(A^{†} A + μ I) x = A^{†} b + μ (u^{(k)} - \frac{m^{(k)}}{μ}) .

The subproblem of step (2) is solved through a shrinkage operation:

u^{(k + 1)} = max {| x^{(k + 1)} + \frac{m^{(k)}}{μ} | - \frac{λ}{μ}, 0} sgn (x^{(k + 1)} + \frac{m^{(k)}}{μ}) .

Several parameters are required by the ADMM: number of sub-iterations of the CG, λ and μ. The more the conditioning number of A^†A differs from 1, the more sub-iterations are required by the CG. This number was set to 15. The parameter λ rules the tradeoff between fidelity and penalty term, for this reason, it should be selected depending on the amount of noise affecting the data. Suited values for λ and μ have been found through several experiments with simulated data. The reconstruction of simulated and real data are, in particular, performed with λ = μ = 1.0. It is important to point out that the trends presented in Sec.5.3 resulted experimentally indipendent from the choice of these two parameters.

References and links

1. H. Schomberg and J. Timmer, “The gridding method for image reconstruction by Fourier transformation,” IEEE Trans. Med. Imag. 14(3), 596–607 (1995). [CrossRef]

2. R. N. Bracewell and A. C. Riddle, “Inversion of fan-beam scans in radio astronomy,” Astrophys. J. 150, 427–434 (1967). [CrossRef]

3. W. N. Brouw, “Aperture synthesis,” Methods in Computational Physics14, 131–175 (Academic, 1975).

4. H. Sedarat and D. Nishimura, “On the optimality of the gridding reconstruction algorithm,” IEEE Trans. Med. Imag. 19(4), 306–317 (2000). [CrossRef]

5. P. J. Beatty, D. G. Nishimura, and J. M. Pauly, “Rapid gridding reconstruction with a minimal oversampling ratio,” IEEE Trans. Med. Imag. 24(6), 799–808 (2005). [CrossRef]

6. JD. O’Sullivan, “A fast sinc function gridding algorithm for Fourier inversion in computer tomography,” IEEE Trans. Med. Imag. 4(4), 200–207 (1985). [CrossRef]

7. J. I. Jackson, C. H. Meyer, D. G. Nishimura, and A. Macovski, “Selection of a convolution function for Fourier inversion using gridding (computerised tomography application),” IEEE Trans. Med. Imag. 10(3), 473–478 (1991). [CrossRef]

8. J. Walden, “Analysis of the direct Fourier method for computer tomography,” IEEE Trans. Med. Imag. 19(3), 211–222 (2000). [CrossRef]

9. S. Matej and I. G. Kazantsev, “Fourier-based reconstruction for fully 3D PET: optimization of interpolaton parameters,” IEEE Trans. Med. Imag. 25(7), 845–854 (2006). [CrossRef]

10. F. Marone and M. Stampanoni, “Regridding reconstruction algorithm for real-time tomographic imaging,” J. Synchrotron Radiat. 19, 1029–1037 (2012). [CrossRef] [PubMed]

11. J. Fessler and B. Sutton, “Nonuniform fast Fourier transforms using min-max interpolation,” IEEE Trans. Signal Process. 51(2), 560–574 (2003). [CrossRef]

12. S. Matej, J. A. Fessler, and I. G. Kazantsev, “Iterative tomographic image reconstruction using Fourier-based forward and back-projectors”, IEEE Trans. Med. Imag. 23(4), 401–412 (2004). [CrossRef]

13. F. Arcadu, M. Nilchian, A. Studer, M. Stampanoni, and F. Marone, “A forward regridding method with minimal oversampling for accurate and efficient iterative tomographic algorithms,” IEEE Trans. Imag. Process. 25(3), 1207–1218 (2016). [CrossRef]

14. S. R. Zhao and H. Halling, “A new Fourier method for fan-beam reconstruction,” Nucl. Scie. Symp. and Med. Imag. Conf. Record 2, 1287–1291 (1995).

15. Y. Zhang and J. A. Fessler, “Fourier-based forward and back-projectors in iterative fan-beam tomographic image reconstruction,” IEEE Trans. Med. Imag. 25(5), 582–589 (2006). [CrossRef]

16. S. Shaller, T. Flohr, and P. Steffen, “An efficient Fourier method for 3D Radon inversion in exact cone-beam CT reconstruction,” IEEE Trans. Med. Imag. 17(2), 244–250 (1998). [CrossRef]

17. P. Toft, “The Radon Transform: Theory and Implementation,” Ph.D. Thesis, (Denmark Technical University, 1996).

18. T. Weitkamp, A. Diaz, C. David, F. Pfeiffer, M. Stampanoni, P. Cloetens, and E. Ziegler, “X-ray phase imaging with a grating interferometer,” Opt. Express 13(16), 6269–6304 (2005). [CrossRef]

19. F. Pfeiffer, T. Weitkamp, O. Bunk, and C. David, “Phase retrieval and differential phase contrast imaging with low brilliance X-ray sources,” Nature Phys. 2, 258–261 (2006). [CrossRef]

20. A. Momose, S. Kawamoto, I. Koyama, Y. Hamaishi, K. Takai, and Y. Suzuki, “Demonstration of X-ray Talbot interferometry,” Jpn. J. Appl. Phys. 42(7B), 866–868 (2003). [CrossRef]

21. F. Pfeiffer, O. Bunk, P. Kraft, and E. F. Eikenberry, “High resolution brain tumor visualization using three-dimensional X-ray phase contrast tomography,” Phys. Med. Biol. 52(23), 6923–6930 (2007). [CrossRef] [PubMed]

22. T. Koehler, B. Brendel, and R. Proksa, “Iterative reconstruction for differential phase contrast imaging using spherically symmetric basis functions,” Med. Phys. 38(8), 4542–4545 (2011). [CrossRef]

23. S. Horbelt, M. Liebling, and M. Unser, “Discretization of the Radon transform and of its inverse by spline convolutions”, IEEE Trans. Med. Imag. 21(4), 363–376 (2002). [CrossRef]

24. Z. Qi, J. Zambelli, N. Bevins, and G. Chen, “A novel method to reduce data acquisition time in differential phase contrast computed tomography using compressed sensing,” Proc. SPIE7258, Medical Imaging (2009). [CrossRef]

25. M. Nilchian, C. Vonesch, P. Modregger, M. Stampanoni, and M. Unser, “Fast iterative reconstruction of differential phase contrast X-ray tomograms,” Opt. Express 21(5), 5511–5528 (2013). [CrossRef] [PubMed]

26. D. Hahn, P. Thibault, A. Fehringer, M. Bech, T. Koehler, F. Pfeiffer, and P. B. Noel, “Statistical iterative reconstruction algorithm for X-ray phase-contrast CT,” Sci. Rep. 5, 10452 (2015). [CrossRef] [PubMed]

27. F. Natterer and F. Wübbeling, Mathematical Methods in Image Reconstruction (SIAM, 2001). [CrossRef]

28. J. F. Kaiser, “Digital Filters,” in Systems Analysis by Digital Computer, (John Wiley and Sons, 1974) pp. 218–285.

29. H. B. Dwight, Table of Integrals and Other Mathematical Data (Macmillan Company, 1947).

30. R. C. Gonzalez and R. E. Woods, Digital Image Processing, 3rd Edition (Prentice Hall, 2008).

31. Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image Quality Assessment: From Error Visibility to Strutural Similarity,” IEEE Trans. Imag. Process. 13(4), 600–612 (2004). [CrossRef]

32. L. Shepp and B. F. Logan, “The Fourier Reconstruction of a Head Section,” IEEE Trans. Nucl. Sci. 21(3), 21–43 (1974). [CrossRef]

33. S. Boyd, N. Parikh, E. Chu, B. Peleato, and J. Eckstein, “Distributed optimization and statistical learning via the alternating direction method of multipliers,” Found. Trends Mach. Learn. 3(1), 1–122 (2011). [CrossRef]

(a)
	HF-DBRP	HF-DBP
SNR	27.37	24.69
CNR	9.99	9.18

(b)
	HF-DBRP	HF-DBP
SNR	24.83	23.57
CNR	8.22	8.19

(c)
	HF-DBRP	HF-DBP
SNR	22.65	22.24
CNR	7.23	7.15

	HF-DBRP	HF-DBP
805 views × 512 pix.	0.23 s	28.93 s
1608 views × 1024 pix.	0.97 s	2.40 · 10² s
3217 views × 2048 pix.	4.23 s	1.20 · 10³ s

(a)
	ADMM-DFRP	ADMM-DRT
PSNR	31.42	30.76
MSSIM	0.74	0.67

Fast gridding projectors for analytical and iterative tomographic reconstruction of differential phase contrast data

Abstract