Expand this Topic clickable element to expand a topic
Skip to content
Optica Publishing Group

Suppressing multi-material and streak artifacts with an accelerated 3D iterative image reconstruction algorithm for in-line X-ray phase-contrast computed tomography

Open Access Open Access

Abstract

In-line X-ray phase-contrast computed tomography typically contains two independent procedures: phase retrieval and computed tomography reconstruction, in which multi-material and streak artifacts are two important problems. To address these problems simultaneously, an accelerated 3D iterative image reconstruction algorithm is proposed. It merges the above-mentioned two procedures into one step, and establishes the data fidelity term in raw projection domain while introducing 3D total variation regularization term in image domain. Specifically, a transport-of-intensity equation (TIE)-based phase retrieval method is updated alternately for different areas of the multi-material sample. Simulation and experimental results validate the effectiveness and efficiency of the proposed algorithm.

© 2022 Optica Publishing Group under the terms of the Optica Open Access Publishing Agreement

1. Introduction

X-ray phase-contrast imaging (PCI) is a promising non-invasive imaging technique. In contrast to conventional absorption-based X-ray imaging, PCI is based on alterations of the phase that arises from the refraction of X-rays passing through the objects to generate image contrast. The interaction between X-rays and the object can be described with its complex refractive index n = 1-δ+iβ, where δ is the refractive index decrement responsible for the phase information and β is the absorption index responsible for the absorption information. For hard X-ray imaging of materials with low density and low-Z (e.g., biological soft tissues), the variations of δ are approximately 3 orders of magnitude larger than those of β, which enables PCI to provide the superior contrast performance at a lower absorbed dose compared to conventional absorption-based X-ray imaging [1]. Due to its high sensitivity and its potential to provide superior image contrast between materials with similar attenuation properties, PCI attracts much attention for decades [2,3]. Over the past two decades, a variety of PCI techniques have been developed to transform phase information into measurable intensity variations on the detector, such as crystal interferometry imaging [4], analyzer-based imaging [5], in-line PCI (IL-PCI) [6], and grating-based imaging [7]. Among these techniques, IL-PCI is the simplest to implement in practice because it omits any additional apparatus in the beam path [8]. By extending IL-PCI to computed tomography (CT), it develops the IL-PCCT, which can provide the high-resolution three-dimensional (3D) visualization of weakly absorbing objects [9]. Nowadays, the applications of IL-PCCT mainly consist of two aspects: 1) it can be used as a high-resolution non-invasive microscopic imaging technique for biomedical and pre-clinical researches [10,11]; 2) it can serve as a 3D virtual histology method of human soft tissues for clinical studies [12,13]. In general, IL-PCCT requires 3D image reconstruction from a series of two-dimensional (2D) phase-contrast projections acquired at various tomographic angles. The reconstruction procedure for IL-PCCT usually consists of two independent inverse problems: phase retrieval and CT reconstruction [14]. Phase retrieval aims to extract the phase-shift information from the phase-contrast projections. CT reconstruction aims to yield the 3D refractive index decrement δ distribution from phase-shift information.

In principle, multiple phase-contrast projections, which are acquired at different sample-to-detector distances (SDDs) or at different energies for same tomographic angle, can be used to accurately retrieve the phase-shift information for multi-material or heterogeneous objects [15,16]. However, these methods are difficult to arrange experimentally and pose a high radiation dose to the biological tissues [17]. Therefore, it is very attractive to retrieve the phase-shift information from a single phase-contrast projection. Currently, many single-distance phase retrieval algorithms have been developed by imposing different restrictions onto the samples [18,19]. Among them, the homogeneous transport of intensity equation (TIE-Hom) method is one of the most widely used phase retrieval algorithms, which imposes an assumption on the homogeneity of the object directly in the TIE model [20]. The so-called homogeneity assumption is that the object consists of only a single material and thus the attenuating and phase-shifting properties of the sample are coupled, meaning that the δ and β are in a certain ratio, i.e., γ=δ/β, where γ is a constant. For many samples containing multiple materials, the assumption holds reasonably well. However, for some multi-material biological samples in which each material has a large difference in γ values (e.g., biological samples containing both weakly and strongly absorbing tissues), due to their violation of the homogeneity assumption, only the phase-shift information of the single-material area can be correctly recovered using TIE-Hom. Thus, artifacts will appear locally or as a whole in CT images. On one hand, when the γ value is suitable for the area of weak absorbing material, it will work well with them using TIE-Hom while causing blurring in the area of strong absorbing material. On the other hand, when the γ value is suitable for the area of strong absorbing material, it will work well with them, while the phase-contrast fringes will appear in the area of weak absorbing material. In this study, the above two types of artifacts are called multi-material artifacts. To solve the multi-material artifact problem, Beltran et al. [21,22] and Ullherr et al. [23] designed single-material phase retrieval method for multi-material samples containing three or more materials. However, the multiple reconstructions that come from each pair of materials are required to be spliced together manually, which is typically difficult for many multi-material biological samples in practice due to the complex textures and edge details.

The risk of X-ray radiation dose to biological samples has always been a concern. To decrease the radiation dose to an acceptable level, reducing the number of projections (e.g., sparse-view projections) is a feasible strategy for IL-PCCT. However, it will inevitably increase the data inconsistency, which can introduce streak artifacts into CT images. Compared with analytic reconstruction (AR) algorithms (such as filtered back-projection (FBP)), iterative reconstruction (IR) algorithms can model the geometric prior of data acquisition and strengthen physical constraints (such as total variation (TV)), and therefore can provide better reconstruction performance based on sparse-view projections, which attracted great attention for low-dose CT reconstruction [24]. In the field of IL-PCCT, many IR algorithms have been proposed to reduce the X-ray radiation dose [2529], which can promote the applications of IL-PCCT for imaging biological samples. Among these algorithms, an iterative FBP (IFBP) algorithm [28,29] was developed for the low-dose IL-PCCT reconstruction of single-material biological soft tissues, in which an advanced filtering technique and a regularization technique based on local feature priors were exploited to further improve CT reconstruction models. In addition, a deep learning-based sparse-view CT reconstruction framework for IL-PCCT was also developed, which could be used to decrease the X-ray radiation dose for single-material biological soft tissues [30]. However, these newly developed reconstruction algorithms are two-step approaches for IL-PCCT, and they mainly focused on the reduction of radiation dose in the CT reconstruction procedure and did not consider the multi-material artifacts problem in the phase retrieval procedure.

Typically, phase retrieval and CT reconstruction in IL-PCCT are performed independently, and this often ignores the influence of phase retrieval on CT reconstruction, such as multi-material artifacts. The combination of phase retrieval and CT reconstruction can help to design a more accurate IL-PCCT reconstruction algorithm, and the combined algorithms have made great progress. In 2013, Kostenko et al. [31] first proposed a combined algorithm by coupling the phase retrieval method based on the contrast transfer function (CTF) and the algebraic reconstruction algorithms, and their results demonstrated that the one-step approach was superior to the conventional two-step approach. Later, Tian et al. [32] designed a forward model by combining TIE with a discretized projection operator in the Fourier domain. Additionally, they also developed an iterative reconstruction method to obtain the phase information of weakly absorbing objects. In 2014, Ruhlandt et al. [33] used a Gerchberg-Saxton-type algorithm with a consistency constraint coupled with an algebraic reconstruction technique (ART) [34] to reconstruct the complex refractive index. In 2020, Hehn et al. [35] developed a novel IL-PCCT reconstruction model by combining a TIE-based phase retrieval method and a statistical iterative reconstruction technique into one step, which included models for the X-ray source and the detector, and it was important for low-dose IL-PCCT reconstruction based on laboratory X-ray sources. However, these algorithms did not consider multi-material samples that could not meet the homogeneity assumption, and therefore they cannot address multi-material artifacts in IL-PCCT reconstruction images.

Inspired by the abovementioned work, to suppress multi-material artifacts and streak artifacts simultaneously, we propose an accelerated 3D iterative image reconstruction algorithm, which mainly consists of three processes. First, TIE-Hom method is updated alternately for different areas of the multi-material sample to obtain the corresponding phase shifts and, thus correcting the multi-material artifacts. Then, IFBP algorithm combined with 3D TV regularization is used to reconstruct CT images and suppress streak artifacts. Finally, CT images in the previous process is forward projected to generate phase-contrast projections, and the data fidelity term based on the Euclidean distance between the generated phase-contrast projections and original phase-contrast projections is built in raw projection domain, which is utilized to correct reconstruction errors. In the above processes, the phase retrieval and CT reconstruction are merged into one step, and this can help to solve multi-material artifacts problem and streak artifacts problem simultaneously. To evaluate the performance of the proposed algorithm, simulations on a 3D microvasculature phantom and IL-PCCT experiments on a rat maxilla sample were performed.

The contributions of this study are summarized as follows:

  • 1. We propose a novel one-step approach that can jointly solve multi-material artifacts problem in phase retrieval step and streak artifacts problem in CT reconstruction step, and additionally, it can provide superior reconstruction performance for IL-PCCT compared with the conventional two-step approaches.
  • 2. In the proposed algorithm, the IFBP is utilized to reduce the computational cost and accelerate the CT reconstruction speed, which can significantly improve the reconstruction efficiency of the one-step approach.
  • 3. Considering that the proposed algorithm can serve as an effective tool to suppress multi-material artifacts and streak artifacts simultaneously, it has potential to promote the imaging application of IL-PCCT for multi-material biological samples.

The rest of this paper is organized as follows. In Section 2, we will briefly introduce the forward model of IL-PCCT in 2.1 and conventional two-step algorithms in 2.2, and then, the proposed one-step algorithm will be explained in 2.3 and 2.4. In Sections 3 and Sections 4, the results and analysis of simulation and real IL-PCCT experiments are provided, respectively. Finally, conclusion and discussion are drawn in Section 5.

2. Methods

2.1. Forward model of IL-PCCT

In IL-PCCT, an object illuminated by spatially coherent X-rays can be described with its 3D complex refractive index: $n(x,y,z) = 1 - \delta (x,y,z) + \textrm{i}\beta (x,y,z)$, where $\delta (x,y,z)$ is the refractive index decrement, $\beta (x,y,z)$ is the absorption index, and (x, y, z) denote spatial coordinates. The wave-object interaction can then be represented as the object transmission function, and it is defined as follows:

$${T_\theta }(x,y) = \exp [{ - {B_\theta }({x,y} )\textrm{ + }i{\varphi_\theta }(x,y)} ],$$
$${B_\theta }({x,y} )= k\int {\beta ({{\mathbf r}_ \bot },z)dz} ,$$
$${\varphi _\theta }(x,y) ={-} k\int {\delta ({{\mathbf r}_ \bot },z)dz} ,$$
where $\theta$ is the projection angle, ${T_\theta }\textrm{(}x,y)$ is the object transmittance function, ${{\mathbf r}_ \bot } = (x,y)$ denote the coordinates of the projection plane, z describes the propagation axis, ${B_\theta } = (x,y)$ denotes the absorption function, ${\varphi _\theta } = (x,y)$ represents the phase-shift function, k = 2π/λ is the wavenumber of the X-ray, and λ is the wavelength of the X-ray. Based on the absorption function and phase-shift function, the measured intensity recorded by the detector at a propagation distance z = D can be modeled as the squared modulus of the exit wave. The recorded intensity $I_\theta ^D$ can be described as follows:
$$I_\theta ^D({{\mathbf r}_ \bot }) = |{T_\theta }(x,y) \ast {P_\theta }^D({{\mathbf r}_ \bot }){|^2},$$
where * is the convolution operator.${P_\theta }^D({{\mathbf r}_ \bot })$ represents the Fresnel propagator, and its expression can be written as follows:
$${P_\theta }^D({{\mathbf r}_ \bot }) = \frac{{{e^{i2\pi D/\lambda }}}}{{i\lambda D}}{e^{i\frac{\pi }{{\lambda D}}{{|{{{\mathbf r}_ \bot }} |}^2}}}.$$

2.2. Conventional two-step algorithms for IL-PCCT

2.2.1. Phase retrieval

TIE-Hom can be performed to recover phase-shift information from a single detected intensity map by solving the TIE. The TIE can be written in the following form:

$${\nabla _ \bot } \cdot (I_\theta ^D({{\mathbf r}_ \bot }){\nabla _ \bot }{\varphi _\theta }({{\mathbf r}_ \bot })) ={-} k\frac{\partial }{{\partial z}}I_\theta ^D({{\mathbf r}_ \bot }),$$
where ${\nabla _ \bot } = (\frac{\partial }{{\partial x}},\frac{\partial }{{\partial y}})$ is the gradient operator in the plane ${{\mathbf r}_ \bot } = (x,y)$.

Assuming that a homogeneous sample is illuminated with parallel X-rays of approximately uniform intensity, its phase-shift information ${\varphi _\theta }({{\mathbf r}_ \bot })$ can be obtained by the TIE-Hom [20], and the formula can be written as follows:

$${\varphi _\theta }({{\mathbf r}_ \bot }) = \frac{1}{2}\gamma \ln \left[ {{{\cal F}^{ - 1}}\left( {\frac{{{\cal F}(I_\theta^D({{\mathbf r}_ \bot }))}}{{\pi \gamma \lambda D({\mathbf k}_ \bot^2) + 1}}} \right)} \right].$$
where $\gamma$ represents the ratio between δ(x, y, z) and β(x, y, z), ${{\mathbf k}_ \bot }$ denotes the coordinates of ${{\mathbf r}_ \bot }$ in the Fourier domain, ${\cal F}$ and ${{\cal F}^{ - 1}}$ are the 2D Fourier transform and its inverse operators, respectively.

2.2.2. CT reconstruction

After the phase-shift information ${\varphi _\theta }$ is retrieved by using the TIE-Hom, the δ information can be obtained by applying the standard FBP algorithm, and it can be written as follows:

$$F ={-} {k^{ - 1}}\int_0^\pi {{\varphi _\theta } \ast vd} \theta ,$$
where * is the convolution operator, ${\varphi _\theta }$ is the phase-shift information obtained by the TIE-Hom, v defines the ramp filter, and F denotes the 3D δ distribution of the illuminated object, i.e., F=δ (x, y, z).

The forward projection model of CT imaging can be approximately represented as the following discrete linear system, which is described as follows:

$$AF = \varphi ,$$
where A is the system matrix that represents the X-ray parallel beam forward projection, and $\varphi$ is the phase-shift function.

To solve the discrete linear system represented in Eq. (9), the optimization-based strategy is often adopted. In this study, Eq. (9) is transformed into a unconstrained optimization problem, and the objective function of the CT reconstruction task can be formulated as follows:

$$F = \arg \mathop {\min }\limits_F \frac{1}{2}\parallel (AF - \varphi )\parallel _2^2 + \; \alpha R(F),$$
where $\frac{1}{2}\parallel (AF - \varphi )\parallel _2^2$ denotes the data fidelity term, R (…) denotes the regularization term, and $\alpha$ is the regularization parameter that controls the trade-off between the data fidelity term and the regularization term.

2.3. Proposed one-step algorithm for IL-PCCT

Conventional two-step algorithms cannot correct the influence of phase retrieval on CT reconstruction because they mainly establish the data fidelity term in the phase-shift domain. In this study, the data fidelity term is built in the raw projection domain to reconstruct the δ distribution more accurately, and an additional constraint is added to the underdetermined linear system to accurately reconstruct F. Therefore, the objective function of the IL-PCCT reconstruction task can be rewritten as the following formula:

$${F^ \ast } = \arg \mathop {\min }\limits_F \frac{1}{2}\parallel \sum\limits_{\theta \textrm{ = }0}^\pi {H_\theta ^D(P_\theta ^D(AF) - I_\theta ^D)} \parallel _2^2 + \; \alpha R(F),$$
where $P_\theta ^D({\cdot} )$ denotes the Fresnel propagation (refer to Eq. (4)), $H_\theta ^D({\cdot} )$ denotes the phase retrieval operator for TIE-Hom (refer to Eq. (7)), and ${F^ \ast }$ represents the optimal solution to the problem in Eq. (9). Considering that the area of the same material in sample has the piecewise smoothness property on CT images, 3D TV is adopted as the regularization term in this study. Its expression is as follows:
$${R_{3DTV}}(F) = \sum\limits_{i,j,k} {\sqrt {{{({\nabla _1}F)}^2} + {{({\nabla _2}F)}^2} + {{({\nabla _3}F)}^2}} } ,$$
where ${\nabla _1}F = {F_{i,j,k}} - {F_{i - 1,j,k}}$, ${\nabla _2}F = {F_{i,j,k}} - {F_{i,j - 1,k}}$ and ${\nabla _3}F = {F_{i,j,k}} - {F_{i,j,k - 1}}$ are the gradients in the horizontal, vertical, and CT axial directions, respectively. Therefore, the objective function of the proposed one-step algorithm for IL-PCCT can be expressed as:
$${F^ \ast } = \arg \mathop {\min }\limits_F \frac{1}{2}\parallel \sum\limits_{\theta \textrm{ = }0}^\pi {H_\theta ^D(P_\theta ^D(AF) - I_\theta ^D)} \parallel _2^2 + \; \alpha {R_{3DTV}}(F).$$

2.4. Implementation of the proposed algorithm

According to the alternating direction method of multipliers (ADMM) [36], Eq. (13) can be transformed into two subproblems: the data fidelity term and the TV regularization term. The optimal solution obtained by alternate iterations of two subproblems is equivalent to the solution of the original problem.

2.4.1. Iterative CT reconstruction

The data fidelity term can be solved using the algebraic reconstruction algorithms (e.g., simultaneous algebraic reconstruction technique (SART). Considering the insufficient projection data problem caused by sparse-view projections, the IR algorithm is adopted in this study. SART is one of the most-widely used IR algorithms, and its formula can be expressed as follows:

$$F_j^t = F_j^{t - 1} + \frac{{{\lambda _t}}}{{\sum\limits_{i = 1}^N {{a_{i,j}}} }}\sum\limits_{i = 1}^N {\frac{{{H^D}(I_i^D) - \sum\limits_{n = 1}^M {{a_{i,n}}F_n^{t - 1}} }}{{\sum\limits_{n = 1}^M {a_{i,n}^{}} }}{a_{i,j}}} ,j = 1,2,3,\ldots ,M,$$
where t denotes the iterative number, ${a_{i,j}}$ denotes the contribution of the j-th pixel to the i-th X-ray, and ${\lambda _t}$ is a relaxation parameter in the t-th iteration.

SART can provide the high-quality reconstruction performance under sparse-view projections condition, but it has high computational cost and slow reconstruction speed [37]). The Landweber is a simple IR algorithm to solve the data fidelity term in Eq. (10) [38], and the IFBP is the improved version of the Landweber, which uses FBP as a backprojection operator to replace the ${A^T}$ in the Landweber formula and thus can significantly reduce computational cost and improve reconstruction speed [28,29]. Thus, the IFBP is used an accelerated iterative reconstruction strategy in this study. The IFBP can be expressed as follows:

$${F^t} = {F^{t - 1}} - \rho \; \mathrm{{\cal R}}(\sum\limits_{\theta \textrm{ = }0}^\pi {H_\theta ^D(P_\theta ^D(A{F^{t - 1}})} - I_\theta ^D)),$$
where ρ is the step size, $\mathrm{{\cal R}}$ represents the FBP operator, as expressed in Eq. (8).

2.4.2. 3D TV regularization

To solve the 3D TV regularization problem, the gradient descent method is adopted [39,40]. The formula is expressed as follows:

$$F_{3DTV}^{t + 1} = F_{IFBP}^t - {\lambda ^\ast }\frac{{\partial {R_{3DTV}}(F)}}{{\partial F}},$$
where ${\lambda ^\ast }$ is a parameter that can control the gradient descent.$F_{3DTV}^{t + 1}$ and $F_{IFBP}^t$ represent the (t+1)-th and t-th regularization results, respectively. The $\frac{{\partial {R_{3DTV}}(F)}}{{\partial {F_{i,j,k}}}}$ can be calculated as follows:
$$\begin{aligned} \; \frac{{\partial {R_{3DTV}}(F)}}{{\partial {F_{i,j,k}}}} &\approx \frac{{({{F_{i,j,k}} - {F_{i - 1,j,k}}} )+ ({{F_{i,j,k}} - {F_{i,j - 1,k}}} )+ ({{F_{i,j,k}} - {F_{i,j,k - 1}}} )}}{{\sqrt {\xi + {{({{F_{i,j,k}} - {F_{i - 1,j,k}}} )}^2} + {{({{F_{i,j,k}} - {F_{i,j - 1,k}}} )}^2} + {{({{F_{i,j,k}} - {F_{i,j,k - 1}}} )}^2}} }}\\&\textrm{ }\; - \frac{{({{F_{i + 1,j,k}} - {F_{i,j,k}}} )}}{{\sqrt {\xi + {{({{F_{i + 1,j,k}} - {F_{i,j,k}}} )}^2} + {{({{F_{i + 1,j,k}} - {F_{i + 1,j - 1,k}}} )}^2} + {{({{F_{i + 1,j,k}} - {F_{i + 1,j,k - 1}}} )}^2}} }}\\&\textrm{ }\; - \frac{{({{F_{i,j + 1,k}} - {F_{i,j,k}}} )}}{{\sqrt {\xi + {{({{F_{i,j + 1,k}} - {F_{i - 1,j + 1,k}}} )}^2} + {{({{F_{i,j + 1,k}} - {F_{i,j,k}}} )}^2} + {{({{F_{i,j + 1,k}} - {F_{i,j + 1,k - 1}}} )}^2}} }}\\&\textrm{ }\; - \frac{{({{F_{i,j,k + 1}} - {F_{i,j,k}}} )}}{{\sqrt {\xi + {{({{F_{i,j,k + 1}} - {F_{i - 1,j,k + 1}}} )}^2} + {{({{F_{i,j,k + 1}} - {F_{i,j - 1,k + 1}}} )}^2} + {{({{F_{i,j,k + 1}} - {F_{i,j,k}}} )}^2}} }}, \end{aligned}$$
where $\xi$ is a positive number introduced to avoid the denominator becoming zero, and $\xi$ was set to ${10^{ - 8}}$ in the study.

2.4.3. Image update

In the iterative process for the alternate update of TIE-Hom, the reconstructed images can be updated as follows:

$${F_T} = {F_T} + \sigma {F_{T - 1}},$$
where T is the iteration number in the alternate update of TIE-Hom, $\sigma$ denotes the tuning parameter that can control the ratio between ${F_T}$ and ${F_{T - 1}}$, ${F_T}$ and ${F_{T - 1}}$ are the T-th and (T-1)-th iteration results, respectively.

2.4.4. Min-Max normalization

To ensure the stability and accuracy of reconstruction, the Min-Max normalization is employed, and it can be expressed as follows:

$${F_T} = ({{{({F_{T - 1}})}_{\max }} - {{({F_{T - 1}})}_{\min }}} )\times \frac{{({{F_T} - {{({F_T})}_{\min }}} )}}{{{{({F_T})}_{\max }} - {{({F_T})}_{\min }}}} + {({F_{T - 1}})_{\min }}.$$
where ${({F_T})_{\max }}$ and ${({F_T})_{\min }}$ represent the maximum and minimum values of ${F_T}$. ${({F_{T - 1}})_{\max }}$ and ${({F_{T - 1}})_{\min }}$ denote the maximum and minimum values of ${F_{T - 1}}$, respectively.

2.4.5. Overall workflow

In this study, the proposed CT reconstruction framework based on Eq. (14) is called proposed SART-TV (pSART-TV), and the proposed accelerated CT reconstruction framework based on Eq. (15) is called proposed IFBP-TV (pIFBP-TV). Here, pIFBP-TV is the accelerated version of pSART-TV, and it is the final proposed framework. To explain the pIFBP-TV algorithm more clearly, its corresponding iterative process is given, as shown in Fig. 1. For a multi-material biological sample composed of materials with a large difference in the δ/β ratio (γ value), TIE-Hom is first used for phase retrieval to obtain the δ information in the initial projection data ${I^D}$, where γ is set as γ1 that is suitable for one material, and then IFBP-TV is used for CT reconstruction. Next, forward projection is performed based on the obtained CT data to generate phase-contrast projection ${I_1}^D$, and the projection residual between ${I_1}^D$ and ${I^D}$ is calculated. We continue to use the TIE-Hom for phase retrieval of the projection residual. At this time, to correct and suppress the multi-material artifacts generated in the initial stage, γ is set as γ2 that is suitable for another material. Then, IFBP-TV is used to obtain the corresponding CT data. Finally, it is fused with the previous CT data to update the CT images. In the above procedures the phase retrieval and CT reconstruction are combined into one iterative reconstruction process, which helps to simultaneously achieve the purpose of suppressing multi-material artifacts and and low-dose CT reconstruction. By setting the total iteration number as the stopping criterion, and performing the above procedures alternately, high-quality reconstruction results can be finally obtained.

 figure: Fig. 1.

Fig. 1. Flowchart of the pIFBP-TV algorithm.

Download Full Size | PDF

2.5. Parameter selection

The conventional two-step algorithms, such as FBP with TIE-Hom (referred to as FBP) and SART-TV with TIE-Hom (referred to as SART-TV), are used as the comparison algorithms to evaluate the reconstruction performance of the pIFBP-TV. In addition, pSART-TV is also used as a comparison algorithm to compare the reconstruction efficiency. In the pIFBP-TV algorithm, there are seven main parameters: γ1 for the initial stage and γ2 for the iteration stage, the iteration number T in the alternate updations of TIE-Hom, the relaxation parameter σ, the parameter ${\lambda ^\ast }$ in 3D TV regularization term, the iteration number t and the step size ρ in the IFBP process. The settings of γ1 and γ2 depend on the γ values of strong and weak absorbing materials in multi-material biological samples. In the simulation, the δ value and β value of the materials can be found at http://henke.lbl.gov/optical_constants/getdb2.HTML, and then the γ1 and γ2 values can be calculated. In the real IL-PCCT experiment, the settings of γ1 and γ2 values are obtained from previous experiments. Among them, the weak absorbing tissues of the biological sample is generally set to 500-1000, and the strong absorbing tissues is generally set to 100-400. The other parameters adopt a greedy strategy to determine the optimal values one by one [41]. In the simulation experiment, the peak signal-to-noise ratio (PSNR) is used as an evaluation metric to find the optimal value of each parameter in turn while keeping other parameters unchanged. In the real IL-PCCT experiment, the signal-to-noise ratio (SNR) is used in the same way to determine the optimal parameter set. In the simulation, the parameters are set as follows: (1) pIFBP-TV: T = 2, σ = 2, ${\lambda ^\ast } = 0.3$, t = 20, ρ = 0.1; (2) pSART-TV: T = 2, σ = 2, ${\lambda ^\ast } = 0.2$, t = 25; and (3) SART-TV: ${\lambda ^\ast } = 0.2$, t = 25. In the real IL-PCCT experiment, the parameters are set as follows: (1) pIFBP-TV: T = 2, σ = 0.5, ${\lambda ^\ast } = 0.05$, t = 5, ρ = 0.1; (2) pSART-TV: T = 2, σ = 2, ${\lambda ^\ast } = 0.05$, t = 8; (3) SART-TV: ${\lambda ^\ast } = 0.05$, t = 8.

2.6. Evaluation

The PSNR, the structural similarity index (SSIM), the contrast-to-noise ratio (CNR), the SNR and the modulation transfer function (MTF) are used to evaluate the reconstruction performance of all algorithms in this study. They can be expressed as follows:

  • (1) The PSNR is a traditional measure of image quality, and a larger PSNR value between the reconstructed image and the reference image means better quality. It is defined as follows:
    $$MSE(x,y) = \frac{1}{{M \times M}}\sum\limits_{i = 1}^M {\sum\limits_{j = 1}^M {{{({x_{i,j}} - {y_{i,j}})}^2}} } ,$$
    $$PSNR(x,y) = 10{\log _{10}}(\frac{{Pea{k^2}}}{{MSE}}),$$
    where MSE is the mean square error function, x denotes the reference image, y denotes the reconstructed image, $M \times M$ is the size of x and y, ${x_{i,j}}$ and ${y_{i,j}}$ represent the pixel intensities of x and y in some pixel (i, j), respectively, and Peak represents the largest pixel intensity in the normalized image, i.e., Peak was set to 255 in our study.
  • (2) The SSIM is used to assess the similarity between the reconstructed image and the reference image, and a larger SSIM value means a higher similarity. The SSIM is defined as follows:
    $$SSIM(x,y) = \frac{{(2{u_x}{u_y} + {c_1})(2{\sigma _{xy}} + {c_2})}}{{(u_x^2 + u_y^2 + {c_1})(\sigma _x^2 + \sigma _y^2 + {c_2})}},$$
    where ${u_x}$ and ${u_y}$ are the mean values of x and y, respectively; $\sigma _x^2$ and $\sigma _y^2$ denote the variances of x and y, respectively; ${\sigma _{xy}}$ is the covariance of x and y; ${c_1}$ and ${c_2}$ are two constants used to stabilize the division operation with weak denominator, i.e. ${c_1} = {({k_1}L)^2}$, ${c_2} = {({k_2}L)^2}$; in our study, L was set to 255, ${k_1}$ and ${k_2}$ were set to 0.01 and 0.03, respectively.
  • (3) The CNR is adopted to evaluate the contrast between two different tissues, and a larger CNR value represents a higher contrast [42]. The CNR is defined as follows:
    $$CNR({y_1},{y_2}) = \frac{{{u_{{y_1}}} - {u_{{y_2}}}}}{{\sqrt {0.5(\sigma _{{y_1}}^2 + \sigma _{{y_2}}^2)} }},$$
    where ${y_1}$ and ${y_2}$ represent local regions of two different tissues in the same reconstructed image, respectively.${u_{{y_1}}}$ and ${u_{{y_2}}}$ denote the mean intensity (signal) values of ${y_1}$ and ${y_2}$, respectively; and $\sigma _{{y_1}}^{}$ and $\sigma _{{y_2}}^{}$ are the standard deviations of ${y_1}$ and ${y_2}$, respectively.
  • (4) The SNR is a no-reference quality estimator for images, and a larger SNR value indicates the better image quality [42]. The SNR is defined as follows:
    $$SNR({y_1}) = 20{\log _{10}}(\frac{{{u_{{y_1}}}}}{{{\sigma _{{y_2}}}}}).$$
    where ${y_1}$ and ${y_2}$ represent two local regions of the tissue and background in the same reconstructed image, respectively. ${u_{{y_1}}}$ is the mean value of ${y_1}$, and $\sigma _{{y_2}}^{}$ is the standard deviations of ${y_2}$.
  • (5) The MTF is used to evaluate the resolution of CT images. For the MTF calculation, an edge spread function (ESF) is first obtained by utilizing a sharp edge phantom. A line spread function (LSF) is then obtained by differentiating ESF. Finally, the MTF curves can be obtained from the fast Fourier transform (FFT) of the LSF. Further details about the MTF are provided in literature [43].

3. Simulation

3.1. Simulation data

To evaluate the reconstruction performance of the pIFBP-TV, a 3D microvasculature numerical phantom is designed for the simulation study. Figure 2(a) and Fig. 2(b) show the 110th CT image and surface-rendering image of the phantom, respectively. In the phantom, less dense water and dense water is used to simulate the parts inside and outside the microvasculature, respectively. Polypropylene (PP) is used to simulate the heterogeneous tissues. The distribution areas of the three materials are represented in Fig. 2(a), and they are marked with numbers “1”, “2” and “3”, respectively. The above materials are defined by δ and β values under an energy of 24 keV. The corresponding numerical information can be found at http://henke.lbl.gov/optical_constants/getdb2.HTML. The specific parameters of the phantom are shown in Table 1. The volume of the 3D microvasculature phantom is 512×512×200 voxels, and the size of the voxel pitch is 9 µm×9 µm ×9 µm. In the simulation, the simulated X-ray energy is 24 keV, and the SDD is adjusted to 20 cm. The detector is modeled as a straight-line array with 724 bins, and the size of the reconstructed images is 512×512 pixels. The total number of reconstructed slices is 200. Sixty-four projections are selected from the full projection dataset (512) with equal angle intervals to simulate sparse-view projections. The two-step algorithms, such as FBP and SART-TV, are reconstructed under the conditions γ1 = 3512.2 and γ2 = 2072.5, respectively. pSART-TV and pIFBP-TV are performed with a γ value of 3512.2 in the initial stage and 2072.5 in the iterative stage. All experiments are carried out in MATLAB 2017b on a workstation computer equipped with an Intel Xeon E5-2640 CPU at 2.4 GHz. The 3D surface-rendering images are generated from a sequential series of CT slices using Amira 6.3.0 software (Visage Imaging, Berlin, Germany).

 figure: Fig. 2.

Fig. 2. The 3D microvasculature numerical phantom. (a) The 110th slice of the 3D microvasculature phantom. (b) The surface-rendering image of the 3D microvasculature phantom for 3D visualization. The areas “1” and “2” represent less dense and dense water materials, respectively, and the area “3” represents polypropylene material. The upper right corner of (a) is the reference profile based on the pixel value of the purple line. The scale bar represents 1 mm.

Download Full Size | PDF

Tables Icon

Table 1. The compositions in the 3D microvasculature numerical phantom.

3.2. Simulation results

The surface-rendering images and the 110th slices of the 3D microvasculature phantom in which all slices were reconstructed using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms under 64 projections are shown in Fig. 3. The 110th slice of the 3D microvasculature phantom reconstructed using the pIFBP-TV algorithm under 64 projections without (TIE-Hom) phase retrieval procedure was provided in Supplement 1. When the set value of γ is suitable for PP, i.e., γ = γ1, the reconstruction image using FBP is overall blurred under the 512 projections, and its corresponding 3D image shows the loss of structural information, especially in and around the PP regions, as shown by the purple arrows (Fig. 3(a) and Fig. 3(e)). Under the 64 projections, the reconstruction result of FBP not only has serious blurring but also exhibits severe streak artifacts (Fig. 3(b) and Fig. 3(f)). SART-TV can effectively suppress streak artifacts, but the blurring still remains in the CT images. Its corresponding 3D image can also display the loss and destruction of 3D structures due to the blurring (Fig. 3(c) and Fig. 3(g)). When the set value of γ is suitable for water, i.e., γ = γ2, phase contrast fringes appear at the edge and middle of the PP regions in the reconstructed image using FBP with the 512 projections. It is also evident that the PP regions are damaged by phase contrast fringes in the 3D image, as shown by blue arrows (Fig. 3(i) and Fig. 3(m)). Under the 64 projections, the reconstructed image using FBP is severely affected by streak artifacts, and the false structures formed by the streak artifacts also appear in the corresponding 3D image (Fig. 3(j) and Fig. 3(n)). SART-TV can suppress streak artifacts, but phase contrast fringes around PP regions still exist (Fig. 3(k) and Fig. 3(o)). Figure 3(d) and Fig. 3(l) are the CT images reconstructed by pSART-TV and pIFBP-TV under the 64 projections, respectively. Figure 3(h) and Fig. 3(p) are their corresponding 3D images. Compared with the FBP and SART-TV, pSART-TV and pIFBP-TV can not only effectively remove streak artifacts but also have a superior suppression effect on multi-material artifacts. The corresponding 3D images show that the pSART-TV and pIFBP-TV can suppress the influence of the two types of artifacts on the 3D structures and have better visual effects.

 figure: Fig. 3.

Fig. 3. The surface-rendering images and the 110th slices of the 3D microvasculature phantom where all slices were reconstructed using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms under 64 projections. (a, b, c) are the 110th slices reconstructed using FBP with 512 projections, FBP with 64 projections, and SART-TV with 64 projections, respectively, in which the γ value is set to 3512.2. (i, j, k) are the 110th slices reconstructed by FBP under 512 projections, FBP under 64 projections, and SART-TV under 64 projections, respectively, when the γ value is set to 2072.5. (d) and (l) are the 110th slices reconstructed using pSART-TV and pIFBP-TV under 64 projections, respectively. (e-h) and (m-p) are the surface-rendering images corresponding to the reconstruction results of the above methods. The upper right corners of (a-d) and (i-l) are the profiles located at the pixel position labeled with the purple line in Fig. 3(a). The lower right corners of (a-d) and (i-l) are the enlarged images of the green rectangle marked in (a). The purple arrows and the blue arrows mark the multi-material artifacts generated when the γ value is set to 3512.2 and 2072.5, respectively. The reconstructed images are normalized with a range of [0, 255] and the display window for CT images is [0, 255].

Download Full Size | PDF

To better evaluate the quality of the reconstruction images at the edges of the structures, the intensity profiles based on the same position in the reconstructed images are drawn in the upper right corner of Figs. 3(a-d) and Figs. 3(i-l), which is marked with purple lines in Fig. 2(a). The results show that the intensity profiles obtained by pSART-TV and pIFBP-TV each have a better smoothness and fewer fluctuations in and around the PP regions, which are closer to the reference profiles. In addition, the regions of interests (ROIs) in the reconstructed images, as marked by the green rectangle in Fig. 3(a), are enlarged in the lower right corner of Figs. 3(a-d) and Figs. 3(i-l). It is obvious that compared with the other algorithms, pSART-TV and pIFBP-TV can greatly suppress multi-material artifacts, and the reconstructed structures are clearer.

3.3. Assessments

Table 2 shows the PSNR and SSIM values, as well as the reconstruction time of the images for different methods under 64 projections. Compared with FBP and SART-TV, pIFBP-TV has the highest PSNR value and pSART-TV has the highest SSIM value. pSART-TV and pSART-TV have the similar reconstruction quality. However, pIFBP-TV has the shorter reconstruction time than pSART-TV, indicating that it can improve the reconstruction efficiency compared with pSART-TV.

Tables Icon

Table 2. Quantitative results regarding the reconstructed images for different methods.a

3.4. Validation of effectiveness

To analyze the performance of the conventional two-step algorithms for eliminating multi-material artifacts, the experiments are performed. Figures 4 (a-e) show the 110th CT images reconstructed by FBP with the full projection data, in which the phase retrieval are performed using TIE-Hom with δ/β ratios 2072.5, 2400, 2800, 3100 and 3512.2, respectively. As seen in the corresponding CT images and the enlarged images, regardless of the δ/β ratios, the blue and purple rectangles in the images show that multi-material artifacts still remain in the reconstruction results. The reconstruction results of FBP show the highest PSNR value when the δ/β ratio is 2800 and the highest SSIM value when the δ/β ratio is 2400 (Fig. 4(f) and Fig. 6(g)). The quantitatively assessed values are all lower than those of the quantitative results of pSART-TV and pIFBP-TV, as shown in Table 2. The above experimental results indicate that the conventional two-step algorithms cannot effectively eliminate multi-material artifacts. The qualitative and quantitative results prove that pSART-TV and pIFBP-TV can effectively suppress multi-material artifacts and improve the image quality. Specifically, the visual and quantitative assessments for the pSART-TV and the pIFBP-TV are very close, but the latter has higher reconstruction efficiency compared to the former. To analyze the robustness to noise of the proposed algorithm, the Gaussian noise with different noise levels was introduced into the projection data and the corresponding image reconstructions were provided in Supplement 1.

 figure: Fig. 4.

Fig. 4. The effect of different δ/β on the results of FBP reconstruction. (a-e) are the 110th CT images reconstructed by FBP under full projection data after phase retrieval with δ/β set to 2072.5, 2400, 2800, 3100 and 3512.2, respectively. (h) is the curve of PSNRs corresponding to the above CT images. (p) is the curve of SSIMs corresponding to the above CT images. The blue and purple rectangles in (a) mark the regions of interest (ROIs), and the corresponding enlarged images are located in (f) and (g), respectively. The scale bar represents 1 mm.

Download Full Size | PDF

4. Real data experiment

4.1. Data acquisition

In this study, IL-PCCT reconstruction was performed using an ex vivo rat maxilla sample provided by the Tianjin Medical University School of Stomatology. The experiment was approved by the Research Ethics Committee of the Tianjin Medical University School of Stomatology. The corresponding IL-PCCT data was acquired from the X-ray imaging and biomedical application beamline (BL13W1) of Shanghai Synchrotron Radiation Facility (SSRF) in China. The beamline can be produced by a 16-pole wiggler source on a 3.5 GeV storage ring and monochromated by a double-crystal monochromator with silicon 111 crystals tuning into an energy range from 8 keV to 72.5 keV. The details of the BL13W1 beamline at SSRF were provided online (http://www.sinap.ac.cn/e-ssrf/beamlines/bl13w1). In this study, the photon energy was adjusted to 33 keV, the exposure time was set to 5 ms per projection, and the SDD was adjusted to 0.9 m. A charge-coupled device (CCD) camera (VHR1:1; Photonic Science Ltd., Robertsbridge, UK) with a 36 mm (horizontal) 5 mm (vertical) field of view (FOV) was adopted as the imaging detector, with an effective pixel size of 9 µm × 9 µm, and the size of the projection image is 3992 ×513 pixels. In addition, 10 dark images (dark signal in the absence of photons) and 10 flat images (with no sample in the beam) were used to perform a dark-field correction and a flat-field correction, respectively [44]. A total of 820 projections at equidistant angles from 0° to 180° were collected, and 205 projections were selected at equal angle intervals to simulate the sparse-view CT reconstruction experiment. In this study, the TIE-Hom method was used for phase retrieval. According to previous experiments and experience, the appropriate γ value of soft tissue (gingiva, etc.) and hard tissue (tooth, etc.) in the maxilla sample was 1000 and 100, respectively. Thus, FBP with TIE-Hom and SART-TV with TIE-Hom were reconstructed under the conditions where γ1 = 1000 and γ2 = 100. The γ values of pSART-TV and pIFBP-TV were set as 1000 in the initial stage and 100 in the iterative stage for reconstruction. The volume of the reconstructed CT data (with redundant data having be cropped to fit the extent of the sample) was 1991×1991×482 voxels. In this experiment, only the 135th reconstructed image was displayed for comparison and analysis. 470 reconstructed images containing the maxilla sample information were selected to generate the 3D volume-rendering images.

4.2. Experimental results

The volume-rendering images and the 135th slices of the rat maxilla sample in which all slices were reconstructed using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms under 205 projections are shown in Fig. 5. The 135th slice of the rat maxilla sample reconstructed using the pIFBP-TV algorithm under 205 projections without (TIE-Hom) phase retrieval procedure was provided in Supplement 1. When the γ value is set to 1000, the 135th slices reconstructed by FBP with 820 projections, FBP with 205 projections and SART-TV with 205 projections are respectively shown in Figs. 5(a-c), and the corresponding 3D images of those are shown in Figs 5 (e-g), respectively. When the γ value is set to 100, the 135th slices reconstructed by FBP with 820 projections, FBP with 205 projections and SART-TV with 205 projections are respectively shown in Figs. 5(i-k), and the corresponding 3D images of those are shown in Figs. 5(m-o), respectively. The 135th slice reconstructed by pSART-TV and pIFBP-TV with 205 projections are respectively displayed in Fig. 5(d) and Fig. 5(l), and the corresponding 3D images are shown in Fig. 5(h) and Fig. 5(p), respectively. The lower left corners located at Figs. 5(a-d) show the enlarged images of the ROIs (as marked with the blue rectangle in Fig. 5(a)) in reconstructed images, respectively. In addition, the enlarged images of the ROIs (as marked with yellow rectangle in Fig. 5(i)) in reconstructed images are shown in the lower left corners of Figs. 5(i-l), respectively. Correspondingly, the intensity profiles drawn with the same position (as labeled with the purple dotted line in blue rectangular ROIs) in reconstructed images are presented in the lower right corner of Figs. 5(a-d), and the intensity profiles drawn with the same positions (as labeled with the purple dotted line in the yellow rectangular ROIs) in reconstructed images are presented in the lower right corner of Figs. 5(m-o). The lower right corners of Figs. 5(e-h) and Figs. 5(m-p) present the enlarged images of the ROIs for the corresponding 3D images, respectively.

 figure: Fig. 5.

Fig. 5. The volume-rendering images and the 135th slices of the rat maxilla sample where all slices were reconstructed using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms under 205 projections. (a, b, c) are the 135th slices reconstructed by FBP under 820 projections, FBP under 205 projections, and SART-TV under 205 projections, respectively, when the γ value is set to 1000. (i, j, k) are the 135th slices reconstructed by FBP under 820 projections, FBP under 205 projections, and SART-TV under 205 projections, respectively, when the γ value is set to 100. (d) and (l) are the 135th slices reconstructed using pSART-TV and pIFBP-TV under 205 projections, respectively. (e-h) and (m-p) are the 3D volume-rendering images for the overall slices reconstructed by the above algorithms, and they correspond to (a-d) and (i-l), respectively. The lower left corners of (a-d) and (i-l) show the enlarged images of the blue rectangle marked in (a) and the yellow rectangle marked in (i), respectively. The purple dotted line in the blue and yellow rectangles are used to draw the intensity profiles of the reconstructed images. The lower right corners of (a-d) and (i-l) exhibit the intensity profiles located at the pixel position labeled with the purple dotted lines in the blue rectangle and yellow rectangle. The lower right corners of (e-h) and (m-p) display the enlarged images of the green rectangle marked in (e). The red rectangle area (40 × 40 pixels), the yellow blue rectangle area (40 × 40 pixels) and the blue rectangle area (40 × 40 pixels) in Fig. 5(a) are the background region, the soft tissue region and the dentin region, respectively. These regions were used to calculate CNR and SNR values. The blue and yellow arrows indicate edge details for visual comparison. The reconstructed images were normalized and then their intensity range is [0, 255]. The display window for CT images is [0, 255]. The scale bar represents 1 mm.

Download Full Size | PDF

Under the 205 projections, when γ was set to 1000, the reconstructed results of FBP contained severe streak artifacts and blurring (Fig. 5(b) and Fig. 5(f)). SART-TV can eliminate only streak artifacts, but blurring still exist, and many details are not clear. Many pore structures information were lost in the corresponding 3D image, resulting in over-smooth effects on the surface of an alveolar bone (Fig. 5(c) and Fig. 5(g)). When γ was set to 100, FBP reconstruction results had poor visual effects (Fig. 5(j) and Fig. 5(n)). SART-TV can effectively suppress streak artifacts, but the serious phase contrast fringes still remain. The edges of soft tissues (e.g., gums) and hard tissues (e.g., alveolar bone), were darkened, which would have a substantial impact on subsequent image analysis. The 3D image also showed that the tissues information were affected by phase contrast fringes and were not continuous (Fig. 5(k) and Fig. 5(o)). Compared with the two-step algorithms, pSART-TV and pIFBP-TV could effectively remove not only the streak artifacts introduced by sparse-view projection data but also the multi-material artifacts introduced by the phase retrieval stage (Fig. 5(d) and Fig. 5(l)). The pSART-TV and pIFBP-TV could effectively eliminate low-frequency artifacts and clearly distinguish the structure between enamel and dentin, as shown by the blue arrow in the blue rectangular ROI in Fig. 5(d). Additionally, they could also suppress the phase contrast fringes, as indicated by the yellow arrow in the yellow rectangular ROI in Fig. 5(i). The 3D images showed that pSART-TV and pIFBP-TV achieved ideal visual effects, effectively preserving the structures and details. The tissues information in the images were also continuous and smooth (Fig. 5(h) and Fig. 5(p)).

The same situation can also be observed by the intensity profiles in Figs. 5(a-d) and Figs. 5(m-p). When γ was set to 1000, due to the existence of blurring, the reconstructed results of FBP and SART-TV did not have any fluctuations in the intensity profiles between enamel and dentin, which meant that the enamel and dentin could not be distinguished by FBP and SART-TV. From (Fig. 5(d) and Fig. 5(l)), it can be seen that the pSART-TV and pIFBP-TV could effectively differentiate these different tissues, which suggested that pSART-TV and pIFBP-TV enabled to eliminate the blurring. When γ was set to 100, the intensity profiles of the FBP and SART-TV both had a sharp jump near the boundaries of different tissue regions because of the phase contrast fringes (Fig. 5(i) and Fig. 5(j)). Seen from (Fig. 5(k) and Fig. 5(l)), it could be found that the intensity profiles of pSART-TV and pIFBP-TV did not have the aforementioned jump and they were much smoother than those of FBP and SART-TV, which suggested that pSART-TV and pIFBP-TV enabled to eliminate the phase contrast fringes.

4.3. Results analysis

In this study, MTF curves were employed to evaluate the resolution of reconstructed images, and they were drawn based on the same position in the 135th slices of the rat maxilla sample reconstructed using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms, which were marked with the red line in CT image in Fig. 6 (a). The MTF curves were shown in Fig. 6 (a), and the enlarged image of the red rectangle marked in Fig. 6 (a) was shown in Fig. 6 (b). Additionally, 10 CT slices were randomly selected from the reconstructed slices of the sample to evaluate the image contrast and quality, in which the background region (as marked with red rectangle in Fig. 5(a)), the soft tissue region (as marked with yellow rectangle in Fig. 5(a)) and the dentin region (as marked with the blue rectangle in Fig. 5(a)) were utilized to calculate CNR and SNR values. The CNR values of the dentin regions and the soft tissue regions were calculated based on the Eq. (22), and the results were shown in Fig. 6 (c). The SNR values of the dentin regions and the background regions were calculated based on the Eq. (23), and the results were shown in Fig. 6 (d). As seen in Figs. 6(a-c), when the γ value was close to the soft tissue in the maxilla sample, i.e., γ = γ1, the images reconstructed by the conventional two-step algorithms displayed a high contrast between the hard regions and the soft tissue regions, while there was severe blurring remained in them (i.e., low resolution). When the γ value was close to the hard tissue in the maxilla sample, i.e., γ = γ2, the images reconstructed by conventional two-step algorithms showed clear details, in which the blurring was greatly reduced. However, they had very low contrast between the hard regions and the soft tissue regions (i.e., low image contrast). Compared with the above algorithms, pSART-TV and pIFBP-TV had better MTF and CNR results, which demonstrated that they could better preserve image contrast and resolution simultaneously. The SNR results in Fig. 6(d) showed that pSART-TV and pIFBP-TV could also effectively reduce artifacts and maintain high image quality.

 figure: Fig. 6.

Fig. 6. Quantitative results of the reconstructed images for the rat maxilla sample. (a) The MTF curves of the reconstructed images using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms. The red line marked in the CT image in the upper left corner of (a) was used to calculate MTF. (b) The enlarged image of the red rectangle marked in (a). (c) The CNR values of the images reconstructed by the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms. (d) The SNR values of the images reconstructed by the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms.

Download Full Size | PDF

Table 3 showed the overall reconstruction time for the rat maxilla sample using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms. It could be clearly found that the reconstruction time of pIFBP-TV was greatly reduced compared with pSART-TV, indicating that pIFBP-TV had a lower computational cost and faster reconstruction speed, namely the higher reconstruction efficiency.

Tables Icon

Table 3. Reconstruction time (s) of the different algorithmsa

5. Conclusion and discussion

Multi-material and streak artifacts are two important problems for IL-PCCT. To alleviate these problems simultaneously, we proposed an accelerated 3D iterative reconstruction algorithm, namely pIFBP-TV. In pIFBP-TV, TIE-Hom is updated alternately for different areas of the multi-material sample to obtain the corresponding phase shifts, and then, IFBP algorithm combined with 3D TV regularization is used to reconstruct CT images (i.e., 3D refractive index decrement δ distribution of sample). Since δ is assumed to be proportional to β in this study, the 3D absorption index β distribution of sample can be found by the reconstructed CT images. Finally, the forward projection is performed based on the reconstructed CT images and the obtained 3D absorption index β distribution of sample to generate phase-contrast projections, and the data fidelity term based on the Euclidean distance between the generated phase-contrast projections and original phase-contrast projections is used to improve reconstruction accuracy. It is worth noting that the phase retrieval and CT reconstruction are merged into one step in the proposed algorithm, which can help to jointly suppress multi-material artifacts and streak artifacts. To evaluate the effectiveness and efficiency of the proposed algorithm, the 3D microvasculature phantom and the IL-PCCT data of rat maxilla samples were utilized. In the experiments, the widely used two-step algorithms (e.g., FBP and SART-TV) were used for comparison, and the proposed one-step algorithm pSART-TV was also used for comparison. Simulation and experimental results demonstrated that pIFBP-TV could suppress multi-material and streak artifacts effectively, and reach a balance between image contrast and resolution. Specifically, it had a lower computational cost and faster reconstruction speed compared with pSART-TV, which is conducive to practical application. Overall, the proposed algorithm can serve as an effective tool to solve both phase retrieval problem and low-dose CT reconstruction problem of multi-material objects simultaneously, and thus enabling to promote the wide application and development of PB-PCCT in biomedical and pre-clinical studies.

The innovations of this study are summarized as follows: First, this paper for the first time presents a multi-task 3D iterative reconstruction framework that can suppress multi-material artifacts and streak artifacts simultaneously. Second, the proposed algorithm is a one-step approach that combines the phase retrieval and CT reconstruction into one step, and this one-step approach can help to further improve the reconstruction performance of IL-PCCT compared with the conventional two-step approaches. Third, to the best of our knowledge, it is the first that the TIE-Hom is used to correct the multi-material artifacts by updating alternately for different areas of the multi-material sample.

Herein, we would like to point out the limitations of our study. First, the selections of γ values may be difficult for multi-material samples with unknown compositions. From our point of view, a feasible strategy is to determine the multiple-parameter values one by one which is commonly used in CT reconstruction [41], and it may be helpful for this problem. The strategy first selects a set of γ values at different scales, and then performs visual evaluation and quantitative evaluation (e.g., CNR and SNR) on the reconstructed CT images in order to find the optimal γ values that are suitable for different areas of multi-material samples. Second, the reconstruction time of the proposed algorithm is still relatively long, and it will significantly increase with the size of the sample (i.e., raw projection data size), this is a great problem for the practical application of the proposed algorithm. In our opinions, a feasible strategy is to reconstruct the ROI projection data of the sample (i.e., interior reconstruction [45,46,47]) using the proposed iterative reconstruction framework, which can help to further reduce the reconstruction time and the radiation dose. In future work, the adaptivity of the parameters in the pIFBP-TV algorithm will be researched deeply, and further research will be conducted using the graphics processing unit (GPU)-based speedup technique to greatly reduce the reconstruction time to a reasonable level for practical implementation. Moreover, we will try to incorporate more geometric priors into regularization terms to further improve the performance of the proposed algorithm. In addition, future studies will also be performed by plugging regularization terms or physical model into deep convolutional neural network to optimize the combined one-step reconstruction framework, in which the generalization ability of the proposed framework can also be enhanced.

Funding

National Natural Science Foundation of China (82071922, 82102037, 81671683); Natural Science Foundation of Tianjin City (16JCYBJC28600); Tianjin Municipal Education Commission (2020KJ208).

Acknowledgements

The authors would like to thank the staffs from beamline (BL13W1) of SSRF, China, for their kindly assistance in our experiments. Additionally, the authors would like to thank the anonymous reviewers for their extremely useful suggestions for improving the quality of the paper.

Disclosures

The authors declare that there are no conflicts of interest related to this article.

Data availability

The data and code that support the findings of this study are available from the corresponding author upon reasonable request.

Supplemental document

See Supplement 1 for supporting content.

References

1. X.Z. Wu, H. Liu, and A.M. Yan, “Optimization of X-ray phase-contrast imaging based on in-line holography,” Nucl. Instrum. Methods Phys. Res., Sect. B 234(4), 563–572 (2005). [CrossRef]  

2. A. Momose, T. Takeda, Y. Itai, and K. Hirano, “Phase–contrast X–ray computed tomography for observing biological soft tissues,” Nat. Med. 2(4), 473–475 (1996). [CrossRef]  

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

4. A. Momose, “Phase-sensitive imaging and phase tomography using X-ray interferometers,” Opt. Express 11(19), 2303–2314 (2003). [CrossRef]  

5. A. Bravin, V. Mocella, P. Coan, A. Astolfo, and C. Ferrero, “A numerical wave-optical approach for the simulation of analyzer-based X-ray imaging,” Opt. Express 15(9), 5641–5648 (2007). [CrossRef]  

6. A. Snigirev, I. Snigireva, V. Kohn, S. Kuznetsov, and I. Schelokov, “On the possibilities of X-ray phase contrast microimaging by coherent high-energy synchrotron radiation,” Rev. Sci. Instrum. 66(12), 5486–5492 (1995). [CrossRef]  

7. 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), 6296–6304 (2005). [CrossRef]  

8. C. K. Hagen, P. Maghsoudlou, G. Totonelli, P. C. Diemoz, M. Endrizzi, L. Rigon, R.-H. Menk, F. Arfelli, D. Dreossi, E. Brun, P. Coan, A. Bravin, P. De Coppi, and A. Olivo, “High contrast microstructural visualization of natural acellular matrices by means of phase-based X-ray tomography,” Sci. Rep. 5(1), 18156 (2016). [CrossRef]  

9. R.J. Xuan, X.Y. Zhao, D.D. Hu, J.J. Jian, T.L. Wang, and C.H. Hu, “Three-dimensional visualization of the microvasculature of bile duct ligation-induced liver fibrosis in rats by x-ray phase-contrast imaging computed tomography,” Sci. Rep. 5(1), 11500 (2015). [CrossRef]  

10. A. Bravin, P. Coan, and P. Suortti, “X-ray phase-contrast imaging: from pre-clinical applications towards clinics,” Phys. Med. Biol. 58(1), R1–R35 (2013). [CrossRef]  

11. W. J. Lv, X.Y. Zhao, D.D. Hu, X.H. Xin, L.L. Qin, and C.H. Hu, “Insight into Bile Duct Reaction to Obstruction from a Three-dimensional Perspective Using ex Vivo Phase-Contrast CT,” Radiology 299(3), 597–610 (2021). [CrossRef]  

12. B. Patrycja, M. Sheridan, M Mikkaela, P. Serena, T. Giuliana, D. Christian, Z. Fabrizio, A. Fulvia, D. Diego, F. Jane, P. Zdenka, C. Marian, Q. Harry, D. Matthew, N. Yakov, T. Darren, B. Patrick, and T. Gureyev, “High-Resolution X-Ray Phase-Contrast 3D Imaging of Breast Tissue Specimens as a Possible Adjunct to Histopathology,” IEEE Trans. Med. Imaging 37(12), 2642–2650 (2018). [CrossRef]  

13. M Töpperwien, F. van der Meer, C. Stadelmann, and T. Salditt, “Three-dimensional virtual histology of human cerebellum by X-ray phase-contrast tomography,” Proc. Natl. Acad. Sci. U. S. A. 115(27), 6940–6945 (2018). [CrossRef]  

14. R. C. Chen, L. Rigon, and R. Longo, “Quantitative 3D refractive index decrement reconstruction using single-distance phase-contrast tomography data,” J. Phys. D: Appl. Phys. 44(49), 495401 (2011). [CrossRef]  

15. M. Langer, P. Cloetens, A. Pacureanu, and F. Peyrin, “X-ray in-line phase tomography of multimaterial objects,” Opt. Lett. 37(11), 2151–2153 (2012). [CrossRef]  

16. H. T. Li, F. Schaff, L. Croton, K. S. Morgan, and M. J. Kitchen, “Quantitative material decomposition using linear iterative near-field phase retrieval dual-energy x-ray imaging,” Phys. Med. Biol. 65(18), 185014 (2020). [CrossRef]  

17. S. Mohammadi, E. Larsson, F. Alves, S. Dal Monego, S. Biffi, C. Garrovo, A. Lorenzon, G. Tromba, and C. Dullin, “Quantitative evaluation of a single-distance phase-retrieval method applied on in-line phase-contrast images of a mouse lung,” J. Synchrotron Radiat. 21(4), 784–789 (2014). [CrossRef]  

18. A. Burvall, U. Lundstrom, P. A. C. Takman, D. H. Larsson, and H. M. Hertz, “Phase retrieval in X-ray phase-contrast imaging suitable for tomography,” Opt. Express 19(11), 10359–10376 (2011). [CrossRef]  

19. T. E. Gureyev, T. J. Davis, A. Pogany, S. C. Mayo, and S. W. Wilkins, “Optical phase retrieval by use of first Born- and Rytov-type approximations,” Appl. Opt. 43(12), 2418–2430 (2004). [CrossRef]  

20. D. Paganin, S. C. Mayo, T. E. Gureyev, P. R. Miller, and S. W. Wilkins, “Simultaneous phase and amplitude extraction from a single defocused image of a homogeneous object,” J. Microsc. 206(1), 33–40 (2002). [CrossRef]  

21. M. Beltran, D. Paganin, K. Uesugi, and M. Kitchen, “2D and 3D X-ray phase retrieval of multi-material objects using a single defocus distance,” Opt. Express 18(7), 6423–6436 (2010). [CrossRef]  

22. M. Beltran, D. Paganin, K. Siu, A. Fouras, S. Hooper, D. Reser, and M. Kitchen, “Interface-specific X-ray phase retrieval tomography of complex biological organs,” Phys. Med. Biol. 56(23), 7353–7369 (2011). [CrossRef]  

23. M. Ullherr and S. Zabler, “Correcting multi material artifacts from single material phase retrieved holo-tomograms with a simple 3D Fourier method,” Opt. Express 23(25), 32718–32727 (2015). [CrossRef]  

24. M. Beister, D. Kolditz, and W. A. Kalender, “Iterative reconstruction methods in X-ray CT,” Phys. Medica 28(2), 94–108 (2012). [CrossRef]  

25. D.J. Ji, C.H. Hu, and H. Yang, “Image reconstruction algorithm for in-line phase contrast imaging computed tomography with an improved anisotropic diffusion model,” J. X-Ray Sci. Technol. 23(3), 311–320 (2015). [CrossRef]  

26. S.A. Melli, K.A. Wahid, P. Babyn, J. Montgomery, E. Snead, L.E. Gayed, M. Pettitt, B. Wolkowski, and M. Wesolowski, “A compressed sensing-based reconstruction algorithm for synchrotron source propagation-based X-ray phase contrast computed tomography,” Nucl. Instrum. Methods Phys. Res., Sect. A 806, 307–317 (2016). [CrossRef]  

27. Y.Q. Zhao, M.Y. Sun, D.J. Ji, C.H. Cong, W. J. Lv, Q. Zhao, L.L. Qin, J.J. Jian, X.D. Chen, and C.H. Hu, “An iterative image reconstruction algorithm combined with forward and backward diffusion filtering for in-line X-ray phase-contrast computed tomography,” J. Synchrotron Radiat. 25(5), 1450–1459 (2018). [CrossRef]  

28. Y.Q. Zhao, D.J. Ji, Y.M. Li, X.Y. Zhao, W. J. Lv, X.H. Xin, S. Han, and C.H. Hu, “Three-dimensional visualization of microvasculature from few-projection data using a novel CT reconstruction algorithm for propagation-based X-ray phase-contrast imaging,” Biomed. Opt. Express 11(1), 364–387 (2020). [CrossRef]  

29. M.T. Zheng, Y.Q. Zhao, S. Han, D.J. Ji, Y.M. Li, W.J. Lv, X.H. Xin, X.Y. Zhao, and C. Hu, “Iterative reconstruction algorithm based on discriminant adaptive-weighted TV regularization for fibrous biological tissues using in-line X-ray phase-contrast imaging,” Biomed. Opt. Express 12(4), 2460–2483 (2021). [CrossRef]  

30. S. Han, Y.Q. Zhao, F.Z. Li, D.J. Ji, Y.M. Li, M.T. Zheng, W.J. Lv, X.H. Xin, B.N. Qi, and C.H. Hu, “A dual-path deep learning reconstruction framework for propagation-based X-ray phase-contrast computed tomography with sparse-view projections,” Opt. Lett. 46(15), 3552–3555 (2021). [CrossRef]  

31. A. Kostenko, K. J. Batenburg, A. King, S. E. Offerman, and L. J. van Vliet, “Total variation minimization approach in in-line X-ray phase-contrast tomography,” Opt. Express 21(10), 12185–12196 (2013). [CrossRef]  

32. L. Tian, J. C. Petruccelli, Q. Miao, H. Kudrolli, V. Nagarkar, and G. Barbastathis, “Compressive X-ray phase tomography based on the transport of intensity equation,” Opt. Lett. 38(17), 3418–3421 (2013). [CrossRef]  

33. A. Ruhlandt, M. Krenkel, M. Bartels, and T. Salditt, “Three-dimensional phase retrieval in propagation-based phase-contrast imaging,” Phys. Rev. A 89(3), 033847 (2014). [CrossRef]  

34. R. Gordon, R. Bender, and G. T. Herman, “Algebraic reconstruction techniques (ART) for three-dimensional electron microscopy and X-ray photography,” J. Theor. Biol. 29(3), 471–481 (1970). [CrossRef]  

35. L. Hehn, R. Gradl, M. Dierolf, K.S. Morgan, D.M. Paganin, and F. Pfeiffer, “Model-Based Iterative Reconstruction for Propagation-Based Phase-Contrast X-Ray CT including Models for the Source and the Detector,” IEEE Trans. Med. Imaging 39(6), 1975–1987 (2020). [CrossRef]  

36. 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 (2010). [CrossRef]  

37. X. Deng, Y. Zhao, and H. Li, “Projection data smoothing through noise-level weighted total variation regularization for low-dose computed tomography,” J. X-Ray Sci. Technol. 27(3), 537–557 (2019). [CrossRef]  

38. O. N. Strand, “Theory and methods related to the singular-function expansion and Landweber’s iteration for integral equations of the first kind,” SIAM J. Numer. Anal. 11(4), 798–825 (1974). [CrossRef]  

39. Z. Tian, X. Jia, K. Yuan, T. Pan, and S. B. Jiang, “Low-dose CT reconstruction via edge-preserving total variation regularization,” Phys. Med. Biol. 56(18), 5949–5967 (2011). [CrossRef]  

40. E. J. Candès, J. Romberg, and T. Tao, “Robust uncertainty principles: exact signal reconstruction from highly incomplete frequency information,” IEEE Trans. Inf. Theory 52(2), 489–509 (2006). [CrossRef]  

41. M. Lohvithee, A. Biguri, and M. Soleimani, “Parameter selection in limited data cone-beam CT reconstruction using edge-preserving total variation algorithms,” Phys. Med. Biol. 62(24), 9295–9321 (2017). [CrossRef]  

42. Z. Berit, J. Mead, T. Deck, R. Tiina, G. Clough, R. Boardman, and S. Philipp, “Soft tissue 3D imaging in the lab through optimised propagation-based phase contrast computed tomography,” Opt. Express 25(26), 33451–33468 (2017). [CrossRef]  

43. M. U. Ghani, B. Gregory, F. Omoumi, B. Zheng, A. Yan, X. Wu, and H. Liu, “Impact of a single distance phase retrieval algorithm on spatial resolution in X-ray inline phase sensitive imaging,” Biomed. Spectrosc. Imaging 8(1-2), 29–40 (2019). [CrossRef]  

44. R. Chen, D. Dreossi, L. Mancini, R. Menk, L. Rigon, T. Xiao, and R. Longo, “PITRE: software for phase sensitive X-ray image processing and tomography reconstruction,” J. Synchrotron Radiat. 19(5), 836–845 (2012). [CrossRef]  

45. Y. Ye, H. Yu, and G. Wang, “Exact Interior Reconstruction with Cone-Beam CT,” Int. J. Biomed. Imaging 2007(3), 10693 (2007). [CrossRef]  

46. E. A. Rashed and H. Kudo, “Towards high-resolution synchrotron radiation imaging with statistical iterative reconstruction,” J. Synchrotron Radiat. 20(1), 116 (2013). [CrossRef]  

47. L. Felsner, S. Kaeppler, A. Maier, and C. Riess, “Truncation correction for x-ray phase-contrast region-of-interest tomography,” IEEE Trans. Comput. Imaging 6, 625–639 (2020). [CrossRef]  

Supplementary Material (1)

NameDescription
Supplement 1       Supplement 1

Data availability

The data and code that support the findings of this study are available from the corresponding author upon reasonable request.

Cited By

Optica participates in Crossref's Cited-By Linking service. Citing articles from Optica Publishing Group journals and other participating publishers are listed here.

Alert me when this article is cited.


Figures (6)

Fig. 1.
Fig. 1. Flowchart of the pIFBP-TV algorithm.
Fig. 2.
Fig. 2. The 3D microvasculature numerical phantom. (a) The 110th slice of the 3D microvasculature phantom. (b) The surface-rendering image of the 3D microvasculature phantom for 3D visualization. The areas “1” and “2” represent less dense and dense water materials, respectively, and the area “3” represents polypropylene material. The upper right corner of (a) is the reference profile based on the pixel value of the purple line. The scale bar represents 1 mm.
Fig. 3.
Fig. 3. The surface-rendering images and the 110th slices of the 3D microvasculature phantom where all slices were reconstructed using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms under 64 projections. (a, b, c) are the 110th slices reconstructed using FBP with 512 projections, FBP with 64 projections, and SART-TV with 64 projections, respectively, in which the γ value is set to 3512.2. (i, j, k) are the 110th slices reconstructed by FBP under 512 projections, FBP under 64 projections, and SART-TV under 64 projections, respectively, when the γ value is set to 2072.5. (d) and (l) are the 110th slices reconstructed using pSART-TV and pIFBP-TV under 64 projections, respectively. (e-h) and (m-p) are the surface-rendering images corresponding to the reconstruction results of the above methods. The upper right corners of (a-d) and (i-l) are the profiles located at the pixel position labeled with the purple line in Fig. 3(a). The lower right corners of (a-d) and (i-l) are the enlarged images of the green rectangle marked in (a). The purple arrows and the blue arrows mark the multi-material artifacts generated when the γ value is set to 3512.2 and 2072.5, respectively. The reconstructed images are normalized with a range of [0, 255] and the display window for CT images is [0, 255].
Fig. 4.
Fig. 4. The effect of different δ/β on the results of FBP reconstruction. (a-e) are the 110th CT images reconstructed by FBP under full projection data after phase retrieval with δ/β set to 2072.5, 2400, 2800, 3100 and 3512.2, respectively. (h) is the curve of PSNRs corresponding to the above CT images. (p) is the curve of SSIMs corresponding to the above CT images. The blue and purple rectangles in (a) mark the regions of interest (ROIs), and the corresponding enlarged images are located in (f) and (g), respectively. The scale bar represents 1 mm.
Fig. 5.
Fig. 5. The volume-rendering images and the 135th slices of the rat maxilla sample where all slices were reconstructed using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms under 205 projections. (a, b, c) are the 135th slices reconstructed by FBP under 820 projections, FBP under 205 projections, and SART-TV under 205 projections, respectively, when the γ value is set to 1000. (i, j, k) are the 135th slices reconstructed by FBP under 820 projections, FBP under 205 projections, and SART-TV under 205 projections, respectively, when the γ value is set to 100. (d) and (l) are the 135th slices reconstructed using pSART-TV and pIFBP-TV under 205 projections, respectively. (e-h) and (m-p) are the 3D volume-rendering images for the overall slices reconstructed by the above algorithms, and they correspond to (a-d) and (i-l), respectively. The lower left corners of (a-d) and (i-l) show the enlarged images of the blue rectangle marked in (a) and the yellow rectangle marked in (i), respectively. The purple dotted line in the blue and yellow rectangles are used to draw the intensity profiles of the reconstructed images. The lower right corners of (a-d) and (i-l) exhibit the intensity profiles located at the pixel position labeled with the purple dotted lines in the blue rectangle and yellow rectangle. The lower right corners of (e-h) and (m-p) display the enlarged images of the green rectangle marked in (e). The red rectangle area (40 × 40 pixels), the yellow blue rectangle area (40 × 40 pixels) and the blue rectangle area (40 × 40 pixels) in Fig. 5(a) are the background region, the soft tissue region and the dentin region, respectively. These regions were used to calculate CNR and SNR values. The blue and yellow arrows indicate edge details for visual comparison. The reconstructed images were normalized and then their intensity range is [0, 255]. The display window for CT images is [0, 255]. The scale bar represents 1 mm.
Fig. 6.
Fig. 6. Quantitative results of the reconstructed images for the rat maxilla sample. (a) The MTF curves of the reconstructed images using the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms. The red line marked in the CT image in the upper left corner of (a) was used to calculate MTF. (b) The enlarged image of the red rectangle marked in (a). (c) The CNR values of the images reconstructed by the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms. (d) The SNR values of the images reconstructed by the FBP, SART-TV, pSART-TV and pIFBP-TV algorithms.

Tables (3)

Tables Icon

Table 1. The compositions in the 3D microvasculature numerical phantom.

Tables Icon

Table 2. Quantitative results regarding the reconstructed images for different methods.a

Tables Icon

Table 3. Reconstruction time (s) of the different algorithmsa

Equations (24)

Equations on this page are rendered with MathJax. Learn more.

T θ ( x , y ) = exp [ B θ ( x , y )  +  i φ θ ( x , y ) ] ,
B θ ( x , y ) = k β ( r , z ) d z ,
φ θ ( x , y ) = k δ ( r , z ) d z ,
I θ D ( r ) = | T θ ( x , y ) P θ D ( r ) | 2 ,
P θ D ( r ) = e i 2 π D / λ i λ D e i π λ D | r | 2 .
( I θ D ( r ) φ θ ( r ) ) = k z I θ D ( r ) ,
φ θ ( r ) = 1 2 γ ln [ F 1 ( F ( I θ D ( r ) ) π γ λ D ( k 2 ) + 1 ) ] .
F = k 1 0 π φ θ v d θ ,
A F = φ ,
F = arg min F 1 2 ( A F φ ) 2 2 + α R ( F ) ,
F = arg min F 1 2 θ  =  0 π H θ D ( P θ D ( A F ) I θ D ) 2 2 + α R ( F ) ,
R 3 D T V ( F ) = i , j , k ( 1 F ) 2 + ( 2 F ) 2 + ( 3 F ) 2 ,
F = arg min F 1 2 θ  =  0 π H θ D ( P θ D ( A F ) I θ D ) 2 2 + α R 3 D T V ( F ) .
F j t = F j t 1 + λ t i = 1 N a i , j i = 1 N H D ( I i D ) n = 1 M a i , n F n t 1 n = 1 M a i , n a i , j , j = 1 , 2 , 3 , , M ,
F t = F t 1 ρ R ( θ  =  0 π H θ D ( P θ D ( A F t 1 ) I θ D ) ) ,
F 3 D T V t + 1 = F I F B P t λ R 3 D T V ( F ) F ,
R 3 D T V ( F ) F i , j , k ( F i , j , k F i 1 , j , k ) + ( F i , j , k F i , j 1 , k ) + ( F i , j , k F i , j , k 1 ) ξ + ( F i , j , k F i 1 , j , k ) 2 + ( F i , j , k F i , j 1 , k ) 2 + ( F i , j , k F i , j , k 1 ) 2   ( F i + 1 , j , k F i , j , k ) ξ + ( F i + 1 , j , k F i , j , k ) 2 + ( F i + 1 , j , k F i + 1 , j 1 , k ) 2 + ( F i + 1 , j , k F i + 1 , j , k 1 ) 2   ( F i , j + 1 , k F i , j , k ) ξ + ( F i , j + 1 , k F i 1 , j + 1 , k ) 2 + ( F i , j + 1 , k F i , j , k ) 2 + ( F i , j + 1 , k F i , j + 1 , k 1 ) 2   ( F i , j , k + 1 F i , j , k ) ξ + ( F i , j , k + 1 F i 1 , j , k + 1 ) 2 + ( F i , j , k + 1 F i , j 1 , k + 1 ) 2 + ( F i , j , k + 1 F i , j , k ) 2 ,
F T = F T + σ F T 1 ,
F T = ( ( F T 1 ) max ( F T 1 ) min ) × ( F T ( F T ) min ) ( F T ) max ( F T ) min + ( F T 1 ) min .
M S E ( x , y ) = 1 M × M i = 1 M j = 1 M ( x i , j y i , j ) 2 ,
P S N R ( x , y ) = 10 log 10 ( P e a k 2 M S E ) ,
S S I M ( x , y ) = ( 2 u x u y + c 1 ) ( 2 σ x y + c 2 ) ( u x 2 + u y 2 + c 1 ) ( σ x 2 + σ y 2 + c 2 ) ,
C N R ( y 1 , y 2 ) = u y 1 u y 2 0.5 ( σ y 1 2 + σ y 2 2 ) ,
S N R ( y 1 ) = 20 log 10 ( u y 1 σ y 2 ) .
Select as filters


Select Topics Cancel
© Copyright 2024 | Optica Publishing Group. All rights reserved, including rights for text and data mining and training of artificial technologies or similar technologies.