research papers
A new inline Xray phasecontrast computed tomography reconstruction algorithm based on adaptiveweighted anisotropic TpV regularization for insufficient data
^{a}School of Biomedical Engineering and Technology, Tianjin Medical University, Tianjin 300070, People's Republic of China, ^{b}The School of Science, Tianjin University of Technology and Education, Tianjin 300222, People's Republic of China, ^{c}School of Physics and Information Engineering, Minnan Normal University, 363000 Fujian, People's Republic of China, ^{d}School of Information and Communication Engineering, University of Electronic Science and Technology of China, 610054 Sichuan, People's Republic of China, ^{e}Radiation Oncology Department, Tianjin Medical University General Hospital, Tianjin 300070, People's Republic of China, and ^{f}Liver Research Center, Beijing Friendship Hospital, Capital Medical University, 100050 Beijing, People's Republic of China
^{*}Correspondence email: chunhong_hu@hotmail.com
Inline Xray phasecontrast computed tomography (ILPCCT) 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 ILPCCT will result in a high radiation dose to biological samples, and thus impede the wider use of ILPCCT 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 adaptiveweighted anisotropic total pvariation (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 ILPCCT 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 noisefree and noisy conditions, and thus could be used as an effective approach to decrease the radiation dose for ILPCCT.
Keywords: inline Xray phasecontrast computed tomography; adaptiveweighted anisotropic total pvariation regularization; projection onto convex sets; radiation dose.
1. Introduction
Inline Xray phasecontrast imaging (ILPCI), also called propagationbased imaging, is a powerful nondestructive imaging technique enabling high contrast between materials with similar attenuation properties. In general, ILPCI 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, ILPCI has a big advantage in enhancing the contrast at boundaries, and this is well suited for imaging biological tissues containing fibrelike 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 ILPCI to computed tomography (ILPCCT), CT images can be produced with higher resolution and better contrast in biological soft tissues compared with traditional CT (termed as absorptionbased CT) imaging techniques (Bravin & Coan, 2012), and thus can be used for highresolution threedimensional (3D) visualization of the internal detailed structures in weakly absorbing objects (e.g. biological soft tissues). Currently, ILPCCT has been widely applied to visualize biological soft tissue details and has become one of the most important preclinical 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 ILPCCT. However, due to the Shannon sampling theorem (Jerri, 1977), the FBP requires high completeness of projection data to produce highquality 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 ILPCCT to decrease the radiation dose while maintaining the high quality of reconstructed images. For ILPCCT, one strategy to decrease the radiation dose is to reduce the amount of sampling data, e.g. fewview projections and limitedangle projections. However, fewview projections and limitedangle 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 ILPCCT (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 illposed, and it may typically be unstable with respect to small perturbations of the noise in projection data (e.g. the lowdose 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 TVbased CT reconstruction algorithm mainly relies on the piecewiseconstant assumption of the reconstructed images and the sparsity of the reconstructed images in the discrete gradient transform domain. In this study, for the ILPCCT images of the human cirrhosis sample stained with iodine, where the iodine was used as contrast agent for contrastenhanced hepatic fibrous tissue imaging, they can be approximated by a piecewiseconstant function and, here, the TVbased reconstruction algorithm can be applicable.
Generally, when using the TVbased reconstruction algorithm, due to the piecewiseconstant 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 TVbased reconstruction algorithm, Liu et al. (2012) proposed an adaptiveweighted TV with the projection data constraints (e.g. data consistency and positivity) enforced by projection onto the convex sets (AwTVPOCS) reconstruction algorithm. In the AwTVPOCS 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 AwTVPOCS 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 AwTVPOCS 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 AwTVPOCS often has a limited ability to preserve the complex textures and edges like liver images, and it also has a poor capability to handle largescale artefacts, such as the limitedangle artefacts. Compared with the isotropic TV, the anisotropic TV allows for an improved detection of edge structure information and can also help to handle largescale 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 pvariation (TpV, 0 < p < 1) regularizationbased reconstruction algorithm can perform better in terms of image quality and reconstruction accuracy compared with the TVbased 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 regularizationbased reconstruction algorithm will be extremely valuable for practical CT reconstructions.
Inspired by the abovementioned work, considering both the piecewisesmooth characteristic and the anisotropic edge property of ILPCCT images of biological tissues (e.g. the hepatic fibrous tissue), a novel adaptiveweighted anisotropic TpVPOCS reconstruction algorithm (AwaTpVPOCS) 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 AwaTpVPOCS algorithm, a numerical livercirrhosis (LC) phantom simulation experiment and a synchrotron ILPCCT experiment on a human cirrhosis sample stained with iodine were carried out.
2. Methods
2.1. ILPCI and the phase retrieval
In ILPCI, when the spatially coherent Xray 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 nearfield 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 ILPCI consists of a light source that can provide the spatially coherent Xray 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, ILPCI shows great potential to become clinically applicable.
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 phasecontrast radiographs, taken at two different sampletodetector 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 phasecontrast is possible (Gureyev et al., 2004), and several phase retrieval methods using one SDD phasecontrast have been proposed, e.g. the modified Bronnikov algorithm (MBA) method (Groso et al., 2006), the phaseattenuation duality Bronnikov algorithm (PADBA) (Chen, Rigon & Longo, 2013) etc.
In this study, the PADBA algorithm, which is usually useful for imaging objects containing components with similar Xray attenuation coefficients, was implemented on the projection images using PITRE software to extract quantitative phaseshift information (Chen et al., 2012). The PADBA algorithm is grounded in two assumptions on the material properties: (i) that the imaging object is weakly absorbing and quasihomogeneous; and (ii) the PAD property, i.e. the δ and β parts of the complex 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 phaseshift distribution from the ILPCI of samples was acquired, and the corresponding ILPCCT images could then be obtained by the CT reconstruction.
2.2. The CT iterative reconstruction method for ILPCCT
Following phase retrieval, taking into consideration that the synchrotron radiation is a monochromatic source, the general model of the ILPCCT imaging can be approximately expressed as the following discrete linear system,
where is the projection coefficient matrix that represents the parallelbeam Xray 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 ILPCCT reconstruction is to accurately reconstruct u from g. Mathematically speaking, equation (1) is an illposed inverse problem for the case that the projection data g is insufficient (e.g. fewview projections and limitedangle 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 optimizationbased strategy is often adopted.
In this study, to find the optimal solution u^{*} to the problem in equation (1), we develop an optimizationbased 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 tradeoff 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_{1}, w_{2}, w_{3} and w_{4} are weights added on to , , and , respectively, and can be adaptively calculated using an edge indicator function, namely w_{1} = , w_{2} = , w_{3} = and w_{4} = . 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 Lpquasinorm, i.e. := , and the Lpquasinorm 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 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 AwTVPOCS algorithm (Liu et al., 2012). In this work, a blockiterative SART technique is adopted (Hansen & SaxildHansen, 2012), which has the potential to handle largescale 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 blockiterative 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_{2} 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 suboptimization problems which are then calculated in the manner of alternating iteration,
where k denotes the number of iteration for AwaTpV minimization.
(i) usubproblem. 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 suboptimization problem in equation (9), which can produce the solution in closed form of the abovementioned suboptimization 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 componentwise multiplication operator; all operators here (such as addition, multiplication and division) are componentwise.
(ii) d_{n}subproblem. To solve the suboptimization problem in equation (10), an iterative pshrinkage (IPS) method was herein utilized (Woodworth & Chartrand, 2016). This method could produce the exact solution to the suboptimization 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 suboptimization problems in equation (10) could be written as follows,
2.3. Pseudocode for the AwaTpVPOCS algorithm
In order to explain the AwaTpVPOCS algorithm more clearly, the corresponding pseudocode was provided, as illustrated in Algorithm 1.
.2.4. Lowdose noise model in projections
To analyze the robustness and reliability of the AwaTpVPOCS method in the presence of lowdose noise, the lowdose noise was added to the projections and then the corresponding CT reconstructions were carried out. In this experiment, the lowdose 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_{0} denotes the incoming 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_{0} was set to 1.0 × 10^{5}, 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 AwaTpVPOCS 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 multiparameter algorithms. In general, the best approach for determining the optimal values of the abovementioned five parameters is to find the global optimal value of peak signaltonoise ratio (PSNR) in the fivedimensional 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 abovementioned way, including the parameter α that was defined as the regularization parameter of the AwTVPOCS method. In the synchrotron ILPCCT experiment, the contrasttonoise ratio (CNR) was used as the evaluation measure. The parameters for the simulation and synchrotron ILPCCT 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 (ZellerPlumhoff 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_{1} and c_{2} are two constants to stabilize the division 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.
(iii) RE is defined as follows,
(iv) CNR is defined as follows,
where y_{1} and y_{2} 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_{1} and y_{2}, respectively, and are the variances of y_{1} and y_{2}, 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 AwaTpVPOCS 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 straightline array with 724 bins, and the size of the reconstructed images is 512 × 512 pixels. The 60 uniformly distributed projections over 180° (fullscan range) were used to simulate fewview projections. In addition, the 60 uniformly distributed noisefree and noisy projections scanned from 30° to 120° were used to simulate limitedangle projections and lowdose noisy limitedangle projections, respectively. The lowdose noise model added into limitedangle 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 AwTVPOCS algorithms were used for comparison with the AwaTpVPOCS 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) i54460 CPU at 3.2 GHz and 16 GB RAM.

3.2. Experimental results
In this section, the simulation experiments were performed under three conditions of fewview projections, limitedangle projections and lowdose noisy limitedangle 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 lowdose 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 AwTVPOCS can suppress most of the streak artefacts; however, there are still some artefacts in Figs. 3(g) and 3(k), which were introduced by limitedangle projections or lowdose noise. AwaTpVPOCS achieves the best visual effect in Figs. 3(d,h,i), which show that AwaTpVPOCS eliminates all the streak artefacts, and also eables to better preserve the edge details and suppresses artefacts from limitedangle projections and lowdose noise effectively. The regions of interest (ROIs) were selected and enlarged to improve the visual effect for better comparison, as shown in Fig. 4.
3.3. Assessments
To assess the accuracy of the four reconstruction approaches under fewview, limitedangle and lowdose noisy limitedangle 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 AwaTpVPOCS 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 AwaTpVPOCS approach is obviously the best.

To evaluate the convergence performance of the SART, AwTVPOCS and AwaTpVPOCS approaches under the fewview, limitedangle and lowdose noisy limitedangle conditions, the PSNRbased and REbased convergence curves of the abovementioned approaches are presented, as shown in Figs. 6(a)–6(c). It can be observed that the SART, AwTVPOCS and AwaTpVPOCS algorithms converged before the iterations reach 50 under the fewview condition, but converged before the iterations reach 300 under other two conditions and, obviously, the convergence performance of the AwaTpVPOCS approach is the best.
4. Real experiment on synchrotron ILPCCT 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 ILPCCT image are presented in Fig. 7. In this experiment, the iodine was used as the contrast agent for the contrastenhanced imaging of fibrous tissue in the cirrhosis sample, and its ILPCCT 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 chargecoupled 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 (930view 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 darkcurrent images (the offset signal recorded when no photons hit the detector) were used to remove the dark signal in projections (darkfield correction), while 20 flatfield images (the images with no sample in the beam) were used to correct the pixeltopixel nonuniformity in projections caused by beam's intensity inhomogeneity or nonuniform detector response (flatfield correction) (Chen et al., 2012). After phase retrieval using the PADBA method, the 186view projections were uniformly chosen from the full projection dataset and the 414view 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 fewview CT reconstruction and limitedangle CT reconstruction, respectively. In this experiment, all parameters were chosen for optimal performance, as shown in Table 3.

4.2. Experimental results
Fig. 8 shows reconstructed images of the human cirrhosis sample using the FBP, SART, AwTVPOCS and AwaTpVPOCS algorithms. Figs. 8(a)–8(d) are a reconstructed slice of human cirrhosis sample from 186view projections under fewview conditions using the FBP, SART, AwTVPOCS and AwaTpVPOCS algorithms, respectively. Figs. 8(i)–8(l) are reconstructed slices of human cirrhosis sample with 414view projections under limitedangle condition using the FBP, SART, AwTVPOCS and AwaTpVPOCS 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 fewview CT reconstruction and the limitedangle 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 AwTVPOCS has been improved significantly in comparison with that from the FBP and SART algorithms, and this suggests that the AwTVPOCS algorithms have good capability when it comes to the fewview CT reconstruction and the limitedangle CT reconstruction problems. However, there are a few residual artefacts in the abovementioned results. Figs. 8(d) and 8(l) depict the reconstructed results using the AwaTpVPOCS algorithm, and it can be easily observed that the results of AwaTpVPOCS 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.
4.3. Result analysis
To evaluate the accuracy of the CT images (Fig. 8) reconstructed using the FBP, SART, AwTVPOCS and AwaTpVPOCS algorithms under fewview and limitedangle 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 abovementioned 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 AwTVPOCS and AwaTpVPOCS are close to the reference image, indicating that AwTVPOCS and AwaTpVPOCS have high accuracy in practical CT reconstruction. However, when compared with the profile of the CT images reconstructed using AwTVPOCS, the profile of the CT images reconstructed using AwaTpVPOCS is smoother in the fibrous tissue region and, considering the homogeneous characteristic of the fibrous tissue region, this demonstrates that the AwaTpVPOCS 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 nonoverlapping 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 AwaTpVPOCS in fewview and limitedangle cases are highest; namely, the contrast between the fibrous tissue and the liver parenchyma in the reconstructed images of AwaTpVPOCS 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 AwaTpVPOCS 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.
To test the convergence performance of the SART, AwTVPOCS and AwaTpVPOCS algorithms in practical CT reconstruction under fewview and limitedangle conditions, the RDbased convergence curves of the abovementioned three methods are drawn, as shown in Figs. 10(a) and 10(b), respectively. It can be clearly observed that the SART, AwTVPOCS and AwaTpVPOCS algorithms converged before the iterations reach 30 under fewview and limitedangle conditions, demonstrating that the solutions to the abovementioned three methods in practical CT reconstruction under fewview and limitedangle conditions can converge.
5. Discussion and conclusion
In this study, an AwaTpVPOCS algorithm was proposed for accurate CT reconstruction based on the fewview and limitedangle projections. In the experiment, the LC phantom and the synchrotron ILPCCT data of the human cirrhosis sample were used to evaluate the accuracy and the feasibility of this algorithm, and the FBP, SART and AwTVPOCS algorithms were utilized for comparison. The results confirmed that the AwaTpVPOCS algorithm was an effective approach to decrease the radiation dose for ILPCCT, which had the ability to significantly reduce the artefacts introduced by the fewview and limitedangle sampling, and enabled the edge details to be preserved and the lowdose noise to be suppressed.
Although the AwaTpVPOCS algorithm is a nonconvex optimization method, which can obtain only local optimal solutions in general, it can yield more accurate solutions than the AwTVPOCS 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 ILPCCT, and the accuracy and the feasibility of this algorithm is demonstrated in this research. Taking into account that the small difference in Xray attenuation coefficients between the adjacent soft tissues still constitutes a significant problem for diagnostic interpretation, a limitation for absorptionbased CT in the soft tissues, the phasebased ILPCCT 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, ILPCI experiments carried out on conventional Xray 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 ILPCI 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 ILPCI systems. However, many obstacles still remain for the realization of preclinical or clinical ILPCI 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 ILPCCT while maintaining excellent image quality using newly iterative CT reconstruction algorithms. In this work, a highquality slice of the human cirrhosis sample stained with iodine was reconstructed using the AwaTpVPOCS algorithm from the highly undersampled ILPCCT projections, indicating that the proposed method is a valuable tool for lowdose 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 ILPCCT images and other modal CT images under the assumptions that they are piecewisesmooth. In further research, the graphic processing unit (GPU)based parallel computing techniques will be implemented to increase the computing speed of the AwaTpVPOCS algorithm and, additionally, the adaptivity of the parameters in the AwaTpVPOCS algorithm will be researched. Furthermore, future studies will be carried out to assess whether the AwaTpVPOCS algorithm also applies for in vivo data, along with testing to employ more nonconvex iterative CT reconstruction algorithms to the ILPCI and other PCI methods, such as the specklebased 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 ILPCCT image by means of iodine contrast agent, and it shows good agreement with its pathological section; (ii) a new nonconvex ILPCCT 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 ILPCCT, and promote the wider use of ILPCCT in clinical and biomedical imaging.
Acknowledgements
The authors thank the staff of the BL13W1 beamline of Shanghai Synchrotron Radiation Facility.
Funding information
The following funding is acknowledged: The National Natural Science Foundation of China (grant No. 81671683, 81670545, and 81371549); The Natural Science Foundation of Tianjin City in China (grant No. 16JCYBJC28600); The Foundation of Tianjin university of technology and education (grant No. KJ1201and KJ1736).
References
Andersen, A. H. & Kak, A. C. (1984). Ultrason. Imaging, 6, 81–94. CrossRef CAS PubMed Google Scholar
Bérujon, S. & Ziegler, E. (2016). Phys. Rev. Appl. 5, 044014. Google Scholar
Bérujon, S., Ziegler, E., Cerbino, R. & Peverini, L. (2012). Phys. Rev. Lett. 108, 158102. Web of Science PubMed Google Scholar
Bravin, A. & Coan, P. (2012). Phys. Med. Biol. 57, 2931–2942. Web of Science PubMed Google Scholar
Chartrand, R. (2007). IEEE Signal Process. Lett. 14, 707–710. Web of Science CrossRef Google Scholar
Chartrand, R. (2012). IEEE Trans. Signal Process. 60, 5810–5819. Web of Science CrossRef Google Scholar
Chartrand, R. & Staneva, V. (2010). Inverse Probl. 24, 657–682. Google Scholar
Chen, R.C., Dreossi, D., Mancini, L., Menk, R., Rigon, L., Xiao, T.Q. & Longo, R. (2012). J. Synchrotron Rad. 19, 836–845. Web of Science CrossRef IUCr Journals Google Scholar
Chen, R., Rigon, L. & Longo, R. (2013). Opt. Express, 21, 7384–7399. Web of Science CrossRef CAS PubMed Google Scholar
Chen, Y. & Anastasio, M. A. (2018). Sens. Imaging, 19, 7–26. Web of Science CrossRef PubMed Google Scholar
Chen, Z., Jin, X., Li, L. & Wang, G. (2013). Phys. Med. Biol. 58, 2119–2141. Web of Science CrossRef PubMed Google Scholar
Cloetens, P., Ludwig, W., Baruchel, J., Van Dyck, D., Van Landuyt, J., Guigay, J. P. & Schlenker, M. (1999). Appl. Phys. Lett. 75, 2912–2914. Web of Science CrossRef CAS Google Scholar
Daubechies, I., Defrise, M. & De Mol, C. (2004). Commun. Pure Appl. Math. 57, 1413–1457. Web of Science CrossRef Google Scholar
Gao, H. (2016). Phys. Med. Biol. 61, 7187–7204. Web of Science CrossRef PubMed Google Scholar
Goldstein, T. & Osher, S. (2009). SIAM J. Imaging Sci. 2, 323–343. Web of Science CrossRef Google Scholar
Groso, A., Abela, R. & Stampanoni, M. (2006). Opt. Express, 14, 8103–8110. Web of Science CrossRef PubMed CAS Google Scholar
Gureyev, T. E., Davis, T. J., Pogany, A., Mayo, S. C. & Wilkins, S. W. (2004). Appl. Opt. 43, 2418–2430. Web of Science CrossRef PubMed Google Scholar
Hansen, P. C. & SaxildHansen, M. (2012). J. Comput. Appl. Math. 236, 2167–2178. Web of Science CrossRef Google Scholar
Horng, A., Brun, E., Mittone, A., Gasilov, S., Weber, L., Geith, T., AdamNeumair, S., Auweter, S., Bravin, A., Reiser, M. & Coan, P. (2014). Invest. Radiol. 49, 627–634. Web of Science CrossRef PubMed Google Scholar
Iyer, J., Zhu, N., Gasilov, S., Ladak, H., Agrawal, S. & Stankovic, K. (2018). Biomed. Opt. Expr. 9, 3757–3767. Web of Science CrossRef CAS Google Scholar
Jerri, A. (1977). Proc. IEEE, 65, 1565–1596. CrossRef Web of Science Google Scholar
Ji, D., Qu, G., Hu, C., Liu, B., Jian, J. & Guo, X. (2017). Chin. Phys. B, 26, 93–100. Google Scholar
Labriet, H., Nemoz, C., Renier, M., Berkvens, P., Brochard, T., Cassagne, R., Elleaume, H., Estève, F., Verry, C., Balosso, J., Adam, J. F. & Brun, E. (2018). Sci. Rep. 8, 12491. Web of Science CrossRef PubMed Google Scholar
Larsson, D. H., Vågberg, W., Yaroshenko, A., Yildirim, A. O. & Hertz, H. M. (2016). Sci. Rep. 6, 39074. Web of Science CrossRef PubMed Google Scholar
Liu, Y., Ma, J., Fan, Y. & Liang, Z. (2012). Phys. Med. Biol. 57, 7923–7956. Web of Science CrossRef PubMed Google Scholar
Lohvithee, M., Biguri, A. & Soleimani, M. (2017). Phys. Med. Biol. 62, 9295–9321. Web of Science CrossRef PubMed Google Scholar
Melli, S. A., Wahid, K. A., Babyn, P., Cooper, D. M. L. & Hasan, A. M. (2018). Comput. Med. Imaging Graph. 69, 69–81. Web of Science CrossRef PubMed Google Scholar
Miao, C. & Yu, H. (2015). IEEE Trans. Image Process. 24, 5455–5468. Web of Science CrossRef PubMed Google Scholar
Mohammadi, S., Larsson, E., Alves, F., Dal Monego, S., Biffi, S., Garrovo, C., Lorenzon, A., Tromba, G. & Dullin, C. (2014). J. Synchrotron Rad. 21, 784–789. Web of Science CrossRef IUCr Journals Google Scholar
Nugent, K. A., Gureyev, T. E., Cookson, D. J., Paganin, D. & Barnea, Z. (1996). Phys. Rev. Lett. 77, 2961–2964. CrossRef PubMed CAS Web of Science Google Scholar
Pan, J., Liu, R., Su, Z. & Gu, X. (2013). Signal Process. Image Commun. 28, 1156–1170. Web of Science CrossRef Google Scholar
Peña Fernández, M., Cipiccia, S., Dall'Ara, E., Bodey, A., Parwani, R., Pani, M., Blunn, G., Barber, A. & Tozzi, G. (2018). J. Mech. Behav. Biomed. Mater. 88, 109–119. Web of Science PubMed Google Scholar
Preissner, M., Murrie, R. P., Pinar, I., Werdiger, F., Carnibella, R. P., RZosky, G., Fouras, A. & Dubsky, S. (2018). Phys. Med. Biol. 63, 1–9. Web of Science CrossRef Google Scholar
Sethasine, S., Jain, D., Groszmann, R. J. & GarciaTsao, G. (2012). Hepatology, 55, 1146–1153. Web of Science CrossRef PubMed Google Scholar
Shu, X. & Ahuja, N. (2010). Computer VisionECCV 2010, edited by K. Daniilidis, P. Maragos & N. Paragios, Vol. 6316 of Lecture Notes in Computer Science, pp. 393–404. Berlin: Springer. Google Scholar
Siddon, R. L. (1985). Med. Phys. 12, 252–255. CrossRef CAS PubMed Web of Science Google Scholar
Sidky, E. Y., Chartrand, R., Boone, J. M. & Pan, X. (2014). IEEE J. Transl. Eng. Heal. Med. 2, 1–18. CrossRef Google Scholar
Sidky, E. Y., Chartrand, R. & Pan, X. (2007). IEEE Nucl. Sci. Symp. Conf. Rec. 5, 3526–3530. Google Scholar
Sidky, E. Y., Kao, C. M. & Pan, X. (2006). J. Xray Sci. Technol. 14, 119–139. Google Scholar
Wang, L., Li, X., Wu, M., Zhang, L. & Luo, S. (2013). Biomed. Eng. Online, 12, 1–13. Web of Science PubMed Google Scholar
Wang, Q., Zhang, X., Wu, Y., Tang, L. & Zha, Z. (2017). IEEE Signal Process. Lett. 24, 1686–1690. CrossRef Google Scholar
Wang, Y., Yang, J., Yin, W. & Zhang, Y. (2008). J. Sci. Comput. 1, 248–272. Google Scholar
Wang, Z., Bovik, A. C., Sheikh, H. R. & Simoncelli, E. P. (2004). IEEE Trans. Image. Process. 13, 600–612. Web of Science CrossRef PubMed Google Scholar
Woodworth, J. & Chartrand, R. (2016). Inverse Probl. 32, 4–28. Web of Science CrossRef Google Scholar
Wu, L., Chen, Y., Jin, J., Du, H. & Qiu, B. (2017). J. Electron. Imaging, 26, 053003. Web of Science CrossRef Google Scholar
ZellerPlumhoff, B., Mead, J. L., Tan, D., Roose, T., Clough, G. F., Boardman, R. P. & Schneider, P. (2017). Opt. Express, 25, 33451–33468. CAS Google Scholar
Zhang, L., Zeng, L., Wang, C. & Guo, Y. (2018). J. Inverse IllPosed Probl. 26, 799–820. Web of Science CrossRef Google Scholar
Zhao, Y., Brun, E., Coan, P., Huang, Z., Sztrókay, A., Diemoz, P. C., Liebhardt, S., Mittone, A., Gasilov, S., Miao, J. & Bravin, A. (2012). Proc. Natl Acad. Sci. USA, 109, 18290–18294. Web of Science CrossRef CAS PubMed Google Scholar
Zhao, Y., Sun, M., Ji, D., Cong, C., Lv, W., Zhao, Q., Qin, L., Jian, J., Chen, X. & Hu, C. (2018). J. Synchrotron Rad. 25, 1450–1459. Web of Science CrossRef IUCr Journals Google Scholar
Zuo, W., Meng, D., Zhang, L., Feng, X. & Zhang, D. (2013). Proceedings of the 2013 IEEE International Conference on Computer Vision, 1–8 December 2013, Sydney, NSW, Australia, pp. 217–224. Google Scholar
Zysk, A. M., Garson, A. B., Xu, Q., Brey, E. M., Zhou, W., Brankov, J. G., Wernick, M. N., Kuszak, J. R. & Anastasio, M. A. (2012). Biomed. Opt. Express, 3, 1924–1932. Web of Science CrossRef PubMed Google Scholar
© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.