A new in-line X-ray phase-contrast computed tomography reconstruction algorithm based on adaptive-weighted anisotropic TpV regularization for insufficient data

In-line X-ray phase-contrast computed tomography (IL-PCCT) is a valuable tool for revealing the internal detailed structures in weakly absorbing objects (e.g. biological soft tissues), and has a great potential to become clinically applicable. However, the long scanning time for IL-PCCT will result in a high radiation dose to biological samples, and thus impede the wider use of IL-PCCT in clinical and biomedical imaging. To alleviate this problem, a new iterative CT reconstruction algorithm is presented that aims to decrease the radiation dose by reducing the projection views, while maintaining the high quality of reconstructed images. The proposed algorithm combines the adaptive-weighted anisotropic total p-variation (AwaTpV, 0 < p < 1) regularization technique with projection onto convex sets (POCS) strategy. Noteworthy, the AwaTpV regularization term not only contains the horizontal and vertical image gradients but also adds the diagonal image gradients in order to enforce the directional continuity in the gradient domain. To evaluate the effectiveness and ability of the proposed algorithm, experiments with a numerical phantom and synchrotron IL-PCCT were performed, respectively. The results demonstrated that the proposed algorithm had the ability to significantly reduce the artefacts caused by insufficient data and effectively preserved the edge details under noise-free and noisy conditions, and thus could be used as an effective approach to decrease the radiation dose for IL-PCCT.

1. Introduction

In-line X-ray phase-contrast imaging (IL-PCI), also called propagation-based imaging, is a powerful nondestructive imaging technique enabling high contrast between materials with similar attenuation properties. In general, IL-PCI images contain both absorption information and the second derivative (the Laplacian) of phase information, and the latter can lead to strong contrast outlining the structural boundaries of a sample (edge enhancement). Therefore, IL-PCI has a big advantage in enhancing the contrast at boundaries, and this is well suited for imaging biological tissues containing fibre-like structures, such as hepatic fibrous tissue. In general, structural changes of hepatic fibrous tissue are considered as one of the most important pathophysiological features of cirrhosis, and hence hepatic fibrous tissue can be used to accurately characterize cirrhosis (Sethasine et al., 2012). By extending IL-PCI to computed tomography (IL-PCCT), CT images can be produced with higher resolution and better contrast in biological soft tissues compared with traditional CT (termed as absorption-based CT) imaging techniques (Bravin & Coan, 2012), and thus can be used for high-resolution three-dimensional (3D) visualization of the internal detailed structures in weakly absorbing objects (e.g. biological soft tissues). Currently, IL-PCCT has been widely applied to visualize biological soft tissue details and has become one of the most important pre-clinical imaging techniques (Horng et al., 2014; Iyer et al., 2018; Labriet et al., 2018). Typically, the standard filtered back projection (FBP) algorithm is used for CT reconstruction in IL-PCCT. However, due to the Shannon sampling theorem (Jerri, 1977), the FBP requires high completeness of projection data to produce high-quality CT images, which leads to a long total exposure time and thus a high radiation dose. The high radiation dose will not only pose health risks to subjects for in vivo imaging but also has the risk of causing damage to the structures and biological properties of the samples when used for ex vivo imaging (Fernández et al., 2018), which may affect the accurate structure measurements based on the CT images and even influence the subsequent use of ex vivo biological samples (e.g. pathological analysis). It is important for clinical or biomedical applications of IL-PCCT to decrease the radiation dose while maintaining the high quality of reconstructed images. For IL-PCCT, one strategy to decrease the radiation dose is to reduce the amount of sampling data, e.g. few-view projections and limited-angle projections. However, few-view projections and limited-angle projections will cause projection data insufficiency, and thus result in reconstruction artefacts for FBP. To alleviate this problem, many efforts have been directed to developing the CT iterative reconstruction algorithms for IL-PCCT (Zhao et al., 2012, 2018; Melli et al., 2018).

The algebraic reconstruction technique (ART), a classic CT iterative reconstruction method, has the potential to be an alternative approach to FBP, formulating the reconstruction problem as a discrete linear system and enables reconstruction of better results using the incomplete sampling data compared with FBP. However, when the projection data are incomplete, the discrete linear system will be ill-posed, and it may typically be unstable with respect to small perturbations of the noise in projection data (e.g. the low-dose noise). Therefore, it is rare for ART to be independently used for CT reconstruction. To improve the performance of ART, one possible strategy is to combine ART with the total variation (TV) regularization technique, and this strategy enables accurate image recovery from very few measurements according to the compressive sensing (CS) theory (Sidky et al., 2006). In general, the TV-based CT reconstruction algorithm mainly relies on the piecewise-constant assumption of the reconstructed images and the sparsity of the reconstructed images in the discrete gradient transform domain. In this study, for the IL-PCCT images of the human cirrhosis sample stained with iodine, where the iodine was used as contrast agent for contrast-enhanced hepatic fibrous tissue imaging, they can be approximated by a piecewise-constant function and, here, the TV-based reconstruction algorithm can be applicable.

Generally, when using the TV-based reconstruction algorithm, due to the piecewise-constant assumption for the image, the edges of the reconstructed image can cause oversmoothness and sometimes can result in staircase artefacts, which can seriously degrade the quality of the images and influence the subsequent image analysis (e.g. image segmentation, texture analysis and structure measurement). In order to improve this limitation of the TV-based reconstruction algorithm, Liu et al. (2012) proposed an adaptive-weighted TV with the projection data constraints (e.g. data consistency and positivity) enforced by projection onto the convex sets (AwTV-POCS) reconstruction algorithm. In the AwTV-POCS algorithm, the AwTV minimization was performed by the gradient descent method, and the algebraic reconstruction procedure with nonnegative constraint could be considered as the POCS strategy. The AwTV-POCS algorithm can preserve the edge details of the reconstructed image to a certain extent, and this is achieved by considering the anisotropic edge property among neighbouring image pixels, where the associated weights are computed with an edge indicator function that could be adaptively adjusted via the local image gradient magnitude. However, The AwTV term of the AwTV-POCS algorithm is isotropic; namely, the finite difference operators of the AwTV term in the horizontal and vertical directions are unseparated, which may cause the mutual interference between them and thus affect edge detection (Ji et al., 2017). Therefore, limited by the isotropic property of the AwTV term, the AwTV-POCS often has a limited ability to preserve the complex textures and edges like liver images, and it also has a poor capability to handle large-scale artefacts, such as the limited-angle artefacts. Compared with the isotropic TV, the anisotropic TV allows for an improved detection of edge structure information and can also help to handle large-scale artefacts (Chen, Jin, Li & Wang, 2013; Wang, Zhang et al., 2017). However, the traditional anisotropic TV only seeks the gradient sparsity horizontally and vertically, but fails to enforce the directional continuity in the gradient domain. To overcome the deficiency of the traditional anisotropic strategy, Shu et al. proposed that the diagonal image gradients could be incorporated into the original regularization term [further details about the experimental certification are given by Shu & Ahuja (2010) and Wu et al. (2017)]. In general, the total p-variation (TpV, 0 < p < 1) regularization-based reconstruction algorithm can perform better in terms of image quality and reconstruction accuracy compared with the TV-based reconstruction algorithm (Sidky et al., 2007, 2014; Miao & Yu, 2015). Moreover, recent theoretical results showed that values of p leading to nonconvex optimization problems may be practical for CS applications (Chartrand & Staneva, 2010; Chartrand, 2012). This indicates that the exploitation and improvement of the TpV regularization-based reconstruction algorithm will be extremely valuable for practical CT reconstructions.

Inspired by the above-mentioned work, considering both the piecewise-smooth characteristic and the anisotropic edge property of IL-PCCT images of biological tissues (e.g. the hepatic fibrous tissue), a novel adaptive-weighted anisotropic TpV-POCS reconstruction algorithm (AwaTpV-POCS) is proposed. The proposed algorithm combines the AwaTpV regularization technique with the POCS strategy, along with two diagonal gradients adopted to enhance the directional continuity in the gradient domain. To demonstrate the effectiveness and capability of the AwaTpV-POCS algorithm, a numerical liver-cirrhosis (LC) phantom simulation experiment and a synchrotron IL-PCCT experiment on a human cirrhosis sample stained with iodine were carried out.

2. Methods

2.1. IL-PCI and the phase retrieval

In IL-PCI, when the spatially coherent X-ray beams pass through an object, a phase shift in the beam caused by the object will arise. As the beams with the phase shift of the object propagate in the near-field regime, due to Fresnel diffraction, an interference pattern with `Fresnel fringes' is generated, and thus the phase shift can be transformed into variations in intensity, which is finally recorded by a detector set at a specific distance. In practice, the experimental setup of IL-PCI consists of a light source that can provide the spatially coherent X-ray beams, sample stage and detector, and no other optical components are needed (see Fig. 1). Due to its simple setup and the low stability requirements, IL-PCI shows great potential to become clinically applicable.

Figure 1 Schematic illustration (a) and photograph (b) of the IL-PCI experimental setup at BL13W1 at the Shanghai Synchrotron Radiation Facility (SSRF). Here, the quasi-coherent parallel X-ray beams are generated by the 3.5 GeV electron storage ring and are monochromated by a double-crystal monochromator with Si(111) crystals, and then illuminate the sample on the rotation stage. When the transmitted beams propagate from the sample, the density variation causes phase shifts. Because of Fresnel diffraction, the phase shifts are transformed into detectable intensity variations and subsequently recorded by a detector set at a specific distance. Finally, the projection image of IL-PCI is displayed in the image acquisition system. For tomographic scans, the sample can be rotated from 0° to 180° to collect the projection images from various views.

To acquire the phase information distribution of the sample, we need to extract the phase shift from the raw projection images, and thus phase retrieval is required (Nugent et al., 1996). However, to obtain the accurate phase shift, phase retrieval typically requires at least two phase-contrast radiographs, taken at two different sample-to-detector distances (SDDs) (Cloetens et al., 1999). This method is difficult to arrange experimentally, and causes a complicated registration problem and delivers a high radiation dose to the biological samples (Mohammadi et al., 2014). According to Gureyev's simulation study, phase retrieval from one SDD phase-contrast radiograph is possible (Gureyev et al., 2004), and several phase retrieval methods using one SDD phase-contrast radiograph have been proposed, e.g. the modified Bronnikov algorithm (MBA) method (Groso et al., 2006), the phase-attenuation duality Bronnikov algorithm (PAD-BA) (Chen, Rigon & Longo, 2013) etc.

In this study, the PAD-BA algorithm, which is usually useful for imaging objects containing components with similar X-ray attenuation coefficients, was implemented on the projection images using PITRE software to extract quantitative phase-shift information (Chen et al., 2012). The PAD-BA algorithm is grounded in two assumptions on the material properties: (i) that the imaging object is weakly absorbing and quasi­homogeneous; and (ii) the PAD property, i.e. the δ and β parts of the complex refractive index n = 1 − δ + iβ are proportional to each other; here the β value is related to the absorption information, whereas the δ value is related to the phase information. In this experiment, after some experimental trials, it was found that the reconstructed result using a δ/β value of 400 had a high contrast between the adjacent tissues (e.g. the hepatic fibrous tissue and the liver parenchyma) and enabled the optimal visual effect for edge details and fine features. Therefore, 400 was considered as the best δ/β value. After phase retrieval, the phase-shift distribution from the IL-PCI of samples was acquired, and the corresponding IL-PCCT images could then be obtained by the CT reconstruction.

2.2. The CT iterative reconstruction method for IL-PCCT

Following phase retrieval, taking into consideration that the synchrotron radiation is a monochromatic source, the general model of the IL-PCCT imaging can be approximately expressed as the following discrete linear system,

where is the projection coefficient matrix that represents the parallel-beam X-ray forward projection (Siddon, 1985). The projection data acquired from the detector are given by , and denotes the phase information distribution of the illuminated object. The goal of IL-PCCT reconstruction is to accurately reconstruct u from g. Mathematically speaking, equation (1) is an ill-posed inverse problem for the case that the projection data g is insufficient (e.g. few-view projections and limited-angle projections), which means that the unique solution for u cannot be obtained by directly inverting equation (1). To solve the discrete linear system represented in equation (1), an optimization-based strategy is often adopted.

In this study, to find the optimal solution u^* to the problem in equation (1), we develop an optimization-based method, and its cost function can be formulated as follows,

where is the data fidelity term, is the AwaTpV regularization term, and λ is a regularization parameter that controls the trade-off between the data fidelity term and regularization term. = f_i,j - f_{i - 1,j} and = f_i,j - f_{i,j - 1} are the horizontal and vertical gradients, respectively, which measure the gradient sparsity. = f_i,j - f_{i - 1,j - 1} and = f_{i,j - 1} - f_{i - 1,j} are two diagonal partial gradients, which measure the gradient continuity. w₁, w₂, w₃ and w₄ are weights added on to , , and , respectively, and can be adaptively calculated using an edge indicator function, namely w₁ = , w₂ = , w₃ = and w₄ = . Here, c and σ in weights are scale factors which control the strength of the diffusion. In this study, the pth power of Lp norm was defined as Lp-quasinorm, i.e. := , and the Lp-quasinorm of the image gradients was formulated as a TpV function.

According to the proximal forward–backward splitting optimization algorithm (Gao, 2016), the optimization problem in equation (2) can be transformed into a few subproblems which are calculated for the optimal solution of the original optimization problem in an alternately iterative manner. Here, the equivalent form of the optimization problem in equation (2) can be formulated as the following two subproblems,

where the step for calculating z^t is a gradient descent update with a step size of 1/a_t for the problem

the superscript T represents the transpose operator; the superscripts t and t+1 denote the tth iteration and (t+1)th iteration, respectively.

2.2.1. The POCS strategy

To solve the subproblem in equation (3), the simultaneous algebraic reconstruction technique (SART), a major refinement of the ART, could be utilized herein (Andersen & Kak, 1984), and this procedure with nonnegative constraint was considered as the POCS strategy, which is similar to the POCS implementation of the AwTV-POCS algorithm (Liu et al., 2012). In this work, a block-iterative SART technique is adopted (Hansen & Saxild-Hansen, 2012), which has the potential to handle large-scale linear inverse problems quickly, and is expressed as follows,

where represents the relaxation parameter of the tth iteration, and V and W are the diagonal matrixes with row sums and column sums of A in the diagonal, respectively.

Considering that each component of z ^t is nonnegative, the nonnegative constraint for the block-iterative SART is here imposed. Let z_j ^t be defined as a component of z ^t, and the nonnegative constraint can be written as follows,

To improve the convergence performance of the SART algorithm, the relaxation parameter in each iteration is chosen using the line search scheme, which can be computed as follows,

where denotes the square of the L₂ norm.

2.2.2. The AwaTpV regularization technique

The subproblem in equation (4) represents the AwaTpV regularization problem. Although it usually leads to a nonconvex optimization problem and can obtain only local optimal solutions, it can yield solutions sparser than the TV regularization (Sidky et al., 2007; Chartrand, 2007). Moreover, some studies suggested that the nonconvex optimization problems may be more accurate for practical image inverse problems compared with convex optimization problems, because the latter often resulted in a biased sparsity approximation and inaccurate reconstruction under some practical problems (Wang, Zhang et al., 2017). According to previous work (Pan et al., 2013), the new auxiliary variables corresponding to the weighted image gradient magnitude were introduced, i.e. d_n = , n = 1, 2, 3, 4, and the equivalent form of the subproblem in equation (4) could be written as follows,

where β is a positive parameter that constraints the auxiliary variables d_n close to their corresponding gradients ; denotes the ratio of the parameter λ and a_t.

In this study, to solve the subproblem in equation (8), the split Bregman iteration (SBI) algorithm (Goldstein & Osher, 2009) was implemented, and the dual variables with respect to d_n,n = 1, 2, 3, 4, were introduced, i.e. b_n, n = 1, 2, 3, 4. Then, the subproblem in equation (8) can be split into the following sub-optimization problems which are then calculated in the manner of alternating iteration,

where k denotes the number of iteration for AwaTpV minimization.

(i) u-subproblem. The objective function in equation (9) is a quadratic function that can be solved by a conventional gradient decent approach. In this study, the fast Fourier transform (FFT) method was implemented to solve the sub-optimization problem in equation (9), which can produce the solution in closed form of the above-mentioned sub-optimization problem and enables acceleration of convergence according to previous work (Wang et al., 2008). Then, the solution to u could be obtained as follows,

where and represent the 2D FFT and its inverse operators, respectively; = denote the horizontals, vertical and two diagonal image gradient operators; the superscript * denotes the complex conjugate, the symbol 1 represents a matrix where all the entries are 1, and the symbol is the component-wise multiplication operator; all operators here (such as addition, multiplication and division) are component-wise.

(ii) d_n-subproblem. To solve the sub-optimization problem in equation (10), an iterative p-shrinkage (IPS) method was herein utilized (Woodworth & Chartrand, 2016). This method could produce the exact solution to the sub-optimization problem in equation (10), and it was also guaranteed to converge (Zuo et al., 2013). Here, the IPS method was defined as follows,

where denotes the IPS operator; τ is the threshold; , when p = 1, the IPS degrades as the standard soft thresholding shrinkage (Daubechies et al., 2004).

Therefore, the solution of the sub-optimization problems in equation (10) could be written as follows,

2.3. Pseudocode for the AwaTpV-POCS algorithm

In order to explain the AwaTpV-POCS algorithm more clearly, the corresponding pseudocode was provided, as illustrated in Algorithm 1.

.

2.4. Low-dose noise model in projections

To analyze the robustness and reliability of the AwaTpV-POCS method in the presence of low-dose noise, the low-dose noise was added to the projections and then the corresponding CT reconstructions were carried out. In this experiment, the low-dose noise in projections was defined as a combination of the Poisson noise and Gaussian noise (Liu et al., 2012), and its mathematical model could be formulated as depicted in equation (15),

where represents the simulated noisy intensity measurement of detector unit i at a projection view and I₀ denotes the incoming X-ray intensity, is defined as the logarithmic transform of , m_ie and are the mean value and variance of the Gaussian noise, for detector unit i. In this experiment, the photon number for Poisson noise I₀ was set to 1.0 × 10⁵, while m_ie and were set to 0 and 10, respectively.

2.5. Parameter selection

In this study, there are five parameters involved in the AwaTpV-POCS method: p, β, λ^*, c and σ. These parameters have significant impact on the reconstructed image, and the suitable selection for parameters can guarantee a favourable reconstructed result. However, selecting optimal parameters is a common difficulty for multi-parameter algorithms. In general, the best approach for determining the optimal values of the above-mentioned five parameters is to find the global optimal value of peak signal-to-noise ratio (PSNR) in the five-dimensional parameter space, but this method requires a lot of computation cost and thus is difficult to implement in reality. To alleviate this problem, a greedy strategy to determine the parameters one by one was used (Lohvithee et al., 2017). Although this strategy may acquire a local optimum, it gives an approximately optimal result. In this experiment, the greedy strategy was implemented to select the optimal parameters. For the selection of the p value in the AwaTpV regularization, the different p values were tuned while keeping other parameters fixed to produce different image reconstructions, and the optimal p value could be determined by the highest PSNR value based on the corresponding image reconstructions. The processes for the selection of p value are shown later in Fig. 6(d). The other parameters can also be determined in the above-mentioned way, including the parameter α that was defined as the regularization parameter of the AwTV-POCS method. In the synchrotron IL-PCCT experiment, the contrast-to-noise ratio (CNR) was used as the evaluation measure. The parameters for the simulation and synchrotron IL-PCCT experiment will be subsequently shown in Tables 1 and 3, respectively.

2.6. Quantitative assessment

In this study, the PSNR, structural similarity index (SSIM), relative error (RE), CNR and relative difference (RD) were adopted as quantitative metrics. The PSNR is a traditional measure for image quality, and a larger PSNR value represents better quality. The SSIM is usually used to assess the similarity between the reconstructed and reference image, which yields a value between 0 and 1 that increases with increasing similarity (Wang et al., 2004). The RE can be used to evaluate the reconstruction error, and a smaller RE value indicates more accuracy. The CNR can be adopted to evaluate contrast between the adjacent tissues, and a larger CNR value represents higher contrast (Zeller-Plumhoff et al., 2017). The RD can be used to measure the difference between two reconstructed images of the same reconstruction method in different iterations, and a greater RD value means a larger difference.

(i) PSNR is widely used and defined as follows,

where MSE is the mean square error function, x denotes the reference image, y denotes the reconstructed image, and M ×N is the size of x and y; x_i, j and y_i, j, respectively, represent the pixel intensity of x and y in some pixel (i,j), and Peak represents the largest pixel intensity in the normalized image; in our study, Peak is 255.

(ii) SSIM is defined as follows,

where u_x, u_y are the mean values of x and y, respectively; , denote the variances of x and y, respectively; represents the standard deviation between x and y; c₁ and c₂ are two constants to stabilize the division with weak denominator, i.e. c₁ = (k₁L)², c₂ = (k₂L)²; in our study, L was set to 255, k₁ and k₂ were set to 0.01 and 0.03, respectively.

(iii) RE is defined as follows,

(iv) CNR is defined as follows,

where y₁ and y₂ represent the local pixel regions of the two adjacent tissues in the same reconstructed image, u_y1 and u_y2 denote the mean value of y₁ and y₂, respectively, and  are the variances of y₁ and y₂, respectively.

(v) RD is defined as follows,

where y_t and y_t+1 represent the reconstructed images in tth and (t+1)th iterations, respectively.

3. Simulation experiment

3.1. Simulations

To evaluate the performance of the AwaTpV-POCS algorithm, a numerical LC phantom designed via the pathological section of the fibrous tissue in human cirrhosis sample was employed, as shown in Fig. 2. The LC phantom is a grayscale image with the grey value ranging from 0 to 255 and an overall size of 512 × 512 pixels, which was used to simulate the phase information distribution of the sample. The detector is modelled as a straight-line array with 724 bins, and the size of the reconstructed images is 512 × 512 pixels. The 60 uniformly distributed projections over 180° (full-scan range) were used to simulate few-view projections. In addition, the 60 uniformly distributed noise-free and noisy projections scanned from 30° to 120° were used to simulate limited-angle projections and low-dose noisy limited-angle projections, respectively. The low-dose noise model added into limited-angle projections were introduced in detail in §2.4. The total iteration number was selected as the stopping criterion and could be determined according to the convergence curves, as shown in Fig. 6. Here, the FBP, SART and AwTV-POCS algorithms were used for comparison with the AwaTpV-POCS algorithm, and all parameters were chosen for optimal performance, as shown in Table 1. All experiments were carried out using the MATLAB language on a personal computer equipped with Intel(R) Core(TM) i5-4460 CPU at 3.2 GHz and 16 GB RAM.

Table 1
Parameters for simulation

	Method	p	β	λ*	α	c	σ
Few-view	AwTV-POCS	–	–	–	0.2	0.6	15
	AwaTpV-POCS	0.2	0.8	0.008	–	0.6	15
Limited-angle	AwTV-POCS	–	–	–	0.3	0.7	25
	AwaTpV-POCS	0.2	0.5	0.01	–	0.6	15
Low-dose noisy limited-angle	AwTV-POCS	–	–	–	0.7	0.6	15
	AwaTpV-POCS	0.8	0.5	0.03	–	0.7	25

Figure 2 The LC phantom designed for testing CT reconstruction algorithms. (a) The ground truth. (b) Enlarged image of green rectangle region in (a).

3.2. Experimental results

In this section, the simulation experiments were performed under three conditions of few-view projections, limited-angle projections and low-dose noisy limited-angle projections. The experimental results of the reconstructed images using different CT reconstruction algorithms under different conditions are given, as shown in Fig. 3. Among the reconstructed images, those reconstructed using the FBP algorithm contain the most severe artefacts due to the incomplete projections or low-dose noise, as shown in Figs. 3(a,e,i). It can be seen in Figs. 3(b,f,j) that there are still undesirable artefacts in the SART results. In Figs. 3(c,g,k), it can be observed that AwTV-POCS can suppress most of the streak artefacts; however, there are still some artefacts in Figs. 3(g) and 3(k), which were introduced by limited-angle projections or low-dose noise. AwaTpV-POCS achieves the best visual effect in Figs. 3(d,h,i), which show that AwaTpV-POCS eliminates all the streak artefacts, and also eables to better preserve the edge details and suppresses artefacts from limited-angle projections and low-dose noise effectively. The regions of interest (ROIs) were selected and enlarged to improve the visual effect for better comparison, as shown in Fig. 4.

Figure 3 Reconstructed images of the LC phantom using the FBP (a,e,i), SART (b,f,j), AwTV-POCS (c,g,k) and AwaTpV-POCS (d,h,l) algorithms from 60-view projections under few-view, limited-angle and low-dose noisy limited-angle conditions, respectively. Panels (a)–(d), (e)–(h) and (i)–(l) are reconstructed images under few-view, limited-angle and low-dose noisy limited-angle conditions, respectively. The blue arrow denotes streak artefacts; the green arrow denotes limited-angle artefacts. The pixel values of the above greyscale images were normalized to the range [0, 255]. The display window is [0, 255].

Figure 4 Enlarged regions of the reconstructed images of the LC phantom using the FBP (a,e,i), SART (b,f,j), AwTV-POCS (c,g,k) and AwaTpV-POCS (d,h,l) algorithms from 60-view projections under few-view, limited-angle and low-dose noisy limited-angle conditions, respectively. The magnified images were from the same regions in the reconstructed images, as marked with the green rectangle in Fig. 2(a). The pixel values of the above greyscale images were normalized to the range [0, 255]. The display window is [0, 255].

3.3. Assessments

To assess the accuracy of the four reconstruction approaches under few-view, limited-angle and low-dose noisy limited-angle conditions, the profiles of the same position in the reconstructed images, the position marked with the red line in Fig. 2(a), were utilized, as shown in Fig. 5. It can be easily observed that the profiles of the AwaTpV-POCS approach are closest to the true result under the three aforementioned conditions. In addition, the PSNR, SSIM and RE values of the reconstructed images (Fig. 3) were further calculated (Table 2), and the results show that the quality of the images reconstructed using the AwaTpV-POCS approach is obviously the best.

Table 2
Quality metrics of reconstructed LC phantom images

	Method	PSNR (dB)	SSIM	RE (%)
Few-view	FBP	19.2982	0.3444	29.4000
	SART	26.2997	0.7167	6.2800
	AwTV-POCS	28.0770	0.8887	2.6500
	AwaTpV-POCS	30.5168	0.9268	1.9700
Limited-angle	FBP	18.0591	0.2918	35.4197
	SART	22.6283	0.5902	13.9500
	AwTV-POCS	24.1052	0.7917	10.8306
	AwaTpV-POCS	25.1669	0.8259	7.6684
Low-dose noisy	FBP	16.8145	0.2519	35.3500
limited-angle	SART	22.0546	0.5166	14.0900
	AwTV-POCS	23.0263	0.7026	12.8900
	AwaTpV-POCS	23.8013	0.7314	10.2000

Figure 5 Profiles of the LC phantom images reconstructed by different methods from 60 projections under few-view (a), limited-angle (b) and low-dose noisy limited-angle (c) conditions. The profiles were located at the pixel position labelled with the red line in Fig. 1(a).

To evaluate the convergence performance of the SART, AwTV-POCS and AwaTpV-POCS approaches under the few-view, limited-angle and low-dose noisy limited-angle conditions, the PSNR-based and RE-based convergence curves of the above-mentioned approaches are presented, as shown in Figs. 6(a)–6(c). It can be observed that the SART, AwTV-POCS and AwaTpV-POCS algorithms converged before the iterations reach 50 under the few-view condition, but converged before the iterations reach 300 under other two conditions and, obviously, the convergence performance of the AwaTpV-POCS approach is the best.

Figure 6 Convergence curves of the SART, AwTV-POCS and AwaTpV-POCS approaches under few-view (a), limited-angle (b) and low-dose noisy limited-angle (c) conditions, and the performance of AwaTpV-POCS with respect to different values of p (d). In panels (a)–(c), the left-hand axes plot PSNR with respect to different iteration number, and the right-hand plot RE with respect to different iteration number. The blue arrows and the red arrows mark PSNR-based and RE-based convergence curves, respectively. In panel (d), the red line with rectangles, the green line with diamonds and the blue line with triangles represent the performance of AwaTpV-POCS with respect to different values of p under few-view, limited-angle and low-dose noisy limited-angle conditions, respectively.

4. Real experiment on synchrotron IL-PCCT data

4.1. Data acquisition

This study was approved by the Research Ethics Committee of Beijing Friendship Hospital, Capital Medical University, Beijing, China, and written informed consent was obtained from all patients. In this study, the experimental ex vivo human cirrhosis sample was provided by the Beijing Friendship Hospital, and its histopathologic section and IL-PCCT image are presented in Fig. 7. In this experiment, the iodine was used as the contrast agent for the contrast-enhanced imaging of fibrous tissue in the cirrhosis sample, and its IL-PCCT data were acquired at the BL13W1 beamline at SSRF, Shanghai, China. In this experiment, the energy of the incoming monochromatic photons was adjusted to 33 keV, and the SDD was set to 0.8 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 used as the imaging detector, with an effective pixel size of 9 µm × 9 µm. The full projection dataset (930-view projections) within a 180° CT scan range was collected with an exposure time of 12 ms per projection and a projection image size of 3992 × 513 pixels. In addition, ten dark-current images (the offset signal recorded when no photons hit the detector) were used to remove the dark signal in projections (dark-field correction), while 20 flat-field images (the images with no sample in the beam) were used to correct the pixel-to-pixel nonuniformity in projections caused by beam's intensity inhomogeneity or non-uniform detector response (flat-field correction) (Chen et al., 2012). After phase retrieval using the PAD-BA method, the 186-view projections were uniformly chosen from the full projection dataset and the 414-view projections were uniformly chosen from the projection dataset with respect to CT scan ranging from 10° to 170°. A sinogram with 186 × 3992 pixels and a sinogram with 414 × 3992 pixels could then be generated for few-view CT reconstruction and limited-angle CT reconstruction, respectively. In this experiment, all parameters were chosen for optimal performance, as shown in Table 3.

Table 3
Parameters for the synchrotron IL-PCCT experiment

	Method	p	β	λ*	α	c	σ
Few-view	AwTV-POCS	–	–	–	0.9	0.6	15
	AwaTpV-POCS	0.1	0.2	0. 8	–	0.6	15
Limited-angle	AwTV-POCS	–	–	–	1.2	0.7	25
	AwaTpV-POCS	0.3	0.15	0.7	–	0.6	15

Figure 7 Histopathologic section stained with Sirius red of a human cirrhosis sample (a) and its IL-PCCT image (b). Here, image (b) is the reconstructed result using the FBP algorithm based on the full projection dataset, and it can be utilized as the reference image for reconstructions. Panels (a) and (d) are enlarged images of the green rectangle regions in (a) and (b), respectively. The black arrows mark the fibrous tissue in the human cirrhosis sample. The scale bar represents 500 µm. It can be seen that the IL-PCCT image of hepatic fibrous tissue shows good agreement with its pathological section.

4.2. Experimental results

Fig. 8 shows reconstructed images of the human cirrhosis sample using the FBP, SART, AwTV-POCS and AwaTpV-POCS algorithms. Figs. 8(a)–8(d) are a reconstructed slice of human cirrhosis sample from 186-view projections under few-view conditions using the FBP, SART, AwTV-POCS and AwaTpV-POCS algorithms, respectively. Figs. 8(i)–8(l) are reconstructed slices of human cirrhosis sample with 414-view projections under limited-angle condition using the FBP, SART, AwTV-POCS and AwaTpV-POCS algorithms, respectively. Figs. 8(a) and 8(i) show that the reconstructed slices using FBP have poor image quality and, where the detail structures and edges are severely affected by artefacts, the subsequent image segmentation and image analysis may be influenced. Figs. 8(b) and 8(j) indicate that the SART algorithm has the limited ability to deal with the few-view CT reconstruction and the limited-angle CT reconstruction problems. From Figs. 8(c) and 8(k), it can be clearly observed that the image quality of the reconstructed images using the AwTV-POCS has been improved significantly in comparison with that from the FBP and SART algorithms, and this suggests that the AwTV-POCS algorithms have good capability when it comes to the few-view CT reconstruction and the limited-angle CT reconstruction problems. However, there are a few residual artefacts in the above-mentioned results. Figs. 8(d) and 8(l) depict the reconstructed results using the AwaTpV-POCS algorithm, and it can be easily observed that the results of AwaTpV-POCS have the fewest artefacts and the clearest edge details of the hepatic fibrous tissue, and the best image quality compared with those of other algorithms.

Figure 8 Reconstructed images of the human cirrhosis sample using the FBP (a) and (i), SART (b) and (j), AwTV-POCS (c) and (k), and AwaTpV-POCS (d) and (l) algorithms from 186-view projections under few-view condition and 414-view projections under limited-angle condition, respectively. Panels (e)–(h) are enlarged images in reconstructed images (a)–(d), respectively. Panels (m)–(p) are enlarged images in reconstructed images (i)–(l), respectively. The enlarged images were from the same regions in the reconstructed images, as marked with the green rectangle in Fig. 7(b). The blue rectangle with the size of 80 × 50 marks the fibrous tissue in cirrhosis sample, and the yellow rectangle with the size of 80 × 50 marks the liver parenchyma tissue in cirrhosis sample. The black arrows mark the fibrous tissue in cirrhosis sample. The pixel values of the above greyscale images were normalized to the range [0, 255]. The display window is [25, 225].

4.3. Result analysis

To evaluate the accuracy of the CT images (Fig. 8) reconstructed using the FBP, SART, AwTV-POCS and AwaTpV-POCS algorithms under few-view and limited-angle conditions, the position labelled with the red line in Fig. 8(a), which crossed through the fibrous tissue region, was used. The profiles of the above-mentioned positions in Figs. 8(a)–8(d) and 8(i)–8(l) are shown in Figs. 9(a) and Fig. 9(b), respectively. In this experiment, Fig. 8(b) was chosen as the reference image, and hence its profile was chosen as the reference line. As shown in Figs. 9(a) and Fig. 9(b), the profiles and intensities of the CT images reconstructed using FBP and SART have a large deviation from the reference image, which may lead to image distortion and thus affect the analysis and judgement of the doctors or researchers. Additionally, the profiles and intensities of the CT images reconstructed using AwTV-POCS and AwaTpV-POCS are close to the reference image, indicating that AwTV-POCS and AwaTpV-POCS have high accuracy in practical CT reconstruction. However, when compared with the profile of the CT images reconstructed using AwTV-POCS, the profile of the CT images reconstructed using AwaTpV-POCS is smoother in the fibrous tissue region and, considering the homogeneous characteristic of the fibrous tissue region, this demonstrates that the AwaTpV-POCS performs better to remove artefacts in the fibrous tissue region. In addition, to quantitatively evaluate the contrast between the fibrous tissue and the liver parenchyma, the CNR values of different reconstructed images were measured [shown in Fig. 9(c)] using a pair of non-overlapping ROIs marked with the blue and yellow boxes in Fig. 8(a), and the CNR value of the reference image was chosen as the reference value. From Fig. 9(c), it can be observed that the CNR values of AwaTpV-POCS in few-view and limited-angle cases are highest; namely, the contrast between the fibrous tissue and the liver parenchyma in the reconstructed images of AwaTpV-POCS is highest, and this helps to offer the best visual effect for distinguishing the fibrous tissue from the liver parenchyma. Therefore, the fibrous tissue in the reconstructed images of AwaTpV-POCS is not only accurate but has the fewest artefacts and the best visual effect. This will help to accurately measure and analyze the fibrous tissue for evaluating the degree of the liver fibrosis or cirrhosis.

Figure 9 Profiles of CT images of the human cirrhosis sample reconstructed using the FBP, SART, AwTV-POCS and AwaTpV-POCS algorithms under few-view (a) and limited-angle (b) conditions, and their corresponding CNR results (c).

To test the convergence performance of the SART, AwTV-POCS and AwaTpV-POCS algorithms in practical CT reconstruction under few-view and limited-angle conditions, the RD-based convergence curves of the above-mentioned three methods are drawn, as shown in Figs. 10(a) and 10(b), respectively. It can be clearly observed that the SART, AwTV-POCS and AwaTpV-POCS algorithms converged before the iterations reach 30 under few-view and limited-angle conditions, demonstrating that the solutions to the above-mentioned three methods in practical CT reconstruction under few-view and limited-angle conditions can converge.

Figure 10 Convergence curves of the SART, AwTV-POCS and AwaTpV-POCS algorithms under few-view (a) and limited-angle (b) conditions. Panels (c) and (d) are enlarged images of the regions depicted by the red rectangles in (a) and (b), respectively. For the real experiment, the total iteration number of 30 was selected as the stopping criterion.

5. Discussion and conclusion

In this study, an AwaTpV-POCS algorithm was proposed for accurate CT reconstruction based on the few-view and limited-angle projections. In the experiment, the LC phantom and the synchrotron IL-PCCT data of the human cirrhosis sample were used to evaluate the accuracy and the feasibility of this algorithm, and the FBP, SART and AwTV-POCS algorithms were utilized for comparison. The results confirmed that the AwaTpV-POCS algorithm was an effective approach to decrease the radiation dose for IL-PCCT, which had the ability to significantly reduce the artefacts introduced by the few-view and limited-angle sampling, and enabled the edge details to be preserved and the low-dose noise to be suppressed.

Although the AwaTpV-POCS algorithm is a nonconvex optimization method, which can obtain only local optimal solutions in general, it can yield more accurate solutions than the AwTV-POCS regularization with the guaranteed convergence to the local optimal solution. Recently, with the rapid development of the nonconvex optimization method in image processing, its medical application has drawn widespread attention (Zhang et al., 2018; Chen & Anastasio, 2018). To our knowledge, it was the first to implement the iterative CT reconstruction algorithm based on nonconvex optimization to the IL-PCCT, and the accuracy and the feasibility of this algorithm is demonstrated in this research. Taking into account that the small difference in X-ray attenuation coefficients between the adjacent soft tissues still constitutes a significant problem for diagnostic interpretation, a limitation for absorption-based CT in the soft tissues, the phase-based IL-PCCT has outstanding performance in revealing small variations inside soft tissues. It therefore has the potential to be an alternative method to conventional CT for imaging the soft tissues. Recently, IL-PCI experiments carried out on conventional X-ray sources have demonstrated that comparable image quality could be obtained from the benchtop imaging system (Zysk et al., 2012; Larsson et al., 2016). Moreover, the IL-PCI experiments utilizing live animal subjects have made great headway (Wang et al., 2013; Preissner et al., 2018). These progresses may pave the way for the realization of preclinical or clinical IL-PCI systems. However, many obstacles still remain for the realization of preclinical or clinical IL-PCI systems, e.g. the small FOV, motion artefacts due to heartbeat and breathing, high radiation dose, etc. Corresponding solutions to the above problems are currently being researched. This current work aims to reduce the radiation dose of IL-PCCT while maintaining excellent image quality using newly iterative CT reconstruction algorithms. In this work, a high-quality slice of the human cirrhosis sample stained with iodine was reconstructed using the AwaTpV-POCS algorithm from the highly undersampled IL-PCCT projections, indicating that the proposed method is a valuable tool for low-dose CT reconstruction of the fibrous tissue and can then be used to evaluate the degree of the liver fibrosis or cirrhosis. It is particularly worth noting that the proposed method also enables reconstruction of other IL-PCCT images and other modal CT images under the assumptions that they are piecewise-smooth. In further research, the graphic processing unit (GPU)-based parallel computing techniques will be implemented to increase the computing speed of the AwaTpV-POCS algorithm and, additionally, the adaptivity of the parameters in the AwaTpV-POCS algorithm will be researched. Furthermore, future studies will be carried out to assess whether the AwaTpV-POCS algorithm also applies for in vivo data, along with testing to employ more non-convex iterative CT reconstruction algorithms to the IL-PCI and other PCI methods, such as the speckle-based PCI technique (Bérujon et al., 2012; Bérujon & Ziegler, 2016).

In summary, the main features and contributions of the present work are summarized as follows: (i) hepatic fibrous tissue is clearly visualized in the IL-PCCT image by means of iodine contrast agent, and it shows good agreement with its pathological section; (ii) a new nonconvex IL-PCCT reconstruction algorithm that eables the preservation of the edge details for biological tissues is presented, and it can be applicable for insufficient projection data; (iii) the proposed algorithm can be potentially used to decrease radiation dose for IL-PCCT, and promote the wider use of IL-PCCT in clinical and biomedical imaging.