Three-dimensional visualization of microvasculature from few-projection data using a novel CT reconstruction algorithm for propagation-based X-ray phase-contrast imaging

: Propagation-based X-ray phase-contrast imaging (PBI) is a powerful nondestructive imaging technique that can reveal the internal detailed structures in weakly absorbing samples. Extending PBI to CT (PBCT) enables high-resolution and high-contrast 3D visualization of microvasculature, which can be used for the understanding, diagnosis and therapy of diseases involving vasculopathy, such as cardiovascular disease, stroke and tumor. However, the long scan time for PBCT impedes its wider use in biomedical and preclinical microvascular studies. To address this issue, a novel CT reconstruction algorithm for PBCT is presented that aims at shortening the scan time for microvascular samples by reducing the number of projections while maintaining the high quality of reconstructed images. The proposed algorithm combines the ﬁltered backprojection method into the iterative reconstruction framework, and a weighted guided image ﬁltering approach (WGIF) is utilized to optimize the intermediate reconstructed images. Notably, the homogeneity assumption on the microvasculature sample is adopted as prior knowledge, and therefore, a prior image of microvasculature structures can be acquired by a k-means clustering approach. Then, the prior image is used as the guided image in the WGIF procedure to eﬀectively suppress streaking artifacts and preserve microvasculature structures. To evaluate the eﬀectiveness and capability of the proposed algorithm, simulation experiments on 3D microvasculature numerical phantom and real experiments with CT reconstruction on the microvasculature sample are performed. The results demonstrate that the proposed algorithm can, under noise-free and noisy conditions, signiﬁcantly reduce the artifacts and eﬀectively preserve the microvasculature structures on the reconstructed images and thus enables it to be used for clear and accurate 3D visualization of microvasculature from few-projection data. Therefore, for 3D visualization of microvasculature, the proposed algorithm can be considered an eﬀective approach for reducing the scan time required by PBCT.

The standard filtered backprojection (FBP) is an analytical CT reconstruction method based on the Fourier-slice theorem, and it typically consists of ramp filtering and backprojection procedures, where ramp filtering is usually used to avoid the blur effect on the reconstructed image that arises from the backprojection procedure. Due to its simplicity and high computational efficiency, the FBP algorithm has been commonly utilized to reconstruct CT images in SR-PBCT. According to the Nyquist-Shannon sampling theory, the minimal number of tomographic projections required for high-quality FBP reconstruction is approximately πd /(2r), where d is the maximum thickness of the sample and r is the detector pixel pitch. Under the condition of the minimum projection (CMP), when using FBP reconstruction in SR-PBCT, the high-resolution 3D visualization of microvascular networks typically requires a very long scan time (from tens of minutes to a few hours) for the acquisition of sufficient projections. However, the long scan time might impede the wider use of SR-PBCT in microvascular imaging for biomedical and preclinical experimental studies. First, the long scan time will result in high radiation doses, which not only poses 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 [17]. Second, the long scan time will cause motion artifacts, which arise from not only the heartbeat and breathing of living subjects during the scan but also the deformation of fresh tissue without chemical fixatives that occurs when the scan time is more than approximately 15 minutes [4,18]. To shorten the scan time, one possible solution is to reduce the number of projections required for CT reconstruction. However, few-projection sampling will make the projection data nonuniform, and this can introduce streaking artifacts into the CT image. The streaking artifacts can severely degrade the image quality, which may affect subsequent 3D visualization and quantification of the microvasculature. In fact, ramp filtering in the FBP algorithm, regardless of whether other window functions are added, typically has a limited capability to reduce the streaking artifacts.
Recently, many CT reconstruction methods for the few-projection data have been developed. These methods can be mainly divided into two categories. The first category is the deep-learning (DL) based reconstruction methods. For instance, Han et al. [19] proposed a deep residual learning network based reconstruction method via the persistent homology analysis; Jin et al. [20] developed a deep convolutional neural network based reconstruction method that combined the FBP and U-net; Zhang et al. [21] presented a reconstruction method based on combination of DenseNet and Deconvolution. The more information can refer to a review on the subject of applying DL [22]. The DL based reconstruction methods can perform better than many conventional CT approaches, however, they generally rely on the high-quality big data and meticulos training, and this might be difficult to achieve in some special fields, such as the PCI field. The second category is the iterative reconstruction (IR) methods. The IR methods can utilize prior knowledge to reduce noise or artifacts. According to the employed prior knowledge, the IR methods can be grouped into two types: empirical-knowledge and prior-image knowledge based methods. The empirical-knowledge based methods can utilize the known information (e.g., the boundary, shape, density range and the sparsity) of an object to formulate the reconstruction model, such as the total variation (TV) based reconstruction algorithm [23], tensor based reconstruction algorithm [24], low-rank based reconstruction algorithm [25], and so on. The prior-image knowledge based methods explore both image sparsity and similarity by utilizing the high-quality prior images, e.g., the dictionary learning based reconstruction algorithm [26], the prior image constrained compressed sensing (PICCS) algorithm [27], and the guided image filtering (GIF) [28] based simultaneous algebraic reconstruction technique algorithm [29]. The IR methods have the excellent performance for the few-projection CT reconstruction, however, they often require large memory space for very large computational tasks, and they also take a long execution time (typically costing tens of minutes and even a few hours for one CT). Moreover, they commonly contain many parameters; these parameters may be difficult to adjust, and in many cases, the optimal parameter setting for one reconstructed slice cannot show good robustness for all slices of one sample. In addition, for the prior-image knowledge based methods, they typically need the aid of previously scanned high-quality CT images, and these prior images may not be available in practice. The abovementioned problems might impede the wider use of IR methods in SR-PBCT. Therefore, a novel SR-PBCT reconstruction algorithm for 3D visualization of microvasculature by means of few-projection data, which enables combining the advantages of the FBP algorithm (simplicity and high computational efficiency) and IR algorithms (utilize prior knowledge for processing artifacts in iteration), is highly desirable.
The GIF is a guided image-based filtering technique that utilizes the edge information of the guided image to preserve the edges of the processing image, but it suffers from halo artifacts [28]. Li et al incorporated an edge-aware weighting into the GIF to form a WGIF, and it can avoid halo artifacts in the edge-preserving filtering process [30]. Due to its superior performance and fast processing speed, the WGIF has played an important role in many applications, such as noise reduction, detail smoothing and image matting. Inspired by the abovementioned work, we incorporated FBP reconstruction into the iterative reconstruction framework and adopted the WGIF to optimize the intermediate reconstructed images. It is noteworthy that one approximation on biological soft tissue containing microvasculature was utilized, i.e., biological soft tissue can be approximated as a single-material object. Based on this approximation, a binarized image of microvasculature, namely, the prior image of microvasculature structures, can be accurately obtained by a k-means clustering method [31] and then used as the guided image in the WGIF procedure to effectively suppress streaking artifacts and preserve microvasculature structures, in which no previous high-quality prior images are needed. To assess the effectiveness and capability of the proposed algorithm, simulation experiments on the 3D microvasculature numerical (3D-MV) phantom and SR-PBCT reconstruction experiments on a real microvasculature sample were performed.

Phase retrieval for PBI
An object illuminated by spatially coherent X-ray can be described with its 3D complex refractive index distribution: n(x, y, z) = 1 − δ(x, y, z) + iβ(x, y, z).
where δ(x, y, z) represents the phase information, and β(x, y, z) represents the absorption information; (x, y) denotes the coordinates of a projection plane orthogonal to the propagation direction, and z denotes the propagation axis. In PBI, as the spatially coherent X-ray beams propagate from an illuminated object, the density variations of the object cause phase shifts of X-rays. Owing to Fresnel diffraction, the phase shifts are converted into detectable intensity modulations and subsequently recorded by the detector set at a specific distance. In general, the raw projections from PBI consist of both absorption information and the Laplacian of phase information. To obtain the phase-information distribution of the illuminated object, we first needed to extract the phase-shift information from the raw projections, and thus, phase retrieval was required [32]. To acquire the accurate phase shift for each tomographic projection, phase retrieval generally needs at least two phase-contrast projections taken at different SDDs. However, the multiple-SDD phase retrieval methods are usually difficult to arrange experimentally due to the complicated registration problem and the long execution time. In this work, we assumed that the microvasculature sample could be approximated as a single-material object; namely, the object was 'homogeneous', and thus, a so-called 'homogeneous' transport of intensity Equation (TIE-Hom) phase retrieval method [33] that required just a single-SDD projection was used to accurately extract the phase-shift information. The TIE-Hom method can be described by the following formula: where ϕ θ (x, y) defines the phase-shift function at the projection angle θ, and I θ,z=SDD (x, y) denotes the intensity function in the detector at a distance D; F and F −1 are the 2D FFT and its inverse operators, respectively; γ represents the ratio between δ(x, y, z) and β(x, y, z), and λ denotes the wavelength of the X-ray; ξ and η are the coordinates of x and y in the Fourier domain, respectively. In the near-field regime of PBI, the refraction and scattering effects between the X-ray and illuminated object are weak, and hence, the X-ray beam propagation path inside the illuminated object can be assumed to be straight. Therefore, the phase-shift function of the sample can be written as follows: where L θ denotes the linear propagation path inside the sample at the projection angle θ, and δ(x, y, z) denotes the 3D phase-information distribution of the sample. Generally, the SR X-ray source has a monochromatic property, and thus, the forward projection model of phase shift in SR-PBCT, namely, the phase-shift set belonging to the same slice (y = constant) that consists of ϕ θ (x) ranging from 0°to 180°, can be approximately expressed as the following discrete linear system: where A is the system matrix used for modeling the parallel-beam X-ray forward projections [34], and P denotes (λ/2π) times the phase-shift set.

The iterative reconstruction strategy based on the FBP algorithm
Following phase retrieval, according to Eq. (2) and Eq. (3), the 2D phase-information distribution of the sample, i.e., the SR-PBCT images of the sample, can be acquired by CT reconstruction.
Here, the solution for the reconstruction problem of the 2D phase-information distribution of the sample using FBP was formulated as follows: where δ(x, z) denotes the SR-PBCT image, P θ (ω) represents the Fourier transform of (λ/2π)ϕ θ (x), and |ω| defines the ramp filter in the Fourier domain. For simplicity, a so-called FBP operator was utilized, and Eq. (5) can be rewritten as following: where δ is an abbreviation for δ(x, z).
In this study, an iterative reconstruction strategy based on the FBP algorithm [35] was used to solve the inverse problem in Eq. (4), and the solution can be described as follows: where α is the step size, and the superscripts t and t + 1 represent the t-th iteration and (t + 1)th iteration of the recursive scheme, respectively.

The WGIF algorithm based on a prior image
In practice, the solution expressed in Eq. (7) would lead to poor performance for the few-projection CT reconstruction problem because it neither builds an effective artifacts model nor utilizes prior knowledge for artifacts suppression in the reconstruction strategy. Generally, the biological soft tissue surrounding the microvasculature can be considered a single-material object, i.e., the biological soft tissue is 'homogeneous'. Moreover, owing to the phase-based contrast generation mechanism, the interfaces between the biological soft tissue and microvasculature typically show a high image contrast. Therefore, a binarized image of microvasculature, where the flat regions represent prior knowledge of the homogeneous biological soft tissue and the background (typically the air) for smoothing and, additionally, the interfaces between the highlight and shadow regions represent prior knowledge of the microvascular structures for preserving, can be accurately generated by the k-means clustering approach. In this study, the binarized image of the microvasculature was considered the prior image. The GIF is a local filtering-based edge-preserving smoothing technique, and it is assumed to be a local linear model between the guided image G and the filtering output ∧ δ; namely, the abovementioned prior image can be used as the guided image of the microvasculature image. In principle, the local linear model between G and ∧ δ can ensure that the output ∧ δ has an edge only if the guided image G has an edge, and it can be formulated as follows: where a k and b k are linear coefficients assumed to be constant in window ω k centered at pixel k, and l is the pixel index.
To determine the linear coefficients a k and b k , a constraint is added to the filtering input δ and the filtering output ∧ δ; the constraint is expressed as follows: where the output ∧ δ represents homogeneous areas with sharp edges, and e l represents unwanted components, such as noise or artifacts.
Based on the constraint in Eq. (9), the values of a k and b k can be determined by minimizing a cost function E(a k , b k ) in the window ω k , which is defined as follows: where ε is a regularization parameter that penalizes the large a k . In general, the cost function E(a k , b k ) can be utilized to minimize the difference between the filtering input δ and the filtering output ∧ δ while maintaining the linear model in Eq. (8). However, as a local optimization-based function, E(a k , b k ) cannot preserve sharp edges, and it usually produces halo artifacts near some edges. To overcome this deficiency, the WGIF was presented that incorporated an edge-aware weighting into the GIF [30]. The WGIF computed the edge-aware weightings in the guided image G, and then, these weightings were used for the better preservation of edges in the filtering input δ. Here, the edge-aware weighting W G k is defined by utilizing local variances of 4×4 windows of all pixels as follows: where M represents the total number of pixels in the guided image G, and ε 0 represents a small number used to avoid singularities. In this study, ε 0 is set to 1.0×10 −6 ; (σ G,1 k ) 2 defines the variance in the guided image G in the 4×4 window centered at the pixel k, and (σ G,1 l ) 2 represents the variance in the guided image G in the 4×4 window centered at pixel l.
According to Eq.(10) and Eq.(11), the cost function ∧ E(ak, b k ) of the WGIF can be written as follows: The ∧ E(ak, b k ) is a linear ridge regression model [36], and its solution can be given as follows: where u k and σ k define the mean value and variance, respectively, of the guided image G in window ω k ; |ω k | represents the number of pixels in window ω k , and additionally, Following that a k and b k were computed for all windows ω k in the filtering input δ, the final filtering output ∧ δ can be obtained as follows: where _ a l = 1 the mean values of all windows overlapping l.

Steps of the proposed CT reconstruction algorithm
The proposed algorithm mainly consists of two steps per iteration: the FBP step and the WGIF step. Before the FBP step, the two images were first acquired: one was the CT reconstruction image using FBP from the inpainted projection data that was obtained by interpolation between undersampled angles based on CMP, and the other was the few-projection CT reconstruction image using FBP that was processed by the WGIF. In the filtering process of the WGIF, the binarized image of the few-projection CT reconstruction image generated by the k-means clustering method was used as the guided image. Considering that the guided image was not updated, it was named a static guided image. By fusing the abovementioned two images, a good initial guess was obtained, and this was important for the convergence of the proposed algorithm. Then, the FBP step was used to acquire few-projection CT reconstruction images with artifacts. Following the FBP step, the WGIF step was performed to reduce the artifacts in CT reconstruction images obtained by the FBP step and improve image quality. It is worth mentioning that the guided image was updated per iteration in the WGIF step, which was termed a dynamic guided image. Subsequently, the resultant image from the WGIF step was reprojected based on CMP to produce the forward projection data, and then, the projection residual was modeled as the inpainted projection data subtracting the forward projection data. Finally, the projection residual flowed into the FBP step for the next iteration, and this procedure was updated based on Eq. (7). The FBP step and the WGIF step alternated during the iteration until the stop criterion was met. The steps of the proposed algorithm are better illustrated in Fig. 1.  Fig. 1. Flowchart of the proposed CT reconstruction algorithm. The symbol + represents the addition operator, the symbol − represents the subtraction operator, and the arrows with dotted lines represent iterative processes.

Pseudocode for the proposed CT reconstruction algorithm
To explain the proposed algorithm more clearly, the corresponding pseudocode is presented, as shown in Algorithm 1.

Quantitative assessment
In this study, the mean structural similarity index (MSSIM) [37], signal-to-noise ratio (SNR) [38], relative loss rate and relative difference (RD) were used as quantitative metrics. The MSSIM is widely utilized to evaluate the similarity between a reconstructed image and a reference image, and it is usually between 0 and 1, which increases with increasing similarity. The SNR is a traditional measure of image quality, and a larger SNR value indicates better image quality. The relative loss rate is used to measure the relative volume loss between the reconstructed microvasculature and the reference microvasculature, and a smaller value indicates more accuracy. The RD is used to assess the difference between two reconstructed images in adjacent iterations, and the smaller RD value represents the smaller difference.
(i) MSSIM is defined as follows, where x and y are the reference image and the reconstructed image, respectively; x j and y j are the jth image patch of x and y, respectively; u x j and u y j are the mean values of x j and y j , respectively; σ 2 x j and σ 2 y j are the variances of x j and y j , respectively; σ x j y j denotes the standard deviation between x j and y j ; c 1 and c 2 represent two constants to stabilize the division with weak denominator, i.e., c 1 = (k 1 L) 2 , c 2 = (k 2 L) 2 ; T is the total number of image patch. In this study, the T was set to 729, L was set to 255, k 1 and k 2 was set to 0.01 and 0.03, respectively; (ii) SNR is defined as follows, where S sample denotes the mean value of the pixel values in the local sample region of the reconstructed image y, and σ background the standard deviation of the pixel values in the local background region of the reconstructed image y.
(iii) Relative loss rate can be defined as follows, where Vessel GT denotes the microvasculature volume calculated from the reference images, and Vessel R represents the microvasculature volume calculated from the reconstructed images.
(iv) RD can be defined as follows, where y t and y t+1 denote the reconstructed images in t-th and (t + 1)-th iterations, respectively.

Simulations
To assess the performance of the proposed algorithm, a 3D-MV phantom designed via the real SR-PBCT images of the microvasculature in a mouse-liver sample was employed, as shown in Fig. 2. The 3D-MV phantom consists of three different types of compositions, such as air, teflcn and water, where the air and teflcn parts located in the field of view (FOV) were used to simulate the parts inside and outside the microvasculature and, additionally, the water located outside the FOV was used as the background for computing the SNR. The pixel values of air, teflcn and water in the 3D-MV phantom were set to 0, 7.02×10 −7 , and 3.69×10 −7 , respectively, which represented the corresponding phase values of the abovementioned three types of compositions at an energy of 25 keV (the energy dependence of the complex refractive index can be obtained for a variety of chemical compounds under http://henke.lbl.gov/optical_constants/getdb2.html). The volume of the 3D-MV phantom was 216×216×100 voxel 3 , and the size of the voxel pitch was 3.25×3.25×3.25 um 3 . The detector is modeled as a straight-line array with 306 bins, and the size of the reconstructed images was 216×216 pixel 2 . The total number of reconstructed slices was 100. According to CMP, the number of the full projection required for FBP was approximately 340. In this study, 68 and 34 uniformly distributed projections over acquisition orbital extents (θ) ranging from 0°to 180°(full-scan range) were used as few-view projections. The experiments with CT reconstruction were performed using the MATLAB programming language on a workstation computer equipped with an Intel (R) Xeon (R) E5-2640 CPU at 2.4 GHz and 128 GB of RAM. Finally, the 3D surface-rendering images were generated from a sequential series of CT slices using Amira 6.3.0 software (Visage Imaging, Berlin, Germany).

Noise-free reconstructed results
In this section, the optimal parameters for the proposed algorithm, such as α and ε, were determined by comparing both the quantitative metrics and the visual assessment of the reconstructed images from the different parameter values, and hence, three sets of parameters were obtained as follows: (i) one set for 340 projections: α = 1.2, ε = 0.01; (ii) another set for 68 projections: α = 1.0, ε = 0.01; and (iii) the other set for 34 projections: α = 0.4, ε = 0.01. The total iteration number t max was set to 5, which was utilized as the stop criterion for the proposed algorithm. Figure 3 shows the surface-rendering images of the 3D-MV phantom in which all slices were reconstructed using the FBP and the proposed algorithm under the noise-free 340 projections, 68 projections and 34 projections. Additionally, Fig. 3 also shows the 24th reconstructed slices using the FBP and the proposed algorithm under the noise-free 340 projections, 68 projections and 34 projections, which were used for comparing tomographic sections. Figure 3(b) and Fig. 3(d) clearly show that the visual effect of the CT image reconstructed using the proposed algorithm under the 340 projections was almost the same as that reconstructed using the FBP, and the same was true for their corresponding 3D images (see Fig. 3(b) and Fig. 3(d)). As the number of projections decreased, the streaking artifacts on the CT images reconstructed using the FBP increased (see Fig. 3(f) and Fig. 3(j)), and the false structures on the corresponding 3D images, which are typically formed by the streaking artifacts on the CT images, also increased (see Fig. 3(e) and Fig. 3(i)). Obviously, compared with the reconstructed results of the FBP under the few-projection conditions, those of the proposed algorithm showed a better visual effect, which suggests that the proposed algorithm has the ability to reduce the streaking artifacts on the CT images (see Fig. 3(h) and Fig. 3(l)) and thereby enables the suppression of the false structures on their corresponding 3D images (see Fig. 3(g) and Fig. 3(k)). To further assess the performance of the proposed algorithm, the profiles along the same position in the reconstructed images, the position labeled with the red line in Fig. 3(b), were drawn, as shown in the upper right corner of the CT images in Fig. 3. When the number of projections decreased, the profiles of the reconstructed images using the FBP began fluctuating in the homogeneous areas and near the boundaries (see Fig. 3(b, f, j)), which means that the observation for some boundaries and fine features were affected. However, as the number of projections decreased, the profiles of the reconstructed images using the proposed algorithm remained smooth in the homogeneous areas and near the boundaries (see Fig. 3(d, h, l)), which demonstrates that the proposed algorithm has the ability to suppress the streaking artifacts in the homogeneous areas and preserve the boundaries (i.e., the microvasculature structures).

Noisy reconstructed results
In this study, a combination of Poisson noise and Gaussian noise [39], which was used to simulate the statistical noise in the projections, was added into the phase-shift sets for analyzing the robustness and reliability of the proposed algorithm in the presence of statistical noise. The parameters of the noise model were set as follows: the photon number of Poisson noise I 0 was set to 1.0×10 5 , while the mean value m G and variance σ 2 G of the Gaussian noise were set to 0 and 5, respectively. In this experiment, the optimal values of the parameters, such as α and ε, were selected by comparing both the quantitative metrics and the visual assessment of the reconstructed images based on the different parameter values. Finally, three sets of parameters were acquired as follows:   projections and 34 projections, as shown in Fig. 4(a, e, i) and Fig. 4(c, g, k). In addition, Fig. 4 also depicts the 24th reconstructed slices using the FBP and the proposed algorithm under the noisy 340 projections, 68 projections and 34 projections, as shown in Fig. 4(b, f, j) and Fig. 4(d,  h, l). Compared with the reconstructed slices using the FBP under the noisy few-projection conditions (see Fig. 4(b, f, j)), the reconstructed slices using the proposed algorithm under the noisy 340 projections, 68 projections and 34 projections (see Fig. 4(d, h, l)), had less noise and fewer artifacts, which achieves a better visual effect for observing the microvasculature structures. Obviously, it can be found that the surface-rendering images of the 3D-MV phantom where all slices were reconstructed using the FBP under the noisy few-projection conditions (see Fig. 4(a, e, i)) were severely affected by considerable false structures, which led to a poor visual effect for distinguishing the microvasculature structures. However, the surface-rendering images of the 3D-MV phantom where all slices were reconstructed using the proposed algorithm under the noisy few-projection conditions (see Fig. 4(c, g, k)) showed a superior visual effect for observing the microvasculature structures, where the false structures were effectively suppressed. Additionally, from the profiles of the reconstructed images in Fig. 4(b, f, j) and Fig. 4(d, h, l), it can be clearly seen that the profiles of the latter are much smoother in the homogeneous areas and nearby the boundaries than those of the former, which demonstrates that the proposed algorithm enables suppression of the noise and artifacts, as well as enables preservation of the boundaries.

Assessments
In this section, both the image quality and the accuracy of the 3D-MV phantom where all slices were reconstructed using the FBP and the proposed algorithm were assessed via quantitative metrics, such as MSSIM, SNR and relative loss rate. The quantitative metrics for the 3D-MV phantom were measured based on all slices (100 slices), where all slices were reconstructed using the FBP and the proposed algorithm under the noise-free and noisy few-projection conditions, respectively. The measured quantitative metrics are presented as the mean ± standard deviation, as displayed in Fig. 5. Under the noise-free and noisy few-projection conditions, the MSSIM values and the SNR values from the reconstructed slices using the proposed algorithm are larger than those of the reconstructed slices using the FBP. The higher MSSIM and SNR values confirmed that the proposed algorithm effectively reduced the noise and artifacts.
Due to the existence of false structures, it is difficult to accurately measure the microvasculature volume. Therefore, both the reconstructed microvasculature without the false structures and the reference microvasculature (see Fig. 2(a)) were used to compute the relative loss rate, and the corresponding results are shown in Fig. 5(e, f). From Fig. 5(e, f), it can be observed that the proposed algorithm remains small relative loss rate values (<3%) of microvasculature structures under the noise-free and noisy few-projection conditions, which demonstrates that the proposed algorithm enables a high reconstruction accuracy for microvasculature.

The influences of the step size α on the reconstructed images
In this section, to analyze the influences of the step size α, as depicted in Eq. (7), the different reconstructed images based on α values ranging from 0 to 1.8 were evaluated by the MSSIM, as shown in Fig. 6. From Fig. 6, it can be observed that the noise and the artifacts in the reconstructed images increased as the α values increased, and the edges (i.e., the microvasculature structures) were more accurate and clearer. However, if the α values became too small, the imaged edges and fine features were severely deformed, and if the α values became too large, the imaged edges and fine features were severely destroyed by the noise and the artifacts. Therefore, the step size α can serve as a balance parameter that controls the trade-off between the structure fidelity of the edges or fine features and the smoothness of the artifacts or noise. Since the step size α plays a critical role in the CT reconstruction process, the selection of the optimal α will be important. In this study, the optimal α values were selected by comparing both the MSSIM and the visual assessment of the reconstructed images based on different α values.

Convergence analysis and computational time
To assess the convergence performance of the proposed algorithm under the noise-free and noisy few-projection conditions, the MSSIM-based and SNR-based convergence curves were presented, as shown in Fig. 7. The proposed algorithm converged near the 5th iteration. Therefore, in this study, the total number of iterations t max was set to 5, and it was used as the stop criterion for the proposed algorithm. In addition, the reconstructed time of all slices from the 3D-MV phantom using both the FBP and the proposed algorithm under the noise-free and noisy few-projection conditions were computed, as presented in Table 1. The computational time of the proposed algorithm was approximately two minutes, which was much longer than that of the FBP, but this was acceptable for 3D reconstruction tasks.

Sample preparation
The animal experiments were performed in accordance with the guiding principles for the care and use of laboratory animals approved by the Research Ethics Committee of the Beijing Friendship Hospital, Capital Medical University, China. In this study, a healthy female C57BL/6 mouse (7 weeks old, weighing 18-20 g) used for SR-PBCT experiment was euthanized with a dose of 800 mg/kg of sodium pentobarbital, and then an entire liver lobe was removed from the mouse. After removal, the entire liver lobe was rinsed by phosphate buffer and perfused with saline 0.9% to clear the blood vessels, finally fixed in 10% neutral buffered formalin solution for 4 months. Before imaging, the entire liver lobe was dehydrated using an ethanol series (50%, 70%, 80%, 95%, 95%, 100% and 100%), and then it was transferred to the imaging facility in motion-proof container to avoid risk of tissue damage during shipping.

PBI data acquisition
The PBI data of the entire liver lobe were collected at the BL13W1 beamline at the Shanghai Synchrotron Radiation Facility, Shanghai, China. In this experiment, the energy of the incoming monochromatic photons was set as 14 keV, and the SDD was adjusted to 0.2 m. A charge-coupled device camera (Photonic Science, UK) with a 6.656 mm (horizontal)×4.264 mm (vertical) FOV was adopted as the imaging detector, and its effective pixel pitch was 3.25×3.25 µm 2 . According to CMP, the reference number of the projection required for FBP was approximately 1740. For the SR-PBCT imaging, 1200 projections at equidistant angles from 0°to 180°were collected. The size of the projection image was 2048×1312 pixel 2 , and the exposure time per projection was 1.5 s. The total scan time was approximately 40 minutes. In addition, ten dark current images (without the X-ray beams) were utilized to reduce the dark signals in projections, and twenty flat-field images (with the X-ray beams on, but without the sample) were utilized to correct the pixel-to-pixel nonuniformity in projections [40]. In the few-projection SR-PBCT experiments, 400 and 200 projections were uniformly chosen from 1200 projections recorded by the detector. The phase retrieval was performed using the TIE-Hom method, and the γ was set to 200. After phase retrieval, 1200, 400 and 200 projections were used for the few-projection CT reconstruction. Finally, 550 slices containing the microvasculature information were chosen from the 1312 reconstructed slices and then used for the 3D visualization of the microvasculature.

Experimental results
In this experiment, the optimal parameters for α and ε were selected by comparing both the quantitative metrics and the visual effect of the reconstructions from different parameter values. Finally, three sets of parameters were acquired as follows: (i) one set for 1200 projections: α = 1.0, ε = 0.10; (ii) another set for 400 projections: α = 1.0, ε = 0.13; and (iii) the other set for 200 projections: α = 0.8, ε = 0.14. The total number of iterations t max was set to 8. Figure 8 shows the surface-rendering images and the 318th slices of the entire liver lobe in which all slices were reconstructed using the FBP and the proposed algorithm from 1200 projections, 400 projections and 200 projections. Figure 8(a, e, i) depict the surface-rendering images of the entire liver lobe in which all slices were reconstructed using the FBP from 1200 projections, 400 projections and 200 projections, respectively. As the number of projections decreased, it can be obviously observed that the false structures on the corresponding 3D images increased, and the same was true for the streaking artifacts on the corresponding CT images (see Figs. 8(b, f, j)). The 3D reconstruction process for microvasculature in which all slices were reconstructed using the FBP approach from 200 projections can be observed in Visualization 1, where the formation process of the false structures can also be seen. Figure 8(c, g, k) depict the surface-rendering images of the entire liver lobe in which all slices were reconstructed using the proposed algorithm from 1200 projections, 400 projections and 200 projections, respectively. Compared with Figs. 8(a, e, i), Figs. 8(c, g, k) have clearer 3D microvasculature structures, where the false structures were greatly reduced, and thus enabled the subsequent structure segmentation, measurement and analysis for the 3D microvasculature. Additionally, compared with Figs. 8(b, f, j), Figs. 8(d, h, l) presented the smoother tissue areas and clearer edge details, which indicated that the proposed algorithm was capable of suppressing the streaking artifacts and preserving the microvasculature structures under the few-projection condition. The 3D reconstruction process for microvasculature in which all slices were reconstructed using the proposed approach from 200 projections can be observed in Visualization 2, where clear 3D microvasculature structures without false structures can also be observed.

Result analysis
To provide a better visual effect for the result analysis, the areas in the reconstructed images, as marked by the green rectangle in Fig. 8(b), were selected as the regions of interest (ROIs) and then enlarged. The enlarged images of ROIs are presented for comparison, as shown in Figs. 9(a-f). Figures 9(a, b, c) are the enlarged ROIs of the reconstructed images using the FBP approach from 1200 projections, 400 projections and 200 projections, respectively. Figures 9(d, e, f) are the enlarged ROIs of the reconstructed images using the proposed approach from 1200 projections, 400 projections and 200 projections, respectively. Compared with Figs. 9(a, b, c), it is obvious that Figs. 9(d, e, f) provided the more superior visual effect for microvasculature, where the streaking artifacts were effectively smoothed, and the edges of microvasculature were well preserved, while some imaged blood vessels of the former were severely destructed by the streaking artifacts and the subsequent structure segmentation, measurement and analysis for the microvasculature were significantly affected. Moreover, the profiles along the same position in enlarged ROIs, the position marked with the red line in Fig. 9(a), were drawn, as presented in Figs. 9(g, h, i). Figures 9(g, h, i) represent the profiles of the reconstructed images using the FBP approach and the proposed approach from 1200 projections, 400 projections and 200 projections, respectively. From Figs. 9(g, h, i), it can be seen that the profiles from the proposed algorithm are much smoother in the tissue regions (homogeneous areas) and near the edges (boundaries) than those from the FBP approach, which demonstrates that the proposed algorithm is capable of smoothing the streaking artifacts, as well as preserving edges (i.e., the microvasculature structures).

Convergence analysis and reconstructed time
In this section, the quantitative metric RD, as depicted in Eq. (20), was used to evaluate the convergence performance of the proposed algorithm. Figure 10 presents the RD-based convergence curves of the proposed algorithm using 1200 projections, 400 projections and 200 projections. It could be clearly observed that the proposed algorithm converged near the 8th iteration under different few-projection conditions. Thus, in this experiment, the total iteration number t max was set to 8. Additionally, under different few-projection conditions, the  reconstructed time for 550 slices of the entire liver lobe based on the FBP and the proposed algorithm were measured, as shown in Table 2. The reconstructed time of the proposed algorithm was approximately twenty hours, which was acceptable for 3D reconstruction tasks.

Discussion and conclusion
As the widely used FBP algorithm fails to reconstruct high-quality SR-PBCT images of microvasculature under the few-projection conditions, we developed a new few-projection SR-PBCT reconstruction algorithm for microvasculature imaging. In the experiments, the 3D-MV phantom and SR-PBI data of the mouse-liver sample were utilized to assess the accuracy and feasibility of the proposed algorithm, and the FBP algorithm was used for comparison. The simulation experiment results confirmed that the proposed algorithm was an effective method for suppressing the streaking artifacts as well as preserving the edges under the few-projection conditions and that it had robustness and reliability in the presence of statistical noise. The SR-PBCT experiment results demonstrated that the proposed algorithm had the ability to effectively suppress the streaking artifacts and well preserve the edges of microvasculature under the few-projection conditions, and the same was true for the corresponding 3D images, where the false structures could be effectively suppressed, and the 3D microvasculature structures could be clearly presented. Overall, the proposed algorithm is significant for SR-PBCT reconstruction of microvasculature from few-projection data, and it is also practical for the 3D visualization of microvasculature since the data acquisition time is significantly reduced. Generally, SR-PBCT enables high-resolution and high-contrast 3D visualization of microvasculature and has the potential to serve as an adjunct to histopathological sections for observing microvascular changes. However, the main limitation of the SR-PBCT method is that the SR source requires an expensive synchrotron facility, which impedes its wider use in preclinical and biomedical microvascular imaging studies. Recently, transferring the PBCT method onto laboratory X-ray sources, such as microfocus X-ray sources, has been an active area of research. Compared with the SR source, the laboratory X-ray sources typically have lower photon flux [9,41], which results in longer scan times for collecting data (typically a few hours). For laboratory X-ray source-based PBCT technology, few-projection PBCT reconstruction algorithms will also be highly desirable, and they have great potential to make both fresh and chemically unfixed biological tissues available for this technology [4]. Additionally, the laboratory X-ray source-based PBCT technology may also be useful for the intraoperative scanning of excised microvascular tissues and even available for in vivo microvascular imaging of small animals, where the few-projection PBCT reconstruction algorithms will have important significance due to the long scan time and the high radiation dose. In summary, we believe that this work will facilitate the development and application of PBCT technology in microvasculature imaging.
In this study, the homogeneity assumption on the microvasculature sample was utilized. Based on this approximation, accurate phase information was obtained by the TIE-Hom method, and a binarized image acquired by the k-means clustering method was used as the prior image of microvasculature structures. For some multimaterial samples, the TIE-Hom method would still be applicable if the regions of different materials could be approximated as homogeneous, and the k-means clustering method would also be available in this case. Thus, the proposed algorithm has the potential to reconstruct other single-material or multimaterial samples, and this will be investigated in our future work. Furthermore, GPU-based parallel computing techniques will be implemented to improve the reconstructed speed of the proposed algorithm, and the adaptivity of the parameters in the proposed algorithm will also be researched. In addition, future studies will also be performed to test whether the proposed algorithm applies for laboratory X-ray source-based PBCT technology and in vivo imaging data.

Funding
National Natural Science Foundation of China (81671683, 81670545, 81371549); Natural Science Foundation of Tianjin City in China (16JCYBJC28600); The Foundation of Tianjin university of technology and education (KJ12-01, KJ17-36).