## Abstract

In-line X-ray phase-contrast computed tomography (IL-PCCT) can produce high-contrast and high-resolution images of biological samples, and it has a great advantage with regard to imaging the microstructures and morphologies of fibrous biological tissues (FBTs). Filtered back projection (FBP) is widely used in ILPCCT. However, it requires long scanning times and high radiation doses to produce high-quality CT images, and this restricts its applicability in biomedical and preclinical studies on FBTs. To solve this problem, a novel IL-PCCT reconstruction algorithm is proposed to decrease the radiation dose by reducing the number of projections and reconstruct high-quality CT images of FBTs. The proposed algorithm incorporates the FBP method into the iterative reconstruction framework. Considering the area types and anisotropic edge properties of FBTs, a discriminant adaptive-weighted total variation model is introduced to optimize the intermediate reconstructed images. A fibrous phantom simulation and real experiment were performed to assess the performance of the proposed algorithm. Simulation and experimental results demonstrated that the proposed algorithm is an effective IL-PCCT reconstruction method for FBTs with incomplete projection data, and it has a great ability to suppress artifacts and preserve the edges of fibrous structures.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Fibrous structures are major components of fibrous biological tissues (FBTs), such as fibrous liver tissue and breast tissue, and they are closely related to the initiations and progressions of some diseases. For example, collagen fiber is one of the most important pathological features in the development of liver fibrosis, cirrhosis and other diseases, and revealing its structural changes has great significance for the diagnosis and treatment of these diseases [1]. Analyzing the fibrous structure and morphology of FBT is of great value for pathological investigations in biomedical and preclinical applications. A nondestructive imaging technique is important for visualization and quantitative studies on FBTs. X-ray imaging is a widely used technique for structure visualization of biological samples. However, due to the limited spatial resolution and contrast [2,3], the conventional absorption-based X-ray imaging technique has a poor ability to distinguish weakly absorbing tissues (e.g., biological soft tissues). X-ray phase-contrast imaging (PCI) generates image contrast from X-ray refraction and scattering, and it has the potential to overcome the drawback of conventional absorption-based X-ray imaging with regard to observing biological soft tissues [4]. Generally, PCI has approximately three orders of magnitude higher sensitivity than conventional absorption-based X-ray imaging for elements with low atomic numbers, such as carbon, hydrogen and oxygen [5]. Among PCI techniques, in-line PCI (IL-PCI) represents in practice the simplest way to implement since no extra optics is required. As a combination of IL-PCI and computed tomography (CT), IL X-ray phase-contrast CT (IL-PCCT) can achieve high-resolution three-dimensional visualizations of biological soft tissues [6]. Owing to the remarkable advancements in synchrotron radiation (SR) X-ray sources and detector techniques, the spatial resolution of IL-PCCT has reached a submicron level [7], so it can effectively reveal the complex microstructures inside biological samples. Therefore, it has a great advantage for imaging and analyzing the microstructures of FBTs.

However, IL-PCCT suffers from the problems of requiring a long scanning time and a high radiation dose. These not only cause health risks to subjects of *in vivo* imaging, but also lead to damage to the structures and biological properties of *ex vivo* biological samples [8], which may impact the quantitative analysis and subsequent use. It is hence necessary to use a shorter scanning time and a lower radiation dose to reconstruct numerically accurate CT images. To reduce the CT scanning time and radiation dose of a given sample, a feasible solution is to decrease the number of projections, e.g., by using few-view and limited-angle projections [9]. In general, streak artifacts are yielded when the projection data are incomplete using traditional analytic reconstruction (AR) algorithms, such as filtered back projection (FBP) [10,11]. Iterative reconstruction (IR) algorithms can effectively solve this problem [12–14]. They can use optimization methods to suppress artifacts and noise and produce high-quality reconstructed images. Among these optimization methods, total variation (TV) has been widely used. Considering the piecewise constant or sparse source distribution, Sidky *et al.* introduced the concept of TV minimization and first proposed the TV-projection onto convex sets (POCS) algorithm, which can reconstruct high-quality CT images under few-view and limited-angle conditions [15,16]. However, the assumption of TV that the signal is piecewise smooth makes TV suffer from the oversmoothing in image edges. As a result, many improved TV methods have been proposed. Tian *et al.* proposed a TV-based edge preserving (EPTV) model [17], which can preserve edges by introducing different weights in the TV term from edges and constant areas of the image. Different from the EPTV model, Liu *et al.* considered the anisotropic edge property of an image and proposed an adaptive-weighted TV (AwTV) model [18]. In AwTV, the associated weights are expressed as an exponential function that can be adaptively adjusted along with the local image-intensity gradient to preserve edge details. However, this model does not take the area differences in the image into account and only uses a weighting function to process and optimize the areas without discriminant ability.

Generally, FBTs can be divided into two area types: flat areas and complex areas including fibrous structures. It is important to achieve the synergistic effect of edge preservation and non-edge smoothing by processing the two types of areas separately. Among the techniques for doing this, the combination of area discrimination and image filtering is a valuable research direction. Following that the application of that combination, the areas are discriminated, and the appropriate weighting functions are applied to the flat area (with small intensity gradients) and complex area (with large intensity gradients). During this process, the diffusion intensity of the flat area can be kept at a high level to promote smoothing and suppress artifacts, and the diffusion intensity of the complex area can be kept at a low level to effectively preserve the edge details. The P-M model is an anisotropic diffusion filtering algorithm proposed by Perona and Malik [19], where the two types of diffusion coefficients provided can produce similar filtering effects, but the diffusion mechanisms are different. Among them, one type preferentially diffuses high-contrast boundaries, i.e., complex areas, and the other type preferentially diffuses non-boundary areas, i.e., flat areas.

In this study, a novel discriminant AwTV-iterative FBP (DAwTV-IFBP) algorithm is proposed for FBT reconstruction from an incomplete projection in IL-PCCT. IR algorithms have good reconstruction performances for incomplete projection data at the cost of increased computational costs [20]. Therefore, we incorporate FBP into the IR framework. The proposed algorithm not only has the advantage of a low computational cost (from FBP) but also preserves the advantages of IR algorithms, such as immunity to few projections and the addition of prior information [21]. At the same time, considering the area types and anisotropic edge property of FBT, a DAwTV regularization method is introduced to optimize the intermediate reconstructed image. It establishes an area discrimination method based on the TV model. First, based on the block strategy, the information entropy parameter, which can reflect the average information content in the image, is used to discriminate between the flat and complex areas of the image [22], and then the two types of diffusion functions in the P-M model are used as the weighting functions of TV to achieve targeted adaptive weighting for the image-intensity gradients of the two areas to effectively suppress artifacts and preserve the edge details of the complex areas. Finally, a fibrous phantom simulation and an IL-PCCT experiment on a human cirrhosis sample are performed to assess the performance of the DAwTV-IFBP algorithm.

## 2. Methods

#### 2.1 IL-PCI and phase retrieval

In IL-PCI, a density difference within the sample causes a change in the X-ray phase shift when the spatially coherent X-ray beams pass through a sample. As the X-rays carrying the phase shift information propagate in the near-field regime, the phase shift is converted into light intensity information and measured by the detector due to the Fresnel diffraction effect [23]. The schematic experimental setup of IL-PCI based on SR is shown in Fig. 1. Because of its advantages of using simple devices and having an easy implementation procedure, IL-PCI possesses high potential in the applications of biomedicine and other fields. When parallel monochromatic X-rays illuminate a sample perpendicularly, the complex refractive index of the sample can usually be described as follows:

where (*x*,

*y*) denotes the coordinates of the projection plane,

*z*denotes the propagation axis, and $\delta (x,y,z)$ and $\beta (x,y,z)$ denote the phase information and the absorption information, respectively.

For IL-PCI, the measured light intensity information contains absorption information and the Laplacian of the phase information. Before CT reconstruction is performed, it is necessary to use the phase retrieval method to extract phase-shift information from the projection images. The phase-attenuation duality Paganin (PAD-PA) method is one of the most widely used phase retrieval methods in the field of biomedical imaging, and it can accurately extract phase-shift information from biological samples. This study uses PAD-PA to solve the problem of phase retrieval. The weakly absorbing sample meets the PAD property [24], i.e. $\delta (x,y,z) = \gamma \beta (x,y,z)$, where $\gamma $ is a constant, and the PAD-PA method can be described by the following formula:

*θ*, $I_\theta ^D(x,y)$ denotes the intensity function measured by the detector at a distance

*D*,

*λ*denotes the X-ray wavelength, $(\tau ,\textrm{ }\upsilon )$ are the coordinates of (

*x*,

*y*) in the Fourier domain, and ${\cal F}$ and ${{\cal F}^{ - 1}}$ are the 2D Fourier transform and its inverse operators, respectively.

#### 2.2 DAwTV-IFBP algorithm for IL-PCCT

Synchrotron radiation is a monochromatic source, and thus, the general model of CT imaging can be approximately represented as the following discrete linear system after phase retrieval:

where*A*denotes the system matrix, $\varphi$ denotes the measured projection data, and $\delta (x,z)$ is the phase CT image to be reconstructed. In the following, we use the abbreviation $\delta$ to represent $\delta (x,z)$.

This study aims to reconstruct $\delta$ from the projection data $\varphi$ and the system matrix *A*. In practice, Eq. (3) is known as an ill-posed problem due to the insufficient projection data resulting from the few-view or limited-angle conditions. To solve the discrete linear system represented in Eq. (3), an optimization-based strategy is often adopted. In this study, Eq. (3) is transformed into a constrained optimization problem, and the objective function of the phase CT reconstruction task can be formulated as follows:

*R*denotes the regularization term, and ${\rho ^{\ast }}$ is a penalty parameter.

For FBTs, the edge details are important information. In this work, a DAwTV regularization method is proposed to alleviate the problem of oversmoothed edges, which are present in the traditional TV regularization method. The DAwTV of the reconstructed image $\delta$, i.e. ${||\delta ||_{DAwTV}}$, is defined as:

Considering that a FBT has both flat areas and complex areas including fibrous structures, this study intends to effectively smooth the flat areas and preserve the edge details of the complex areas. Therefore, it is particularly important to discriminate between the area types and select appropriate weighting functions to address the two types of areas.

In terms of area discrimination, this study is based on a block strategy, dividing an image of size $M \times M$ into *L* block areas of size $N \times N$, where $\sqrt L \times N = M$. Two conditions must be ensured between the block areas: (1) First, they are independent and do not overlap; (2) The sum of the block areas can cover the entire original image. This block strategy reflects the independence and integrity of the block areas. After the image is divided into blocks, the information entropy of each block area is calculated, and the calculation formula is expressed as follows:

*H*represents the information entropy of the block area;

*k*denotes the number of gray levels, and $P(s)$ denotes the probability of gray level

*s*being present in the block area. The information entropy commonly reflects the texture transformation of the block area [25]. The information entropy becomes high when the texture is rich, and an area with high information entropy can be considered a complex area; otherwise, it can be considered a flat area. Here, we use

*u*as a criterion for area discrimination, and it is defined as follows: where ${H_l}$ denotes the information entropy of the $l$th block area, ${H_{mean}}$ is the average value of the information entropy from

*k*block areas, and

*ε*is an adjustable parameter that is used to adjust the size of the criterion in the area discrimination, its value range is 0.8-1. Based on Eq. (7), the block area can be considered a complex area if $H \ge u$ and a flat area if $H\lt u$.

After discriminating the area, targeted optimization is performed for each of the two areas. The P-M model contains two diffusion functions, i.e., ${C_1}(\textrm{g}) = \exp [ - {(g/\sigma )^2}]$ and ${C_2}(g) = \frac{1}{{1 + {{(g/\sigma )}^2}}}$, where *g* is the magnitude of the gradient and $\sigma$ is a scale parameter. Among them, ${C_1}$ preferentially diffuses complex areas, and ${C_2}$ preferentially diffuses flat areas. They are both used as the weighting functions of the DAwTV model, and they are applied to the image-intensity gradient of the two areas. Here, an adaptive selection method for the weighting function based on the block area is introduced. If the block area is a complex area, the ${C_1}$ function is used as the weighting function to reduce the diffusion intensity and preserve the edge details; if the block area is a flat area, the ${C_2}$ function is used as the weighting function to enhance the diffusion strength and promote smoothing. The formula for adaptively selecting the weighting function is as follows:

Using DAwTV as the regularization term, Eq. (4) can be converted to the following minimization problem:

*t*denotes the number of iterations. During each iteration, a trial solution constrained by the data fidelity term is first obtained in Eq. (11). Then, Eq. (12) adopts the DAwTV term to effectively remove the noise and artifacts in the trial solution and improve image quality.

### 2.2.1 Iterative reconstruction strategy based on the FBP algorithm

For the subproblem in Eq. (11), ${A^T}(A{\delta ^t} - \varphi )$ can be seen as the error image reconstructed by the projection residual. By superimposing the error image on the current reconstructed image, the projection residual of the reconstruction result can be reduced, making the result closer to the real image. However, the computational cost of ${A^T}$ is high, which not only requires more computer memory but also leads to a longer reconstruction time. Considering that the FBP algorithm is simple and efficient and that the iterative algorithm has the advantage of high reconstruction accuracy, the IFBP algorithm is utilized in this study. FBP is regarded as a back-projection operator, and it is used to replace ${A^T}$. We define *Q* as a FBP operator, and Eq. (11) can be expressed as:

### 2.2.2 Optimization process based on DAwTV

The subproblem in Eq. (12) represents the DAwTV regularization problem. To solve this problem, the gradient descent method is implemented [17,27]. The formula is as follows:

where $\rho$ denotes the ratio of the parameters ${\rho ^{\ast }}$ and $\mu$. The functional variation of DAwTV, i.e., $\frac{{\partial {{||\delta ||}_{DAwTV}}}}{{\partial {\delta _{i,j}}}}$, is calculated as follows:### 2.2.3 Overall workflow

Two steps are involved in the implementation of the DAwTV-IFBP algorithm. In the first step, an initially estimated image is updated iteratively to fulfill the data constraints using the IFBP strategy. An intermediate image is yielded by DAwTV regularization for further updating in the second step. Before the IFBP step, FBP is performed on the actual undersampled projection data and inpainted projection data obtained by interpolation between the undersampled angles [28]. By fusing the two reconstructed images, a good initial image is obtained. The IFBP step is used to acquire CT reconstruction images with artifacts. Following the IFBP step, the DAwTV step is performed to effectively reduce the artifacts and improve image quality. Subsequently, the updated image from the DAwTV step is reprojected to produce the forward projection data, and the projection residual is obtained by subtracting the forward projection data and the inpainted projection data. Then, FBP is used to reconstruct the error image of the projection residual, and the error image is superimposed with the previous updated image. After repeated iterations, the projection residual of the reconstructed image is gradually reduced, and the reconstruction accuracy of the image is improved. To explain the DAwTV-IFBP algorithm more clearly, the corresponding pseudocode is given, as shown in the Appendix in Sec. 6.1.

#### 2.3 Parameter selection

In the DAwTV-IFBP algorithm, there are 5 main parameters: $\alpha$, $\rho$, ${\sigma _1}$, ${\sigma _2}$ and *N*. These parameters have significant influences on the reconstruction result, so the optimal selection of their values is important. In this study, a greedy strategy is used to determine the optimal parameter values one by one [29]. In the simulation experiment, the peak signal-to-noise ratio (PSNR) is used as an evaluation metric. As shown in Figs. 9 (a, b), to determine the optimal value of the parameter *N*, different *N* values are selected while keeping other parameters unchanged to obtain different reconstructed images. Then, the corresponding PSNR values are calculated, and the optimal *N* value can be determined by the highest PSNR value. The other parameters are also determined in this manner, including the scale parameter $\sigma$ in the AwTV-POCS and AwTV-IFBP algorithms. In the real IL-PCCT experiment, the contrast-to-noise ratio (CNR) is used in the same way to determine the optimal parameter set.

#### 2.4 Quantitative assessment

Five metrics, the PSNR, structural similarity index (SSIM), edge preservation index (EPI), CNR, signal-to-noise ratio (SNR) and relative error (RE), are used to evaluate the quality of reconstructed images.

The PSNR is a traditional measure of image quality, and a larger PSNR value represents better quality. It is defined as follows:

*MSE*is the mean square error function,

*x*denotes the reference image,

*y*denotes the reconstructed image, and $M \times M$ is the size of

*x*and

*y*; ${x_{i,j}}$ and ${y_{i,j}}$ represent the pixel intensities of

*x*and

*y*in some pixel (

*i*,

*j*), respectively, and

*Peak*represents the largest pixel intensity in the normalized image (it is 255 in our study).

The SSIM is usually used to assess the similarity between the reconstructed and reference images, and a larger SSIM value represents higher similarity [30]. It is defined as follows:

*x*and

*y*, respectively; $\sigma _x^2$ and $\sigma _y^2$ denote the variances of

*x*and

*y*, respectively; ${\sigma _{xy}}$ is the covariance of

*x*and

*y*; and ${c_1}$ and ${c_2}$ are two constants used to stabilize a division operation with a weak denominator.

The EPI is used to assess the edge preservation ability of the model with respect to reconstructed images [31], and a larger EPI value represents a stronger preservation ability. It is defined as follows:

*x*and

*y*, and they are approximated by the 3×3-pixel standard of the Laplace operator; $\overline {\Delta x}$ and $\overline {\Delta y}$ denote the mean pixel values of $\Delta x$ and $\Delta y$, respectively; and $\Gamma (x,y) = \sum\limits_{i,j} {x(i,j) \cdot y} (i,j)$.

The CNR is adopted to evaluate the contrast between adjacent tissues, and a larger CNR value represents higher contrast [32]. It is defined as follows:

The SNR is a traditional measure of image quality, and a larger SNR value indicates better image quality [32]. It is defined as follows:

The RE can be used to evaluate the relative error between the reconstructed and reference images, and a smaller RE value indicates higher accuracy. It is defined as follows:## 3. Simulation

#### 3.1 Simulation data

To evaluate the reconstruction performance of the DAwTV-IFBP algorithm on FBTs, a fibrous phantom with 512×512 pixels was adopted for simulation purposes [33], as shown in Fig. 2(a). The fibrous phantom was a grayscale image with a gray value ranging from 0 to 255, and it was used to simulate the phase information distribution of the FBT. The detector was modeled as a straight-line array of 840 detector bins. 60 uniformly distributed projections at angular ranges from 0° to 180° were selected to simulate few-view projection, and 107 uniformly distributed projections at angular ranges from 10° to 170° were selected to simulate the limited-angle projection. To further evaluate the proposed DAwTV-IFBP algorithm, we compared DAwTV-IFBP with FBP, SART, IFBP, AwTV-POCS and AwTV-IFBP, and all parameters were obtained by optimal selection. According to the convergence curves (Fig. 8), the iterative algorithms adopted the total numbers of required iterations as the stopping criteria. All experiments were carried out in MATLAB 2017b on a PC equipped with an Intel Core i7-6700 CPU at 3.4 GHz and 8 GB RAM.

#### 3.2 Simulation results

### 3.2.1 Few-view reconstruction results

For the few-view condition, the parameters were set as follows: (1) DAwTV-IFBP: $\alpha = 0.1$, $\rho = 0.4$, ${\sigma _1} = 50$, ${\sigma _2} = 10$ and $N = 16$; (2) IFBP: $\alpha = 0.1$; (3) AwTV-POCS: $\rho = 0.5$ and $\sigma = 50$; and (4) AwTV-IFBP: $\alpha = 0.1$, $\rho = 0.4$ and $\sigma = 50$.

The few-view reconstruction results for the fibrous phantom and the corresponding absolute difference images (blue: low error, red: high error) are shown in Fig. 3. Due to the few-view projection data, the result of FBP exhibits severe streak artifacts (Fig. 3(a)). SART and IFBP can only remove some of the streak artifacts (Fig. 3(c) and Fig. 3(e)). Compared with the previous three algorithms, AwTV-POCS, AwTV-IFBP and DAwTV-IFBP can suppress artifacts more effectively (Figs. 3(g, i, k)). However, fewer errors are observed at the reconstruction results of DAwTV-IFBP according to the absolute difference images (Fig. 3(l)). To better compare and analyze the results, the blue and green rectangular regions of interest (ROIs) in Fig. 3(a) are enlarged in Fig. 4. DAwTV-IFBP suppresses most of the artifacts; moreover, the edges and details are better preserved, as shown by the red arrows and red circle (Fig. 4(k) and Fig. 4(w)).

### 3.2.2 Limited-angle reconstruction results

For the limited-angle condition, the parameters were set as follows: (1) DAwTV-IFBP: $\alpha = 0.1$, $\rho = 0.4$, ${\sigma _1} = 100$, ${\sigma _2} = 50$ and $N = 8$; (2) IFBP: $\alpha = 0.1$; (3) AwTV-POCS: $\rho = 0.4$ and $\sigma = 100$; and (4) AwTV-IFBP: $\alpha = 0.1$, $\rho = 0.4$ and $\sigma = 100$.

In Fig. 5, the limited-angle reconstruction results for the fibrous phantom and the corresponding absolute difference images (blue: low error, red: high error) are given. Figures 5 (a, c, e) indicate that the reconstruction results of FBP, SART and IFBP have poor visual effects, where the detailed structures and the edges of the fibrous tissue are severely affected by artifacts. AwTV-POCS can suppress most artifacts, but there are still some artifacts introduced by limited-angle projections (Fig. 5(g)). AwTV-IFBP and DAwTV-IFBP have better abilities to suppress artifacts than the previous four algorithms (Fig. 5(i) and Fig. 5(k)). The ROIs are enlarged in Fig. 6. It can be observed that AwTV-IFBP can eliminate the artifacts effectively, but it blurs the edges of the tissues (Fig. 6(i) and Fig. 6(u)). DAwTV-IFBP suppresses more artifacts while effectively maintaining the edges, as indicated by the red arrow, and the details shown in the red circle are more complete (Fig. 6(k) and Fig. 6(w)). The absolute difference images in Fig. 5 and Fig. 6 show the same results. It is clearly to see that DAwTV-IFBP suppresses most of artifacts and preserves more edges and details.

#### 3.3 Assessments

Figure 7 shows the horizontal profiles of the same position (the position is marked with the red line in Fig. 2(a)) in the reconstructed images. Here, the horizontal profile from the ground truth was given as a reference. The results show that the profiles obtained by the DAwTV-IFBP algorithm are closest to the reference profiles under the two conditions, with better smoothness in the flat area and less fluctuations near the edge of the fibrous structure. The quantitative results are given in Table 1. It can be clearly seen that the DAwTV-IFBP has the highest PSNR, SSIM and EPI values. Furthermore, compared with SART and AWTV-POCS, the DAwTV-IFBP algorithm has a higher reconstruction speed.

As shown in Fig. 8, the RE curves versus the number of iterations were plotted to assess the reconstruction convergence of IFBP, SART, AwTV-POCS, AwTV-IFBP and DAwTV-IFBP. It can be observed that the algorithms converged within 55 iterations under the few-view condition and converged within 300 iterations under the limited-angle condition. Thus, in the simulation, the stopping criteria of the five iterative algorithms were set as 55 iterations under the few-view condition and 300 iterations under the limited-angle condition.

#### 3.4 Selection of parameter N

To analyze the influence of the block area side length *N* on the reconstruction results of DAwTV-IFBP, experiments were performed. Since the size of the fibrous phantom is 512×512 pixels, i.e., $M = 512$, under the few-view and limited-angle conditions, $N$=2, 4, 8, 16 and 32 were selected to obtain the reconstructed images, and the corresponding PSNR and SSIM values were calculated, as shown in Fig. 9. It can be seen that the PSNR reaches a peak at $N = 16$, and it slowly declines as *N* increases further under the few-view condition (Fig. 9 (a)). The same situation also exists under the limited-angle condition, the PSNR reaches its peak value when $N = 8$ and decreases gradually as *N* increases further (Fig. 9 (b)). However, regardless of how the value of *N* selected, the SSIM stays at stable levels under both the few-view and limited-angle conditions (Fig. 9 (c) and Fig. 9 (d)). The results indicate that the selection of *N* only has a great influence on the reconstructed image quality and does not destroy the structures and accuracy values of the reconstructed results. Based on the results, the parameter *N* of the DAwTV-IFBP algorithm was set as 16 under the few-view condition and 8 under the limited-angle condition in the simulation.

## 4. Real experiment on IL-PCCT data

#### 4.1 Data acquisition

In this study, IL-PCCT reconstruction was performed using an ex vivo human viral cirrhosis sample provided by the Beijing Friendship Hospital affiliated with Capital Medical University. The experiment was approved by the Research Ethics Committee of Beijing Friendship Hospital, Capital Medical University, Beijing, China, and written informed consent was obtained from the patients. In this experiment, iodine was used as the contrast agent to further enhance the imaging contrast of fibrous tissue in cirrhosis sample, and the corresponding IL-PCCT data were acquired at the BL13W beamline of SSRF, Shanghai, China. For IL-PCCT imaging, the X-ray energy was set as 33 keV, and the SDD was set as 0.8 m. A charge-coupled device (CCD) camera (Photonic Science, UK) with a 36 mm (horizontal) × 4.5 mm (vertical) field of view (FOV) was adopted as the imaging detector, with an effective pixel size of 9 µm × 9 µm. 949 projections were collected from 0° to 180° at equidistant angles with an exposure time of 12 ms for each projection, and the total imaging time was approximately 25 minutes. In addition, 10 dark-current images and 20 flat-field images were used for dark-field and flat-field correction in the projections, respectively. In this experiment, the PAD-PA method was used for phase retrieval, and its parameter $\gamma $ was set as 400 [34]. A CT image with a size of 2826× 2826 was reconstructed using the FBP algorithm based on the full projection dataset of the 290th layer of the CT dataset, and a region with a size of 512× 512 was selected as our reference image (as shown in Fig. 10(a)). At the same time, 233 projections were selected from the full projection dataset with equal angle intervals to simulate the few-view CT reconstruction experiment, and 414 projections were selected from the projection dataset at angular ranges from 10° to 170° to simulate the limited-angle CT reconstruction experiment. The above data were reconstructed, and the same areas as those in Fig. 10(a) were chosen as our research objects. Among them, all results were obtained after optimal selection of the algorithm parameters. The number of iterations for all iterative algorithms in this experiment was set to 10, thereby guaranteeing that each algorithm could reach convergence.

#### 4.2 Experimental results

### 4.2.1 Few-view reconstruction results

Under the few-view condition, the parameters were set as follows: (1) DAwTV-IFBP: $\alpha = 0.1$, $\rho = 0.3$, ${\sigma _1} = 50$, ${\sigma _2} = 100$ and $N = 9$; (2) IFBP: $\alpha = 0.3$; (3) AwTV-POCS: $\rho = 0.3$ and $\sigma = 50$; and (4) AwTV-IFBP: $\alpha = 0.1$, $\rho = 0.2$ and $\sigma = 50$.

Figures 11(a, b, c, g, h, i) show the few-view reconstruction results, and Figs. 11(d, e, f, j, k, l) are the enlarged images from the same regions in the reconstructed images. Figures 11(a-f) show that the details and edges of the fibrous tissue are seriously affected by streak artifacts in the FBP, SART and IFBP results, and these make it difficult to identify clinically valuable structures. Although AwTV-POCS removes some artifacts, the result is not satisfactory (Fig. 11(g) and Fig. 11(j)). AwTV-IFBP can remove a large number of streak artifacts, but there is still obvious damage to the edges (Fig. 11(h) and Fig. 11(k)). Compared with other algorithms, DAwTV-IFBP can eliminate most of the artifacts, smooth the liver parenchyma and better preserve the edges of the fibrous tissue (Fig. 11(i) and Fig. 11(l)).

### 4.2.2 Limited-angle reconstruction results

Under the limited-angle condition, the parameters were set as follows: (1) DAwTV-IFBP: $\alpha = 0.2$, $\rho = 0.3$, ${\sigma _1} = 150$, ${\sigma _2} = 100$ and $N = 9$; (2) IFBP: $\alpha = 0.2$; (3) AwTV-POCS: $\rho = 0.3$ and $\sigma = 50$; and (4) AwTV-IFBP: $\alpha = 0.2$, $\rho = 0.2$ and $\sigma = 50$.

Figures 12(a, b, c, g, h, i) show the limited-angle reconstruction results, and Figs. 12(d, e, f, j, k, l) are the enlarged images from the same areas in the reconstructed images. It can be seen from Figs. 12(a-f) that the results of FBP, SART and IFBP are seriously affected by limited-angle artifacts, resulting in poor image quality. Compared with the previous three algorithms, the visual effect of AwTV-POCS is improved, but there are still slight artifacts (Fig. 12(g) and Fig. 12(j)). AwTV-IFBP and DAwTV-IFBP can effectively suppress more artifacts, but DAwTV-IFBP has the better ability to maintain the edges of the fibrous tissue (Figs. 12(h, i, k, l)).

#### 4.3 Results analysis

Figure 13 depicts the profiles of the reconstructed images (the profile positions under the few-view and limited-angle conditions are indicated by the red lines in Fig. 11(a) and Fig. 12(a), respectively). The profiles from the ground truth are given as references, and the red arrows indicate the differences in structure details. It can be observed that the profiles generated from DAwTV-IFBP are closest to the references. In addition, Table 2 shows the CNR, SNR values and reconstruction times of the six reconstruction algorithms under the two conditions. The higher CNR and SNR values confirmed that DAwTV-IFBP can effectively distinguish tissues and reduce the noise and artifacts. At the same time, the reconstruction times of DAwTV-IFBP are greatly reduced compared with those of SART and AwTV-POCS, indicating that DAwTV-IFBP has higher computational efficiency.

## 5. Discussion and conclusion

In this study, the DAwTV-IFBP algorithm was proposed for FBT reconstruction in IL-PCCT from few-view and limited-angle projections. The fibrous phantom and the IL-PCCT data of human cirrhosis sample were utilized to assess the performance of this algorithm, and FBP, SART, IFBP, AwTV-POCS, and AwTV-IFBP were used for comparison. The results demonstrated that the DAwTV-IFBP algorithm is an effective method for IL-PCCT reconstruction of FBTs from incomplete projection data, and it has a great ability to suppress artifacts and preserve the edges of fibrous structures.

In the DAwTV-IFBP algorithm, the information entropy parameter was used to discriminate between different areas, but this may not be particularly accurate in the initial stage of reconstruction. Because both flat areas and complex areas contain many artifacts and noise during this stage, they are considered areas with rich texture information. However, with the gradual improvement of image quality in the later stage of the reconstruction method, the discriminant accuracy of the information entropy becomes higher, so the DAwTV process can optimize the image more finely and obtain more ideal results.

The selection of the parameters has an important influence on the performance of DAwTV-IFBP. This algorithm contains many parameters, and the manual selection of these parameters is a complicated process. Furthermore, the optimal selection for each image may be different, making this problem more complicated. Further research will be conducted using the training method or adaptive parameter selection method to simplify the parameter adjustment process. In addition, compared with SART and AwTV-POCS, DAwTV-IFBP greatly reduces the reconstruction time but still has space for optimization. Among the current methods of reducing computational costs, one feasible method is the ordered subset method, which can be considered a parallelization of computational operations. Another popular method is the graphics processing unit (GPU)-based speedup technique [17,35]. In the future, these powerful methods will be utilized to further optimize the performance of DAwTV-IFBP.

## 6. Appendix

## 6.1 Pseudocode

The pseudocode of DAWTV-IFBP algorithm is shown in algorithm 1.

## Funding

National Natural Science Foundation of China (82071922, 81671683, 81670545); Natural Science Foundation of Tianjin City (16JCYBJC28600).

## Acknowledgements

The authors would like to thank the staffs from beamline (BL13W1) of SSRF, China, for their kindly assistance in our experiments.

## Disclosures

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

## References

**1. **S. Sethasine, D. Jain, R. J. Groszmann, and G. Garcia-Tsao, “Quantitative histological-hemodynamic correlations in cirrhosis,” Hepatology **55**(4), 1146–1153 (2012). [CrossRef]

**2. **T. Takeda, A. Momose, Y. Itai, W. Jin, and K. Hirano, “Phase-contrast imaging with synchrotron X-rays for detecting cancer lesions,” Acad. Radiol. **2**(9), 799–803 (1995). [CrossRef]

**3. **F. Pfeiffer, J. Herzen, M. Willner, M. Chabior, and F. Bamberg, “Grating-based X-ray phase contrast for biomedical imaging applications,” Z. Med. Phys. **23**(3), 176–185 (2013). [CrossRef]

**4. **A. Momose, “X-ray phase imaging reaching clinical uses,” Phys. Med. **79**, 93–102 (2020). [CrossRef]

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

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

**7. **L. Xu, R. Chen, Y. Yang, B. Deng, G. Du, H. Xie, and T. Xiao, “Monochromatic-beam-based dynamic X-ray microtomography based on OSEM-TV algorithm,” J. X-Ray Sci. Technol. **25**(6), 1007–1017 (2017). [CrossRef]

**8. **F. M. Peña, C. Silvia, D. Enrico, J. B. Andrew, P. Rachna, P. Martino, W. B. Gordon, H. B. Asa, and T. Gianluca, “Effect of SR-micro CT radiation on the mechanical integrity of trabecular bone using in situ mechanical testing and digital volume correlation,” J. Mech. Behav. Biomed. Mater. **88**, 109–119 (2018). [CrossRef]

**9. **M. Chang, L. Li, Z. Chen, Y. Xiao, L. Zhang, and G. Wang, “A few-view reweighted sparsity hunting (FRESH) method for CT image reconstruction,” J. X-Ray Sci. Technol. **21**(2), 161–176 (2013). [CrossRef]

**10. **J. Bian, J. H. Siewerdsen, X. Han, E. Y. Sidky, J. L. Prince, C. A. Pelizzari, and X. Pan, “Evaluation of sparse-view reconstruction from flat-panel-detector cone-beam CT,” Phys. Med. Biol. **55**(22), 6575–6599 (2010). [CrossRef]

**11. **S. O. Jin, J. G. Kim, S. Y. Lee, and O. K. Kwon, “Bone-induced streak artifact suppression in sparse-view CT image reconstruction,” Biomed. Eng. Online **11**(1), 44 (2012). [CrossRef]

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

**13. **S. A. Melli, K. A. Wahid, P. Babyn, D. M. L. Cooper, and A. M. Hasan, “A wavelet gradient sparsity based algorithm for reconstruction of reduced-view tomography datasets obtained with a monochromatic synchrotron-based X-ray source,” Comput. Med. Imaging Graph. **69**, 69–81 (2018). [CrossRef]

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

**15. **E. Y. Sidky, C.-M. Kao, and X. Pan, “Accurate image reconstruction from few-views and limited-angle data in divergent-beam CT,” J. X-Ray Sci. Technol. **14**(2), 119–139 (2006).

**16. **E. Y. Sidky and X. Pan, “Image reconstruction in circular cone-beam computed tomography by constrained, total-variation minimization,” Phys. Med. Biol. **53**(17), 4777–4807 (2008). [CrossRef]

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

**18. **Y. Liu, J. Ma, Y. Fan, and Z. Liang, “Adaptive-weighted total variation minimization for sparse data toward low-dose X-ray computed tomography image reconstruction,” Phys. Med. Biol. **57**(23), 7923–7956 (2012). [CrossRef]

**19. **P. Perona and J. Malik, “Scale-Space and Edge Detection Using Anisotropic Diffusion,” IEEE Trans. Pattern Anal. Mach. Intell. **12**(7), 629–639 (1990). [CrossRef]

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

**21. **S. Li, Z. Dong, Q. Gan, J. Song, and Q. Yang, “An adaptive regularized iterative FBP algorithm with high sharpness for irradiated fuel assembly reconstruction from few projections in FNCT,” Ann. Nucl. Energy **145**, 107515 (2020). [CrossRef]

**22. **C. E. Shannon, “IEEE Xplore Abstract - A mathematical theory of communication,” Bell Syst. Tech. J. **27**(4), 623–656 (1948). [CrossRef]

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

**24. **X. Wu, H. Liu, and A. Yan, “X-ray phase-attenuation duality and phase retrieval,” Opt. Lett. **30**(4), 379–381 (2005). [CrossRef]

**25. **M. Ferraro, G. Boccignone, and T. Caelli, “Entropy-based representation of image information,” Pattern Recognit. Lett. **23**(12), 1391–1398 (2002). [CrossRef]

**26. **P. L. Combettes and V. R. Wajs, “Signal recovery by proximal forward-backward splitting,” Multiscale Model. Simul. **4**(4), 1168–1200 (2005). [CrossRef]

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

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

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

**30. **Z. Wang, A. C. Bovik, H. R. Sheikh, and E. P. Simoncelli, “Image quality assessment: from error visibility to structural similarity,” IEEE Trans. Image Process. **13**(4), 600–612 (2004). [CrossRef]

**31. **F. Sattar, L. Floreby, G. Salomonsson, and B. Lovstrom, “Image enhancement based on a nonlinear multiscale method,” IEEE Trans. Image Process. **6**(6), 888–895 (1997). [CrossRef]

**32. **B. P. Fahimian, Y. Zhao, Z. Huang, R. Fung, Y. Mao, C. Zhu, M. Khatonabadi, J. J. Demarco, S. J. Osher, and M. F. Mcnitt-Gray, “Radiation dose reduction in medical x-ray CT via Fourier-based iterative reconstruction,” Med. Phys. **40**(3), 031914 (2013). [CrossRef]

**33. **C. M. Li, W. P. Segars, J. Y. Lo, A. I. Veress, and J. T. D. Iii, “Three-dimensional computer generated breast phantom based on empirical data,” Proc. SPIE **6913**, 691314 (2008). [CrossRef]

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

**35. **R. Liu, M. K. Kalra, J. Hsieh, and H. Yu, “Evaluation of GPU-based CT reconstruction for morbidly obese patients,” JSM. Biomed. Imaging. Data. Papers **4**(1), 1008 (2017).