Edge-illumination x-ray phase contrast imaging restoration using discrete curvelet regularization transform

Abstract

This article considers the problem of recovering edge-illumination x-ray phase contrast (EIXPC) images from a set of potentially Poisson noisy projection measurements. The authors cast a recovery as a sparse regularization problem based on Anscombe multiscale variance stabilizing transform (MS-VST) with fast discrete curvelet transform which was applied to simulated edge-illumination x-ray phase contrast images. For accurate modelling, the noise characteristics of the EIXPCi data are used to determine the relative importance of each projection. Two implementations of curvelet sparse regularization transforms were applied, including the unequally-spaced fast Fourier transform and the wrapping-based transform. The algorithms were evaluated in terms of contrast improvement, quality of image restoration, object perceptibility, and peak signal-to-noise ratio. The methods provide nearly optimal solution without excessive memory and recovery time requirement. The performance of the proposed algorithms is demonstrated through a series of complex numerical geometric and anthropomorphic phantom studies. The results of numerical simulations demonstrate that the discrete curvelet transform with MS-VST is fast and robust, and it can effectively improve image quality, preserve and enhance edges and restore lost information while significantly reducing the noise. Additionally, both sparse sampling and decreasing x-ray tube current (i.e. noisy data) lead to the reduction of radiation dose in the x-ray imaging.

Keywords

biomedical imaging curvelet regularization x-ray phase-contrast x-ray imaging

1 Introduction

Edge-illumination x-ray phase contrast imaging (EIXPCi) has great potential in the differentiation of small objects with increased sensitivity and enhanced image quality over absorption methods [1 –4] while significantly reducing the dose to the imaged specimen. Among numerous x-ray phase contrast imaging techniques, EIXPCi allows the use of conventional x-ray sources [5, 6]. However, EIXPCi modality requires the introduction of pre-sample and detector masks, so-called coded apertures, to achieve illumination of the detector pixel fraction based on the edge illumination principle [7]. Talbot interferometry also requires beam modifiers, such as gratings [8], and these additional optical components along the beam pathway increase the overall cost of such a system. The simplest is free-space propagation, since no extra optical elements are needed [9, 10]. In free-space propagation (FSP), an x-ray traversing the sample undergoes phase-shift, diffracts and generates interference fringes at the boundaries of the object. Phase information, which is essential for phase retrieval, is carried by these fringes [11]. However, the quality of the images decreases through the limitations imposed by the imaging system, which restrict precise identification of individual features and details.

In practical cases, reduced image quality is due to the finite size and partial coherence of the x-ray source, the limited spatial resolution of the detector and the inevitable noise of the imaging system components. We can identify two predominant sources of noise in the x-ray image acquisition: 1) the stochastic nature of the photon-counting process at the detectors; 2) the intrinsic thermal and electronic fluctuations of the acquisition devices. For the x-ray imaging, it is more reasonable to model the output of the detectors as a Poisson-distributed random vector.

Higher image contrast and quality assumes the use of a reduced source size (when compared to the source to sample distance) and superior imaging detector [12 –16]. However, such improvements require the utilization of additional and expensive hardware, which so far reveals physical limitations. Although Fresnel/Kirchhoff diffraction theory [17] with modifications provides the most effective ways to model the edge (partial) illumination x-ray phase contrast imaging, it requires high spatial resolution and coherence of the incident x-ray wave field which is impractical. In many medical applications, we are required to reconstruct an object from noisy and undersampled tomographic data [18] corrupted by the uncertainty of noise modeling in the measurements. The restoration methods can serve as a substitute to increase image quality when the hardware improvement is not achievable. Nevertheless, the unsuitability of numerous frequency and spatial domain algorithms for x-ray phase-contrast image quality enhancement corrupted by Poisson noise imposes a numerical challenge due to sampling artifacts, blurring of edges, or excessive computational time [19 –29]. These features mean that the results of such algorithms in EIXPCi applications are unsatisfactory. Wavelet regularization deconvolution transforms can deal with objects with point singularities but fail to deal with singularities along edges. Furthermore, wavelets exhibit high complexity and slow calculation time and are unsuitable for detecting or providing a compact representation of high dimensional geometric structures [30 –34]. Therefore, a fast and robust image restoration routine adapted to reconstruction problems with missing data, especially for EIXPCi, is needed.

This paper presents a study on the application of the discrete curvelet regularization transform (DCRT) algorithms [35 –37] with multiscale variance stabilizing transform (MS-VST) [29, 38], for the first time to the best of our knowledge, to the EIXPCi numerical system [39]. Our approach includes preprocessing the count data by a MS-VST which consists combining a VST with multiscale transform (wavelets or curvelets), yielding asymptotically normally distributed coefficients with known variances. The transform reforms the data so that the noise approximately becomes Gaussian with a constant variance [40].

In particular, both variants based on DCRT were implemented to the EIXPCi model, including the unequispaced fast Fourier transform (USFFT) and wrapping-based transform (WBT). Additionally, these transforms were applied to the FSP x-ray phase contrast imaging to test their usefulness in denoising of the FSP images. The advantage of curvelet regularization over other transforms lies in fact that DCRT constitutes a highly edge-preserving method, considers a multiscale time-frequency local partition and makes use of the direction of geometric features. USFFT and WBT were applied to simulated EIXPCi data and the corresponding results were compared to the fast Daubechies-based symlet 8 wavelet regularization technique.

2 Theory and methods

2.1 Poisson noise model and multiscale variance stabilizing transform

In many acquisition instruments i.e. biomedical detectors, computerized tomography, etc., the noise stems from fluctuations of a counting process [29]. In such a cases, the Poisson distribution represents the noise in the images obtained from the detection of particles, e.g., photons. The signal-to-noise ratio of the images is signal dependent and varies across the image plane. Each observation $Y [k] \sim P (X [k]), \forall k,$ where X [k] is underlying signal (so-called Poisson intensity) [40]. The problem is that now the variance of the noise is equal to its mean X [k], e.g. signal-dependent. Variance stabilizing transforms (VST) such as the Anscombe variance-stabilizing transform $A$ [38] offer a solution for Poisson noise removal. That method stabilizes the variance of the Poisson image to a constant, and the resulting image tends to homoscedastic Gaussian [29]. $Y_{s} [k] = A (Y [k]) = 2 \sqrt{Y [k] + 3 / 8} .$ (1)

After the Anscombe VST, we faced to the additive Gaussian noise. This makes it possible to employ Gaussian-based denoising methods. A drawback of the VST is that when the intensity level of the input image is very low (low signal-to-noise ratios), its performance deteriorates [40].

Multiscale variance-stabilizing transforms combine the VST with the low pass filters implemented in various multiscale multidirection transforms such as curvelets [35]. The low pass filters average out the noise to an extent thereby improving the signal-to-noise ratio [40]. In particular, multiscale multidirection transforms sparsely capture the intrinsic geometry of an image [29]. For Poisson image recovering using Multiscale-VST (MS-VST), we follow the steps described by Zhang et al. [40]. First, we apply the DCRT-based multiscale-VST for converting the image with Poisson noise to approximately Gaussian. Then we perform hard thresholding to detect the DCRT coefficients. The locations of the detected coefficients define a DCRT multiresolution support $M$ for the image [40]. The true intensity image is estimated by minimizing the ℓ₁ norm of the DCRT coefficients, reindexed as a column vector, with the constraint that the resulting image has non-negative pixel values and the coefficients belonging to $M$ are preserved. This is solved by using the Hybrid Steepest Descent (HSD)algorithm [41].

The Poisson curvelet denoising algorithm with MS-VST summary:

Apply the Multiscale-VST+wavelet with J scales to get the stabilized wavelet subbands w_j.

Set B₁ = B_min .

For j = 1 to J do:

Partition the subband w_j with blocks of side-length B_j .

Apply the digital ridgelet transform to each block to obtain the stabilized curvelet coefficients.

Detect the stabilized curvelet coefficients to obtain $M .$

If j mod 2 = 1 then

B_j+1 = 2B_j

else

B_j+1 = B_j

End if.

End for.

Apply the HSD iterations to the curvelet coefficients before getting the final estimate.

2.2 Noise level estimation

Gaussian noise level defined by noise variance ${\tilde{σ}}^{2}$ is estimated using a robust median estimator on the finest scale of the wavelet decomposition of the Fourier shrinkage estimate $\hat{f} (x)$ proposed byDonoho [42] $\tilde{σ} = \frac{median (| w_{J - 1, k}^{(n)} - median (w_{J - 1, k}^{(n)}) |)}{0.6745}$ where J = log ₂ (n) and the true wavelet coefficients w_j-1,k, k = 0,2^j-1-1 of the noisy data.

Such an estimation of noise level is based only on the empirical wavelet coefficients at the highest level of resolution, which consist mostly of noise – often, all but the finest details in the function, are accounted for in coefficients at lower levels.

For additive Gaussian noise ɛ, we have η_j,ℓ,k : = 〈 ɛ, φ_j,ℓ,k 〉, and then ${\tilde{σ}}_{j, ℓ} = \tilde{σ} ∥ φ_{j, ℓ, k} ∥,$

where φ_j,ℓ,k are the transform atoms. If the atoms in the dictionary all have equal unit ℓ₂ - norm, then ${\tilde{σ}}_{j, ℓ} = \tilde{σ}, \forall (j, ℓ) .$ This can be easily implemented for the case where the ℓ₂ - norm is known analytically, as is for the curvelet frame [29].

2.3 Discrete curvelet regularization transform theory

Implementation of discrete curvelet regularization transform to the EIXPCi numerical system is based on the method proposed by Candes et al. [37]. The curvelet regularization uses partitioning of the phase-space followed by the ridgelet transform [43 –45] applied to data precisely localized in terms of space and frequency. Strategies for the implementation of the edge-preserving unequispaced fast Fourier transform and wrapping-based transform provide the opportunity to improve the quality of EIXPCi signal with fast and less redundant methods. These two implementations differ in the way curvelets at a given scale and angle are translated with respect to each other [37]. In USFFT the translation grid is tilted to be aligned with the rotation of the curvelet, thus yielding the most precise discretization of the continuous definition. The wrapping-based version uses the same grid for all angles within each quadrant; hence, each curvelet is given the proper orientation. The curvelet regularizations were achieved by hard-thresholding of the curvelet coefficients at 3σ_j,l for all but the finest scale, where the threshold λ = 4σ_j,l, and σ_j,l is the noise level of the coefficient at scale j and angle l.

Both DCRT algorithms are based, in general, on linear transformations with input Cartesian arrays $f \in L^{2} (ℝ^{2}),$ where f [t₁, t₂], 0 < t₁, t₂ < n (constitute a directional frame), which has anisotropic needle-shaped elements and allows for an optimal sparse representation of objects with discontinuities along smooth curves. Output as curvelet coefficients c_j,l,k as integral over frequency plane can be defined as $c (j, l, k) : = \int_{ℝ^{2}} \bar{f (x) φ_{j, l, k} (x)} d x = \frac{1}{(2 π)^{2}} \int \hat{f} (ω) U_{j} (R_{θ_{l}} ω) e^{i 〈 x_{k}^{(j, l)}, ω 〉} d ω .$ (2)

Hence the discrete coefficients $c_{j, l, k}^{D}$ $c_{j, l, k}^{D} : = \sum_{0 \leq t_{1}, t_{2} < n} f [t_{1}, t_{2}] \bar{φ_{j, l, k}^{D} [t_{1}, t_{2}]},$ (3)

Curvelets are defined by three parameters as a function of spatial variable $x = (x_{1}, x_{2}) \in ℝ^{2}$ at scale $2^{- j}, j \in ℕ_{0}$ , an equispaced sequence of rotation angles θ_j,l = 2πl · 2^-⌊j/2⌋, with l = 0, 1, … such that 0 ≤ θ_l < 2π, 0 ≤ l ≤ 2^⌊j/2⌋ - 1 and position $x_{k}^{(j, l)} = R_{θ_{j, l}}^{- 1} (k_{1} \cdot 2^{- j}, k_{2} \cdot 2^{- j / 2}),$ $k = (k_{1}, k_{2}) \in ℝ^{2}$ are a sequence of translation parameters by $φ_{j, l, k} (x) = φ_{j} (R_{θ_{j, l}} (x - x_{k}^{(j, l)})),$ (4) where R_{θ
_j,l} denotes the rotation matrix with angle θ_j,l radians and $R_{θ_{j, l}}^{- 1} = R_{θ_{j, l}}^{T} = R_{- θ_{j, l}}$ its inverse and transpose, $R_{θ} = (\begin{matrix} cos θ_{j, l} & sin θ_{j, l} \\ - sin θ_{j, l} & cos θ_{j, l} \end{matrix})$ and φ_j are smooth functions with compact support on wedges in the Fourier domain. For such an object as f with discontinuities along C²- continuous curves, the best m-term approximation ${\tilde{f}}_{m}$ by curvelet thresholding obeys ${∥ f - {\tilde{f}}_{m} ∥}_{2}^{2} \leq {Cm}^{- 2} {(log m)}^{3},$ while for wavelets the decay rate is only m^-1 and for the Fourier transform the rate is m^-1/2. Unlike the isotropic elements of wavelet bases, the needle-shaped elements of this transform possess very high directional sensitivity and anisotropy. Such elements are very efficient in representing line-like edges; thus, curvelets are suitable for edge-preserving EIXPC image quality improvements.

Candes et al. [37] propose the family of curvelet functions which forms a tight frame of $L_{2} (ℝ^{2})$ hence each function $f \in L_{2} (ℝ^{2})$ has the representation $f = \sum_{j, l, k} 〈 f, ψ_{j, l, k} 〉 ψ_{j, l, k},$ (5) where 〈f, ψ_j,l,k〉 denotes L₂-scalar product of f and ψ_j,l,k . The coefficients 〈f, ψ_j,l,k〉 are curvelet coefficients of function f .

The frequency window is defined in the Fourier domain as $U_{j} (r, θ) = 2^{- 3 j / 4} W (2^{- j} r) V (\frac{2^{⌊ j / 2 ⌋} θ_{j, l}}{2 π})$ (6) where radial window W (r) and angular window V (t) satisfy admissibility conditions and are defined by $\sum_{j = - \infty}^{\infty} W^{2} (2^{j} r) = 1, r \in (3 / 4, 3 / 2);$ $\sum_{l = - \infty}^{\infty} V^{2} (t - l) = 1, t \in (- 1 / 2, 1 / 2);$ where j ≥ j₀ and ⌊j/2⌋ is the integer part of j/2 . The waveform φ_j (x) is defined by its Fourier transform ${\hat{φ}}_{j} (ω) = U_{j} (ω)$ and according to equation (6) U_j (ω) is the window defined in the polar coordinate system with ω a frequency domain variable.

Iterative curvelet thresholding [46] for signal u^k can be denoted as follows $u^{k + 1} = S_{σ, T} (u^{k} + Φ^{T} (f - Φ u^{k})) .$ where Φ corresponds to a low-pass filtering and ψ is a sparsity basis.

Here the shrinkage operator S_σ,T is given by $S_{σ, T} (f) = \sum_{j, l, k} T_{σ} (〈 f, ψ_{j, l, k} 〉) ψ_{j, l, k},$ where T_σ can be taken as a soft shrinkage function defined by a fixed threshold σ > 0, $T_{s, σ} (x) = {\begin{matrix} x - σ, & x \geq σ, \\ 0, & | x | < σ, \\ x + σ, & x \leq - σ, \end{matrix} ∥$ or a hard shrinkage function $T_{h, σ} (x) = {\begin{matrix} x, & | x | \geq σ, \\ 0, & | x | < σ . \end{matrix}$

In general, simple thresholding will not find sparse decompositions when good decompositions exist. A curvelet transform allows the restoration of sparsity by reducing redundancy across scales. In detail, interscale orthogonality can be introduced by means of subband filtering. Thus, different levels of the multiscale ridgelet pyramid are used to represent different subbands of a filter output. Hence, subband decomposition imposes the anisotropy scaling relationship between the width and effective length of the important frame elements width = length².

The translation is direct. However, in digital applications improved numerical and visual results can be obtained by leaving the subbands separate, applying spatial partitioning to each subband and applying the ridgelet transform on each subband separately. This version of the curvelet transform is redundant with redundancy factor equal to 16j + 1 with j scales [29]. Two presented methods, USFFT and Wrapping-based use the same coronization, which were explained in detail by Candes et al. [37].

A. Unequispaced Fast Fourier Transform

The USFFT algorithm allows the curvelet coefficients to be found by irregularly sampling the Fourier coefficients of an image. The USFFT requires about 2^3j/2 coefficients per scale and angle and O (n²) storage and achieves an accuracy of the order of 10^–6 with low running times.

The localizing window for a Fourier transform can be defined as ${\tilde{U}}_{j} (ω) : = ψ_{j} (ω_{1}) V_{j} (ω),$ (7) with concentric tiling where ψ_j (ω₁) = ψ (2^-jω₁) and $ψ (ω_{1}) = \sqrt{φ (ω_{1} / 2)^{2} - φ (ω_{1})^{2}}$ for each θ_l ∈ [- π/4, π/4] ${\tilde{U}}_{j, l} : = ψ_{j} (ω_{1}) V_{j} (S_{θ_{l}} ω) = {\tilde{U}}_{j} (S_{θ_{l}} ω),$ which assures shear-invariant windows at any given scale. Shear matrix S_θ is defined as $S_{θ_{l}} : = (\begin{matrix} 1 & 0 \\ - tan θ & 1 \end{matrix})$

Angles θ_l are unequispaced and ${\tilde{U}}_{j, l}$ define the Cartesian version of U_j (R_{θ
_l}ω) and the Cartesian curvelets ${\tilde{φ}}_{j, l, k} (x) = 2^{3 j / 4} {\tilde{φ}}_{j} (S_{θ_{l}}^{T} (x - S_{θ_{l}}^{- T} b))$ where b : = (k₁ · 2^-j, k₂ · 2^-j/2) takes the discrete values.

Then, curvelet coefficients are given by $c (j, l, k) = \int \hat{f} (ω) {\tilde{U}}_{j} (S_{θ_{l}}^{- 1} ω) e^{i 〈 S_{θ_{l}}^{- T} b, ω 〉} d ω$ (8)

To find the digital version of the coefficients c (j, l, k), it is convenient to pass the shearing operation to $\hat{f}$ then $c (j, l, k) = \int \hat{f} (ω) {\tilde{U}}_{j} (S_{θ_{l}}^{- 1} ω) e^{i 〈 b, S_{θ_{l}}^{- 1} ω 〉} d ω = \int \hat{f} (S_{θ_{l}} ω) {\tilde{U}}_{j} (ω) e^{i 〈 b, ω 〉} d ω$ (9)

Then, the Equation (9) can be solved numerically after the shearing operation to recover a rectangular grid. The two-dimensional discrete Fourier transform $\hat{f} [n_{1}, n_{2}], - n / 2 \leq n_{1}, n_{2} < n / 2$ of a given Cartesian array f [t₁, t₂], 0 ≤ t₁, t₂ < n can be defined by $\hat{f} [n_{1}, n_{2}] = \sum_{t_{1}, t_{2} = 0}^{n - 1} f [t_{1}, t_{2}] e^{- i 2 π (n_{1} t_{1} + n_{2} t_{2}) / n},$

Additionally, non-integer values n₁, n₂ are allowed and some interpolation is necessary. Samples $\hat{f} [n_{1}, n_{2}] = \hat{f} (2 π n_{1}, 2 π n_{2})$ from the interpolating trigonometric polynomial are defined by $\hat{f} (ω_{1}, ω_{2}) = \sum_{0 \leq t_{1}, t_{2} < n} f [t_{1}, t_{2}] e^{- i (ω_{1} t_{1} + ω_{2} t_{2}) / n} .$

Assume ${\tilde{U}}_{j} [n_{1}, n_{2}]$ as supported on a rectangle of length L_1,j and width L_2,j $P_{j} = {(n_{1}, n_{2}) : n_{1, 0} \leq n_{1} < n_{1, 0} + L_{1, j}, n_{2, 0} \leq n_{2} < n_{2, 0} + L_{2, j}},$ (10) where (n_1,0, n_2,0) is the index of the pixel at the bottom-left of the rectangle.

Then, discrete curvelet coefficients are evaluated through unequispaced fast Fourier transform $c_{j, l, k}^{D} = \sum_{n_{1}, n_{2} \in P_{j}} \hat{f} [n_{1}, n_{2} - n_{1} tan θ_{l}] {\tilde{U}}_{j} [n_{1}, n_{2}] e^{i 2 π (k_{1} n_{1} / L_{1, j} + k_{2} n_{2} / L_{2, j})},$ (11)

To sum up, the USFFT algorithm consists of the following steps:

Obtain Fourier samples $\hat{f} [n_{1}, n_{2}], - n / 2 \leq n_{1}, n_{2} < n / 2$ through applying 2D FFT

Obtain sampled values $\hat{f} [n_{1}, n_{2} - n_{1} tan θ_{l}]$ for (n₁, n₂) ∈ P_j, for each scale/angle pair (j, l) resample or interpolate $\hat{f} [n_{1}, n_{2}] .$

Multiply parabolic window ${\tilde{U}}_{j} [n_{1}, n_{2}]$ with the sheared object $\hat{f},$ localizing $\hat{f}$ near the parallelogram with orientation θ_l, and obtain ${\tilde{f}}_{j, l} [n_{1}, n_{2}] = \hat{f} [n_{1}, n_{2} - n_{1} tan θ_{l}] {\tilde{U}}_{j} [n_{1}, n_{2}] .$

Collect the discrete regularization coefficients $c_{j, l, k}^{D}$ for each 2D FFT applied to ${\tilde{f}}_{j, l}$ .

B. Wrapping-based Transform

The Wrapping-based transform using a series of translations and a wrap-around technique gives a faster computation time compared to the USFFT. The choice of the spatial grid to translate curvelets at each scale and angle in the wrapping-based approach assumes using a rectangular grid instead of a tilted grid with a similar Cartesian curvelets definition $c (j, l, k) = \int \hat{f} (ω) {\tilde{U}}_{j} (S_{θ_{l}}^{- 1} ω) e^{i 〈 b, ω 〉} d ω$ (12) where ${bS}_{θ_{l}}^{- T}$ in Equation (9) has been replaced by b ≃ (k₁2^-j, k₂2^-j/2), taking on values on a rectangular grid.

However, in the frequency plane, the window ${\tilde{U}}_{j, l} [n_{1}, n_{2}]$ does not fit in a rectangle of size ∼2^j × 2^j/2, aligned with the axes, in which the 2D inverse FFT could be applied to solve Equation (12). Discretization of the integral over ω leads to a sum over n₁, n₂ which extends beyond the bounds allowed by the 2D IFFT. In principle, it is possible to use the 2D inverse FFT on a larger rectangle instead, which is similar to the discretization of the continuous directional wavelet transform [47]. Unfortunately, such an approach has an important limitation which is the very high oversampling of the coefficients. For fine scale curvelets, a wrapping-based transform requires to computation at the order of 2^2j coefficients per scale and angle and O (n^2.5) storage, which is more than the USFFT method.

Downsampling the grid for each angle and scale allows a subgrid to be obtained which has the same cardinality as in the unequispaced Fourier transform. Downsampling can be performed through periodization of the frequency samples. For each scale j two constants exist L_1,j ∼ 2^j and L_2,j ∼ 2^j/2 where, for every orientation θ_l, it is possible to tile the 2D plane with translates of P_j,l (see Equation (11)) by multiples of L_1,j in the horizontal direction and L_2,j in vertical direction. Hence, periodization of the windowed data $d [n_{1}, n_{2}] = {\tilde{U}}_{j, l} [n_{1}, n_{2}] \hat{f} [n_{1}, n_{2}]$ $Wd [n_{1}, n_{2}] = \sum_{m_{1} \in ℤ} \sum_{m_{2} \in ℝ} d [n_{1} + m_{1} L_{1, j}, n_{2} + m_{2} L_{2, j}]$

Then the wrapped data must be defined as the restriction of Wd [n₁, n₂] to indices n₁, n₂ inside a rectangle with sides of length L_1,j × L_2,j near the origin: 0 ≤ n₁ < L_1,j and 0 ≤ n₂ < L_2,j . Thus, the wrapping method is a reindexing of the data. The Wrapping-based algorithm summary:

Obtain Fourier samples $\hat{f} [n_{1}, n_{2}], - n / 2 \leq n_{1}, n_{2} < n / 2$ through applying 2D FFT to the image f

Obtain the product ${\tilde{U}}_{j, l} [n_{1}, n_{2}] \hat{f} [n_{1}, n_{2}]$ for each scale j and angle ℓ (divide FFT into collection of digital corona tiles)

Wrap ${\tilde{U}}_{j, l} [n_{1}, n_{2}] \hat{f} [n_{1}, n_{2}]$ around the origin and obtain

{\tilde{f}}_{j, l} [n_{1}, n_{2}] = W ({\tilde{U}}_{j, l} \tilde{f}) [n_{1}, n_{2}]

Collect the discrete regularization coefficients $c_{j, l, k}^{D}$ for each 2D inverse FFT applied to ${\tilde{f}}_{j, l}$ .

The shift-invariant block thresholding of the wrapping-based curvelet transform with curvelets at the finest scale was used as an additional strategy to improve the output of the restoration algorithm.

2.4 Wavelet transform theory

In general, signals with singularities such as images can be represented using wavelet regularization. The wavelet transform was implemented with reconstruction filters using the orthogonal wavelet family. Orthogonal wavelets are used to develop the discrete wavelet transform (DWT), which returns a vector of the same length as the input [48]. A large amount of data in such a vector is equal to zero. This corresponds to the fact that it decomposes the data into a set of wavelets that are orthogonal to its translations and scaling. Therefore, it is possible to decompose a signal to the same or lower number of wavelet coefficient spectra as the number of signal data points. The discrete wavelets are used in multi-resolution analysis constituting an orthonormal basis for $L^{2} (ℝ)$ [48, 49]. The Daubechies-based symlet 8 wavelet regularization was introduced to the EIXPCi model to provide reconstruction using a well-known fast algorithm for a comparison with discrete curvelet regularizationmethods.

A. Daubechies-based wavelet regularization

Implementation of wavelet regularization to the EIXPCi simulation was based on the well-known Daubechies wavelets [50], which involve a scaling function φ and a wavelet function ψ . In general, signal x (t) decomposed into low and high frequency components, named as detail and approximation coefficients respectively can be reconstructed for orthonormal basis $L^{2} (ℝ)$ as $x (t) = \sum_{m = 1}^{L} [\sum_{k = - \infty}^{\infty} δ_{m} (k) ψ_{j, k} (t) + \sum_{k = - \infty}^{\infty} α_{l} (k) φ_{j, k} (t)]$ (13) where ψ_j,k (t) is the discrete wavelet and φ_j,k (t) is discrete scaling, δ_m (k) is the detailed signal at scale 2^m and α_l (k) is the approximated signal at scale 2^l . The parameters δ_m (k) and α_l (k) are obtained using the scaling and wavelet filters [51].

As proposed by Donoho and Johnstone [49] thresholding was applied to the given signal by implementing hard and soft thresholding, the so-called shrinkage method. In hard thresholding, the wavelet coefficients below a given value are set to zero; while in soft thresholding, coefficients are reduced by a quantity equal to the threshold value, which is the estimation of the noise level calculated from the standard deviation of the detail coefficient. Analytically, it can be expressed as $HT : {\begin{matrix} y = x, & | x | > λ \\ y = 0, & | x | < λ \end{matrix}$ $ST : {y = sgn (x) (| x | - λ)$ where x is the input signal, y is the signal after threshold and $λ = σ \sqrt{2 log (n)}$ is the threshold value. Although more sophisticated methods have been proposed to overcome the drawbacks of the non-linear methods, soft and hard thresholding are still the most efficient and reliable [42]. The image recovery algorithm using wavelets regularization consists of the following steps:

The discrete wavelet transform should be applied to the noisy dataset y (i) = (y_n) _n∈[1,N] and y (i) is the sum of the noise-free signal f (i) = (f_n) _n∈[1,N] and the noise b = (b_n) _n∈[1,N] denoted as $y (i) = f (i) + b (i), b (i) = σ ɛ (i), i = 0, 1, \dots, n - 1$ where ɛ (i) are independently normal random variables and σ represents the intensity of the noise in y (i),

Computing J estimates ${\hat{f}}_{j}$ of the noise-free high-pass subbands f_j to restore J wavelet subimages $y^{j} = f^{j} + b^{j} = (y_{n}^{j})_{n \in [1, N]}, j \in [1, J]$ ,

Applying the inverse discrete wavelet transform on the high-pass wavelet subimages ${\hat{f}}_{j}$ to obtain an estimate $\hat{f}$ of the noise-free data f, and hence recover a denoised image.

2.5 Numerical system and performance metrics

The data acquisition system is based on the numerical model of the edge illumination x-ray phase-contrast imaging technique including monochromatic x-ray point source of 30 keV energy (for free-space propagation case) and tungsten target x-ray tube source of energy range of 10–100 keV with spot size of 75 μm, pre- and post-sample masks and a flat panel detector. Detailed description of the model and its parameters can be found in [39].

The EIXPC imaging was performed using 60 keV x-ray tube source energy. A conventional x-ray tube source generating diverging beam was applied [39]. The pre-sample mask pitch and aperture were 66.8 μm and 12 μm; for the detector mask, they were 83.5 μm and 20 μm, obtained in an approximately 30 μm thick gold layer. The source-to-sample distance (SSD) of 160 cm with a corresponding sample-to-detector distance (SDD) of 40 cm. However, for free-space propagation SDD and SSD were both switched to 80 cm. In the presented experiment the images were acquired with detector pixel pitch of 85 μm. The detector is designed to incorporate a scalable pixel array. However, images were scaled down to a 1024×1024 pixel array resolution.

The simulation serves as an important tool for evaluating, examining and optimizing latest EIXPCi technology; it provides practical predictions of results without difficulties due to uncontrolled experimental variables. Differential edge-illumination phase images were obtained using phase retrieval method implemented to the numerical simulation [52]. Modeling of the imaging process enables the acquisition of the projections of a numerical phantom of human anatomy as if it were a live patient; hence, it is feasible to execute pre-clinical studies entirely in silico. The advantage of such examinations is that the anatomy of the object (phantom) is known and only simulation provides a “ground truth” from which to quantitatively assess and refine imaging equipment and methods.

Simple geometric wire phantoms, soft tissue bubble and anthropomorphic complex 4D XCAT 2 phantom [53] imitating the human head region were introduced to the imaging system. These virtual phantoms served for comparison of the results of adopted curvelet and wavelet regularization methods for x-ray phase contrast image restoration.

Here, the used x-ray propagation model was based on Monte Carlo code and wave-optics techniques, where the emission and history of each ray were calculated in parallel using graphics processing unit (GPU), implemented with Nvidia CUDA. The synergy of GPU combined with CUDA technology was applied to optimize the level of simulation performance. The calculated multiple values were grouped in the process within the GPU, reducing the number of memory reads [39].

The analysis and restoration procedures of edge illumination x-ray phase contrast images were conducted on a highly specialized workstation, equipped with dual Intel Xeon E5 2.6 GHz processors (12 cores), 32 GB RAM and GPU Nvidia Quadro 600 graphics.

The numerical examples of image quality improvement methods consist of a sparsity analysis with a quantitative study of recovery of a given signal f for the images obtained from the model. The peak signal-to-noise ratio (PSNR) of the best approximation is $PSNR = 10 {log}_{10} (\frac{max (f^{2})}{〈 {| \hat{f} - f |}^{2} 〉}) (dB)$ (14) where $\hat{f}$ is the estimate of the noise-free image and the mean-squared error (MSE) in the space domain is a weighted sum of the MSE of each individual subband defined as $MSE = 〈 {| \hat{f} - f |}^{2} 〉 = \sum_{j = 0}^{J} \frac{N_{j}}{N} 〈 {| {\hat{f}}_{j} - f_{j} |}^{2} 〉$ (15)

For 8-bit images, this is usually max(f²) =255² . The aim is to maximize the PSNR and to minimize MSE (see Equation (16). It is worth noting that the advantage of simulated image is that it is the only case we have a true value of f.

The Human Visual System (HVS) is mainly adapted to extract the structure of objects from the images, which is the principle behind the SSIM [54]. Consequently, a measure of the structural similarity constitutes fine approximation of perceptual image quality. SSIM tries to ignore non-structural distortions, such as local intensity patterns and is defined in the spatial domain as $SSIM (x, y) = \frac{(2 μ_{x} μ_{y} + C_{1}) (2 σ_{xy} + C_{2})}{(μ_{x}^{2} + μ_{y}^{2} + C_{1}) (σ_{x}^{2} + σ_{y}^{2} + C_{2})},$ (16) where C₁ and C₂ are two small positive constants; μ_x and μ_y are the mean of images x and y, respectively; σ_x and σ_y are the standard deviation of images x and y, correspondingly; and σ_xy is simply $σ_{xy} = \frac{1}{n} \sum_{i = 1}^{n} (x_{i} - μ_{x}) (y_{i} - μ_{y}) .$ (17)

The maximum value SSIM can achieve is 1, and it is only achieved when x and y are identical. The value of sigma for the gaussian filter in the local window to avoid blocking artefacts is 1.5 with Gaussian window type and filter width 11 pixels. The values of C₁ = 0.01 and C₂ = 0.03 were set to avoid instability in the calculation of SSIM index [54].

Contrast values presented in this work were computed as $C = \frac{I_{max} - I_{min}}{I_{background}},$ (18) where I_max is the maximum intensity of the positive fringes, is I_min the minimum intensity of the negative ones and I_background is the background intensity outside the detail. This is an adaptation to x-ray images of the classic Born and Wolf [17] formula for fringe visibility.

3 Results

The free-space propagation x-ray phase-contrast was simulated first to test proposed methods in simple setup without apertures [9]. Figure 1(a) gives the result obtained by imaging a simple pseudo-organic phantom containing a relatively weakly absorbing wire of 465 μm diameter. Such a sample allowed observation of positive and negative peaks in the signal; thus, dark fringes were not suppressed by absorption [55]. An image of the wire phantom was obtained with SSD = 80 cm and corresponding SDD = 80 cm. The simulation was performed for an x-ray point source of 30 keV; the noisy images have the PSNR = 20 dB. Figure 1(b) shows Poisson noisy image. Figure 1(c) and (d) show symlet 8 wavelet regularization with soft thresholding, decomposition level 4, maximum sparsity and a wavelet symlet 8 with hard thresholding, respectively. The results were obtained with corresponding PSNR_sym8,ST = 31.60 dB and PSNR_sym8,HT = 31.52 dB separately. Figure 1(e) presents the result obtained with the USFFT method. PSNR_USFFT in Fig. 1(e) is calculated as 35.94 dB. The WBT result is shown in Fig. 1(f), and a PSNR_WBT of 36.07 dB was obtained. It should be noted that the image contrast rises from the original 0.589 visible in Fig. 1(a) to 0.606 in Fig. 1(f). To illustrate the improved visibility of the Fresnel fringes and the edge enhancement, the right edge of the wire in the middle of Fig. 1(a)–(f) is shown on an expanded scale in Fig. 2(a)–(f). Figure 3 presents the profile comparison in the displacement direction. Figure 4(a)–(d) illustrate SSIMs for recovered images and values are 0.9839, 0.9677, 0.9903, and 0.9901 for Fig. 1(c)–(f), respectively.

To further investigate the effect of the mentioned algorithms on phase contrast image restoration, another simulation of a noisy signal was introduced for recovery. A wire-phantom consisting of four wires of 450 μm, 350 μm, 585 μm and 100 μm diameter, respectively, was used. The wires were represented by different atomic Z numbers and absorption properties. SSD = 160 cm and corresponding SDD = 40 cm were set in the experiment and the numerical imaging system was switched to work in an edge-illumination x-ray phase contrast imaging regime. The offset of a post-sample mask was set to achieve 66% illumination of the detector pixels. Symlet 8 wavelet regularization with hard and soft thresholding, 5th level of decomposition with high sparsity, and the unequispaced fast Fourier transform and wrapping based curvelet regularization were separately performed on these noisy differential phase maps, and mixed absorption and phase contrast images. The PSNR of the noisy images was kept constant at 20 dB. Figure 5(c)–(f) show the visual comparison results between these techniques. Figure 5(d) was restored using DWT with symlet 8 wavelets and a strategy based on the 2D stationary wavelet transform (SWT) [56]. The basic idea was to average many slightly different discrete wavelet analyzes to obtain better results, since DWT is not a translation invariant transform. To illustrate the improved visibility of the edges and the quality of the image in that particular region of interest, the edges of two wires (middle and right) of Fig. 5(a)–(f) are shown on an expanded scale in Fig. 6(a)–(f). Furthermore, Fig. 7(a)–(d) present the profile comparison for implemented image restoration methods with an original and noisy image. After symlet 8 regularization (with soft and hard thresholding), the corresponding PSNRs are 23.59 dB (SSIM = 0.9184) and 24.07 dB (SSIM = 0.9363), respectively. The USFFT method obtained a significantly higher PSNR_USFTT = 30.03 dB than the symlet 8 techniques with SSIM = 0.9165. The wrapping-based transform still obtained the best level of the peak signal-to-noise ratio PSNR_WBT = 31.01 dB (SSIM = 0.9598). The wrapping-based regularization and the USFFT method result in similar PSNR values; however, the WBT algorithm was twice as fast as the USFFT. For differential phase maps, the wrapping-based curvelet transform obtained the highest contrast value of 1.26, while simultaneously preserving the edges and avoiding excessive blurring (see Fig. 7(d)). Structural artifacts corresponding to the Fig. 5(c)–(f) are visible on Fig. 8(a)–(d).

Figures 9(c)–(f) demonstrate the visual comparison results between applied transforms. Figure 9(d) was restored using DWT with symlet 8 wavelets and SWT. In addition, the edges of two wires (middle and right) of Fig. 9(a)–(f) are presented on an expanded scale in Fig. 10(a)–(f). Moreover, Fig. 11(a)–(d) show the profile comparison for restoration methods with ground truth and noisy image. After symlet 8 regularization (with soft and hard thresholding), the corresponding PSNRs are 30.91 dB (SSIM = 0.9532) and 31.86 dB (SSIM = 0.9612), respectively. The USFFT method obtained higher PSNR_USFTT = 38.62 dB than the symlet 8 techniques with SSIM = 0.9438. The wrapping-based transform again obtained the best level of the peak signal-to-noise ratio PSNR_WBT = 38.94 dB (SSIM = 0.9591). Among these methods, the wrapping-based curvelet transform obtained the highest contrast value of 1.34, with better edge preserving and without excessive blurring (see Fig. 11(d)). Structural artifacts corresponding to the Fig. 9(c)–(f) are also visible on Fig. 10(a)–(d).

A two-dimensional bubble sample imitating soft tissue was introduced to the EIXPCi simulation to perform test of the unequispaced fast Fourier algorithm and wrapping method based on the discrete curvelet transform. Soft tissue is a relatively weakly absorbing material with an atomic Z≈ 7.4. The simulation was performed for an x-ray tube mean energy of 60 keV. The offset of a post-sample mask was set to achieve only a 20% illumination of the detector pixel. Figure 13(a) shows the original differential phase image. Figure 13(c) and (d) separately present the Wavelet regularization with symlet 8 soft and hard thresholding results. The corresponding PSNRs are 25.41 dB and 25.82 dB. Figure 13(e) shows the result of the Wavelet regularization with symlet 8 hard thresholding with Stationary Wavelet Transform (SWT). PSNR in Fig. 13(e) is calculated as 25.51 dB. The USFFT result is shown in Fig. 13(f), and a PSNR of 32.94 dB is obtained. Figure 13(g) presents the Wrapping-based result with a corresponding PSNR of 33.31 dB. Figure 14 presents the profile comparison for implemented image restoration methods. Structural similarity index maps are shown on Fig. 15(a)–(e) with corresponding SSIM values 0.9573, 0.9477, 0.9601, 0.9402, and 0.9592 respectively.

The EIXPCi model was subsequently switched to work with complex samples; and anthropomorphic head phantom 4D XCAT 2 was imaged. The phantom contains high complexity pseudo biological structures and inhomogeneities with different attenuation properties, which makes it feasible to test image restoration techniques in preclinical conditions. The 4D XCAT 2 phantom is based on non-uniform rational B-spline shapes of mixed absorbing properties with a wide range of atomic Z numbers. The phantom was illuminated by polychromatic 60 keV x-rays. The images of differential phase maps and mixed phase and absorption contrast restoration are presented separately in Figs. 16 and 17 for better clarity. Original image of head phantom is shown in Figs. 16(a) and 17(a). Figure 16(b) shows the Poisson noisy differential phase image and Fig. 17(b) presents the Poisson noisy mixed phase and absorption contrast image. Figure 16(c) and (d) separately present the Wavelet regularization with symlet 8 soft and hard thresholding results. The Stationary Wavelet transform was applied only to the hard thresholding result; the corresponding PSNRs are 25.64 dB and 25.31 dB. The USFFT result is shown in Fig. 16(e), and a PSNR of 28.67 dB is obtained. Figure 16(f) presents the Wrapping-based result with a corresponding PSNR of 31.74 dB. Figure 17(c) and (d) separately show the results of the Wavelet symlet 8 with soft and hard thresholding; the corresponding PSNRs are 32.78 dB and 28.19 dB. For mixed phase and absorption contrast the USSFT result is shown in Fig. 17(e) and Wrapping-based recovery image is presented in Fig. 17(f). The curvelet transforms result PSNRs are 31.96 dB and 36.37 dB respectively. The Structural Similarity Index maps for Fig. 16(c)–(f) are presented in Fig. 18(a)–(d) with corresponding values 0.7971, 0.7993, 0.7918, and 0.8914. The SSIMs for Fig. 17(c)–(f) are shown in Fig. 19(a)–(d) with SSIM values 0.8369, 0.8378, 0.8411, and 0.9248. Contrast values were not calculated for these cases.

4 Discussion

In practical situations, to perform 1D or 2D deconvolution for an experimental x-ray phase-contrast imaging system, the line-spread function (LSF) and full point-spread function (PSF) should be obtained respectively. In case of good isotropic properties of the x-ray source and the detector, the PSF could be obtained as an approximation by rotating the LSF of the system [34, 57]. The PSF could be used for deconvolution procedures to improve image quality and denoise images obtained from the x-ray phase contrast imaging modalities. The deconvolution procedure can partially remove negative effects of finite source size and limited detector resolution. Significant improvements on the hardware detectors and conventional x-ray sources still require a long period of research. The development of such equipment allows an improvement in the quality of the output of the imaging system, but simultaneously leads to sharply increasing costs.

The software implementation of the image restoration technique for differential x-ray phase contrast based on Wiener deconvolution as suggested by Olivo and Speller [58] was still not satisfactory. In a real experimental implementation of the EIXPCi technique, system noise cannot be excluded and direct deconvolution leads to low contrast. The MS-VST + Wavelet regularization methods allow better quality images to be obtained and can deal with ringing artifacts, such as the results presented for SWT, but these methods decrease contrast (profiles in Fig. 7(a) and 11(a) reveal up to a 90% reduction of positive and negative peaks in the signal for weakly absorbing wire) and cause blurring of the edges (see Figs. 1, 5, 9, 13, 16 and 17); thus, a more appropriate image restoration technique is required. In this paper, we mainly presented implementation of a Multiscale Variance Stabilizing Unequispaced fast Fourier transform and MS-VST+Wrapping-based transform to the EIXPCi simulation; corresponding results were compared with the MS-VST+Wavelet symlet 8 regularization (soft and hard thresholding) and Stationary Wavelet Transform technique.

The proposed USFFT and Wrapping-based algorithms were evaluated with simulated data in terms of peak signal-to-noise ratio, perceptibility as SSIMs and phase contrast improvement. Original noisy phase contrast profiles of a 465 μm thick wire were analyzed for free-space propagation case. The phase contrast images of pseudo-organic wires of 450 μm, 350 μm, 585 μm and 100 μm diameter, soft tissue bubble and complex human head 4D XCAT 2 phantoms were used as practical data for the preclinical test cases in edge-illumination x-ray phase contrast. The Wavelet symlet 8 algorithm with soft thresholding and SWT allowed noisy projections to be restored and better visual image quality to be obtained, but suppressed positive and negative peaks in the signal (see Figs. 3, 7, 11, and 14) due to excessive blurring of the edges (see Figs. 1(c), 5(c), 9(c), 13(c), 13(d), 16(c), 16(d), 17(c), and 17(d)). The Wavelet regularization with hard thresholding provides higher contrast and is a better edge preserving method than SWT; nonetheless, the ability to distinguish between certain materials diminishes in the presence of ringing artifacts, and the perceptibility (SSIM) of the details is decreased (see Figs. 16(b) and 17(b)).

Among the tested methods, the Wrapping-based transform results in the highest PSNR, SSIM and best contrast improvement (see Figs. 1(f), 5(f), 9(f), 13(g), 18(d) and 19(d)). The WBT technique allows a slightly higher PSNR and SSIM to be obtained than the USFFT algorithm; however, the Wrapping-based transform exhibits significantly shorter computation times. The comparison of the SSIM images resulting from filtering the phase contrast images for both curvelet and wavelet transforms reveals the presence of different features. The curvelets result in high fidelity in the restoration of the curvy edges, while wavelets perform well only at point singularities (see Figs. 18 and 19). The superiority of curvelet regularization in noise suppression upon the wavelet transform mainly originates from the curvelet shrinkage and produces a Mean Squared Error (MSE) order of magnitude better than what is achieved by other methods. Strictly speaking, the recovery is demonstrably asymptotically near-optimal.

In summary, the implemented sophisticated MS-VST+Curvelet based regularization algorithm outperformed the wavelet techniques in terms of the following advantages: 1) It is efficient and sensitive in detecting and preserving the edges and curvy features in the presence of Poisson noise; 2) The simplicity of the MS-VST allows the computation time to be comparable to that of a Gaussian denoising, which makes that method capable of processing large data sets. The investigated algorithms can serve as an effective approach to quantitative edge-illumination x-ray phase contrast imaging. To obtain better image recovery additional techniques could be used, i.e. cycle spinning for translation invariant wavelet regularization or l-1 iterations [35]. Orthogonal wavelet transforms are not shift-invariant (translation-invariant), which means that the processing of an original and translated image leads to different results. The cycle spinning allows the recovery algorithm to turn into a shift-invariant version when the transformation is applied to several shifted copies of the image and the resulting recovered images are shifted back to the original position, and the results are averaged.

The main advantages of WBT and USFFT methods are the high level of noise suppression and the preservation of the phase contrast edges, while wavelets are efficient in isotropic feature analysis. The DCRT method is also successful to bring out the high frequency details, but contains a few ringing artifacts, which is due to the lesser redundancy and the inherent limitation of that technique in recovering applications. Combining the DCRT with other methods can overcome these drawbacks greatly [59], as well as applying image denoising methods based on effective filtering in 3D transform domain by combining sliding-window transform processing with block matching [60] which is now under investigation. However, both methods suffer from pseudo-Gibbs oscillating artifacts in the discontinuities (see Figs. 3(e), 3(f), 4, 6(d), 6(e), 7, 10, and 11) and the recovery artifacts are visible. In addition, depending on the noise level present in the input images, with increasing values of ${\tilde{σ}}^{2}$ we observed decreasing PSNR values and for noise level standard deviation >50, the PSNR values are very low (<10 dB). The algorithms demonstrated excessive blurring of the edges where visibility of low contrast objects was deteriorated, which for medical imaging leads to the lost of diagnostic information. The redundancy of the curvelet transform implementation is about 2.8 when wavelets are chosen at the finest scale, and 7.2 otherwise. However, curvelets possess an infinite number of directional moments. This property implies that, if the essential support of the curvelet φ_j,k,l lies in a smooth part of f, then the corresponding curvelet coefficient c (j, k, l) will be small, while, if the essential support of φ_j,k,l is aligned with an edge of f, then c(j,k,l) will be significant [37]. Furthermore, the noise in the EIXPCi, depends on the noise in the input images, and on the characteristics of the experimental setup. The real imaging system define how much the refraction of x-ray traversing the object is “amplified” to provide the detected signal. In particular, the noise in the refraction image depends upon many parameters, including the fraction of the beam impinging on the detector mask (illumination level), the propagation distance and the value of the beam distribution at the two edges of the aperture [61].

The EIXPCi technique with MS-VST+Wrapping-based curvelet regularization could be switched to work as a 3D tomography system with incomplete, undersampled data [62] allowing further dose reduction to the imaged specimen while maintaining the quality of reconstruction, which is highly important in the biological imaging. Such an implementation could be highly useful in preclinical testing of the EIXPCi system and is the focus of our future work. Additionally, the curvelet-based transform may perform even better as a deconvolution procedure with known approximation of the PSF of real imaging system. It is worth noting that the presented algorithms could be easily adapted to work with any differential x-ray phase-contrast system, i.e. FSP, EIXPCi or grating interferometry, significantly improving the quality of the output.

5 Conclusions

A curvelet regularization transform with MS-VST using unequispaced fast Fourier method and wrapping-based algorithm has been developed for EIXPC image recovery with Poisson noisy low-dose projection data. Presented techniques allow higher quality image of the specimen, by eliminating the step of mapping EIXPCi data to the Fourier domain. The performance of the proposed method is successfully demonstrated by simulation phantom studies. It is shown that wrapping-based curvelet transform outperforms the well-known wavelet algorithm when dealing with sparse data in the presence of relatively high Poisson noise of the imaging system and reduced imaging quality due to the limited spatial resolution of the detector, finite size and partial coherence of the x-ray source. The results indicate that it is possible to reduce radiation dose without loss of useful information for diagnostics or image-guided radiation therapy (IGRT).

References

Arfelli

, et al., Mammography with synchrotron radiation: Phase–detection techniques, Radiology 215 (2000), 286–293.

and Liu

, Clinical implementation of x-ray phase-contrast imaging: Theoretical foundations and design considerations, Medical Physics 30 (2003), 2169–2179.

Lewis

R.A.

, Medical phase contrast x-ray imaging: Current status and future prospects, Physics in Medicine and Biology 49 (2004), 3573–3583.

Bravin

, et al., X-ray phase-contrast imaging: From pre-clinical applications towards clinics, Physics in Medicine and Biology 58 (2013), R1.

Olivo

and Speller

, Experimental validation of a simple model capable of predicting the phase contrast imaging capabilities of any x-ray imaging system, Physics in Medicine and Biology 51 (2006), 3015–3030.

Olivo

and Speller

, A coded-aperture technique allowing x-ray phase contrast imaging with conventional sources, Applied Physics Letters 91 (2007), 074106.

Munro

P.R.

, et al., A quantitative, non-interferometric X-ray phase contrast imaging technique, Optics Express 21 (2013), 647–661.

Pfeiffer

, Weitkamp

, Bunk

and David

, Phase retrieval and differential phase-contrast imaging with low-brilliance x-ray sources, Nature Physics 2 (2006), 258–261.

Wilkins

, Gureyev

, Gao

, Pogany

and Stevenson

, Phase contrast imaging using polychromatic hard x-ray, Nature 384 (1996), 335–338.

10.

Ghani

, et al., Low dose high energy x-ray in-line phase sensitive imaging prototype: Investigation of optimal geometric conditions and design parameters, Journal of X-ray Science and Technology 23 (2015), 667–682.

11.

and Liu

, X-Ray cone-beam phase tomography formulas based on phase-attenuation duality, Optics Express 13 (2005), 6000–6014.

12.

Nesterets

Y.I.

, Wilkins

S.W.

, Gureyev

T.E.

, Pogany

and Stevenson

A.W.

, On the optimization of experimental parameters for x-ray in-line phase-contrast imaging, Review of Scientific Instruments 76 (2005), 093706.

13.

Toth

, Kieffer

J.C.

, Fourmaux

, Ozaki

and Krol

, In-line phase-contrast imaging with a laser-based hard x-ray source, Review of Scientific Instruments 76 (2005), 083701.

14.

Wong

M.D.

, et al., Image quality and dose efficiency of high energy phase sensitive x-ray imaging: Phantom studies, Journal of X-ray Science and Technology 22 (2014), 321–334.

15.

, Liu

and Yan

, Optimization of X-ray phase-contrast imaging based on in-line holography, Nuclear Instruments and Methods in Physics Research Section B 234 (2005), 563–572.

16.

Chen

, Zheng

, Ai-Min

and Cheng-Quan

, Influence of tube voltage and current on in-line phase contrast imaging using a microfocus x-ray source, Chinese Physics 16 (2007), 2319–2324.

17.

Born

, Wolf

, Principles of Optics, London, Pergamon, 1980.

18.

Natterer

, The mathematics of computerized tomography, Teubner

B.G.

, John Wiley & Sons, 1986.

19.

Zhou

Y.T.

, Chellappa

, Vaid

and Jenkins

B.K.

, Image restoration using a neural network, IEEE Transactions on Acoustics, Speech, and Signal Processing 36 (1988), 1141–1151.

20.

Paik

J.K.

and Katsaggelos

A.K.

, Image restoration using a modified Hopfield network, IEEE Transactions on Image Processing 1 (1992), 49–63.

21.

Rudin

, Osher

and Fatemi

, Nonlinear total variation based noise removal algorithm, Physica D 60 (1992), 259–268.

22.

and Kundu

, Image estimation using fast modified reduced update Kalman filter, IEEE Transactions on Signal Processing 40 (1992), 915–926.

23.

Chambolle

and Lions

, Image recovery via total variation minimization and related problems, Numerische Mathematik 76 (1997), 167–188.

24.

Caselles

, Morel

J.M.

, Tannenbaum

and Sapiro

, Introduction to the special issue on partial differential equations and geometry-driven diffusion in image processing and analysis, IEEE Transactions on Image Processing 7 (1998), 269–273.

25.

Chen

, Chen

and Zhou

, Adaptively regularized constrained total least-square image restoration, IEEE Transactions on Image Processing 9 (2000), 589–596.

26.

Karl

W.C.

, Regularization in image restoration and reconstruction in Handbook of Image and Video Processing, Second Edition, Elsevier, 2005.

27.

Chan

, Esedoglu

, Park

and Yip

, Recent development of total variation image restoration in Handbook of Mathematical Models of Computer Vision, Springer Verlag, 2005.

28.

Kim

, PDE-based image restoration: A hybrid model and color image denoising, IEEE Transactions on Image Processing 15 (2006), 1163–1170.

29.

Starck

J.L.

, Murtagh

and Bijaoui

, Sparse Image and Signal Processing: Wavelets and Related Geometric Multiscale Analysis, Cambridge University Press, Cambridge, 2015.

30.

Banham

M.R.

, Galatsanos

N.P.

, Gonzalez

H.L.

and Katsaggelos

A.K.

, Multichannel restoration of single channel images using a wavelet-based subband decomposition, IEEE Transactions on Image Processing 3 (1994), 821–833.

31.

Lee

and Paik

, Space-frequency adaptive image restoration based on wavelet decomposition, in Proceedings of IEEE Asia Pacific Conference on Circuits and Systems, Institute of Electrical and Electronics Engineers, Seoul, 1996.

32.

Sendur

and Selesnick

, Bivariate shrinkage functions for wavelet-based denoising exploiting interscale dependency, IEEE Transactions on Signal Processing 50(11) (2002), 2744–2756.

33.

Neelamani

, Choi

and Baraniuk

, ForWaRD: Fourier-wavelet regularized deconvolution for illconditioned systems, IEEE Transactions on Image Processing 52 (2004), 418–433.

34.

Zhou

, Gao

, Zhao

and Zhang

, Application of Fourier-wavelet regularized deconvolution for improving image quality of free space propagation x-ray phase contrast imaging, Physics in Medicine and Biology 57 (2012), 7459–7479.

35.

Candes

E.J.

and Donoho

D.L.

, New tight frames of curvelets and optimal representations of objects with piecewise C2 singularities, Communications on Pure and Applied Mathematics 57 (2003), 219–266.

36.

Candes

E.J.

and Donoho

D.L.

, Continuous curvelet transform. II. Discretization and Frames, Applied and Computational Harmonic Analysis 19 (2005), 198–222.

37.

Candes

E.J.

, Demanet

, Donoho

D.L.

and Ying

, Fast discrete curvelet transforms, Multiscale Modeling and Simulation 5 (2006), 861–899.

38.

Anscombe

F.J.

, The transformation of Poisson, binomial and negative-binomial data, Biometrika 35 (1948), 246–254.

39.

Włodarczyk

and Pietrzak

, Analytical reconstructions of intensity modulated x-ray phase-contrast imaging of human scale phantoms, Biomedical Optics Express 6(11) (2015), 4255–4272.

40.

Zhang

, Fadili

J.M.

and Starck

J.L.

, Wavelets, ridgelets, and curvelets for Poisson noise removal, IEEE Transactions on Image Processing 17 (2008), 1093–1108.

41.

Yamada

, The hybrid steepest descent method for the variational inequality problem over the intersection of fixed point sets of nonexpansive mappings in Inherently parallel algorithm in deasibility and optimization and their applications, Elsevier, New York, 2001.

42.

Donoho

D.L.

, Denoising by soft-thresholding, IEEE Transactions on Information Theory 41 (1995), 613–627.

43.

Candes

E.J.

, Harmonic analysis of neural networks, Applied and Computational Harmonic Analysis 6 (1999), 197–218.

44.

Candes

E.J.

and Donoho

D.L.

, Ridgelets: The key to higher-dimensional intermittency?, Philosophical Transactions of The Royal Society A 357 (1999), 2495–2509.

45.

Candes

E.J.

and Donoho

D.L.

, Curvelets – a surprisingly effective nonadaptive representation for objects with edges, edition Rabut

, Cohen

and Schumaker

L.L.

, Curves and Surfaces, Nashville, TN: Vanderbilt University Press, 2000, pp. 105–120.

46.

and Le

F.X.

, Dimet, Deblurring from highly incomplete measurements for remote sensing, IEEE Transactions on Geoscience and Remote Sensing 47 (2009), 792–802.

47.

Vandergheynst

and Gobbers

J.F.

, Directional dyadic wavelet transforms: Design and algorithms, IEEE Transactions on Image Processing 11 (2002), 363–372.

48.

Daubechies

, Orthonormal basis of compactly supported wavelets, Communications on Pure and Applied Mathematics 41 (1998), 909–996.

49.

Donoho

D.L.

and Johnstone

J.M.

, Ideal spatial adaptation by wavelet shrinkage, Biometrika 81 (1994), 425.

50.

Vetterli

and Herley

, Wavelets and filter banks: Theory and design, IEEE Transactions on Signal Processing 40 (1992), 2207–2232.

51.

Mallat

S.G.

, A wavelet tour of signal processing, Academic Press, 1999.

52.

Munro

P.R.

, et al., Phase and absorption retrieval using incoherent X-ray sources, Proceedings of the National Academy of Sciences 109 (2012), 13922–13927.

53.

Segars

W.P.

, Sturgeon

, Mendonca

, Grimes

and Tsui

B.M.

, 4D XCAT phantom for multimodality imaging research, Medical Physics 37 (2010), 4902–4915.

54.

Wang

, Bovik

A.C.

, Sheikh

H.R.

and Simoncelli

E.P.

, Image quality assessment: From error visibility to structural similarity, IEEE Transactions on Image Processing 13 (2004), 600–612.

55.

Hwu

, et al., Coherence-enhanced synchrotron radiology: Simple theory and practical applications, Journal of Physics D: Applied Physics 35 (2002), R105–R120.

56.

Coifman

R.R.

and Donoho

D.L.

, Translation invariant de-noising, Lecture Notes in Statistics 103 (1995), 125–150.

57.

Samei

, Flynn

M.J.

and Reimann

D.A.

, A method for measuring the presampled MTF of digital radiographic systems using an edge test device, Medical Physics 25 (1998), 102–113.

58.

Olivo

and Speller

, Deconvolution of x-ray phase contrast images as a way to retrieve phase information lost due to insufficient resolution, Physics in Medicine and Biology 54 (2009), N347–N354.

59.

and Plonka

, The curvelet transform, IEEE Signal Processing Magazine 27 (2010), 861–899.

60.

Makitalo

and Foi

, Optimal inverse of the Anscombe transformation in low-count Poisson image denoising, IEEE Transactions on Image Processing 20 (2011), 99–109.

61.

Diemoz

P.C.

, et al., Edge illumination X-ray phase-contrast imaging: Nanoradian sensitivity at synchrotrons and translation to conventional sources, Journal of Physics: Conference Series 499 (2014), 012006.

62.

Ying

, Demanet

and Candes

E.J.

, 3D discrete curvelet transform, Optics & Photonics (2005), 591413–591413–11.