1 Introduction

Nowadays, liver fibrosis and its end-stage, cirrhosis, become a major threat to the public’s health. It occurs in many chronic hepatic diseases and may lead to severe complications ^{[1, 2]}. The evaluation of the hepatic fibrosis is critical in the diagnosis and treatment of such hepatic disease. Angiogenesis is one of the prominent characteristic of hepatic fibrosis ^{[3, 4]}. With the development of imaging techniques, particularly the synchrotron radiation based micro-tomography (SR-μCT) ^[5], hepatic vasculature imaging and its quantification are widely used in studying the development of hepatic fibrosis^{[6, 7]}. PCI(Phase-contrast imaging)-based SR-μCT is able to detect minute density variation between biological soft tissues, which allows inspection of vascular microstructures without contrast agent. However, three dimensional images with elaborate hepatic vasculature introduce challenges to the structural analysis and extraction of quantitative information.

For quantification of angiogenesis to evaluate the fibrosis progression, a series of parameters should be extracted from the vasculature, such as micro-vessel density(MVD), the ratio of extreme vessel to volume(E/V) and the ratio of junction point to volume(J/V) ^{[3, 4, 8, 9]}etc. However, vessels in the liver do not distribute uniformly, so the density descriptors limit accurate assessment of angiogenesis. Due to limitations in the sample preparation and the imaging, the intact vasculature in a specimen may not be imaged integrally and three dimensionally, which means that available parameters for angiogenesis assessment cannot ensure the accuracy and reliability. In this paper, based on the tree structure of hepatic vasculature, a quantification is introduced to assess the degree of angiogenesis. Experiments based on a series of collected and simulated liver specimens at different stages of fibrosis are performed to evaluate practicability of the proposed method.

2 Specimen preparation and image acquisition

2.3 Image reconstruction

All the image datasets are reconstructed by the software PITRE^[13]. Phase retrieval was implemented to each projection image, and the stack image was reconstructed using the filtered back projection algorithm. The output 32-bit stack images are converted to 8-bit so as to conduct further image processing and analysis. There are 8 stack images in total. Fig.1 shows the 3D rendering of all the reconstructed images, some of which miss a portion of the specimen so that the vasculature inside the liver is not intact.

Fig.1Full-screen View

Typical reconstructed images for eight pieces of mouse lever samples at different liver fibrosis stages.

3 Quantitative analysis to the vascular tree

In this section we explain the method of tree analysis for the hepatic vasculature in details so that the acquired quantitative parameters are able to evaluate the degree of liver fibrosis efficiently. The procedures of acquiring quantitative information of the vasculatures from the micro-CT images as illustrated in Fig. 2.

Fig. 2Full-screen View

Flowchart of tree analysis to the hepatic vasculature.

3.1 Image preprocessing

Fig. 3(a) is the raw image after CT reconstruction and image preprocessing. Fig. 3(b) is the liver mask obtained by employing a Gaussian filtering to Fig. 3(a) with a high sigma value. The major outline of the liver is isolated from Fig. 3(a) using the liver mask, Fig. 3(d) shows the image processed further by Ostu algorithm in which the vasculature and image background voxels are shown at the same grey value. Fig. 3(e) is a dedicated mask image used to isolate vasculature, which is generated by the convex hull of the liver shape in Fig. 3(c) and used for morphological 3D erosion to confine the vasculature out of the background if necessary. After all these procedures, the slice image of the isolated vasculature is achieved (Fig. 3f).

Fig. 3Full-screen View

Procedures for isolating vasculature from background.

The particle-removing and cavity-filling procedures are as follows.

Cavity is a region where grey values of all voxels are 0, and it is surrounded by object voxels. There is only one background component in an image. All the cavities are filled once the only background component is segmented by region-growing method. Then, the grey value of the whole 3D image is inversed.

All the object components are labeled and counted by ImageJ plugin “3D Objects Counter” ^[14]. The components of smaller sizes than a threshold are removed. And the binary image containing the solid vasculature is ready for further quantification.

3.3 Quantification and vascular tree analysis

The segmented vasculature and its associated skeleton are prepared to extract the statistical information of the vasculature based on its tree structure. This section describes the quantitative analysis of hepatic vasculature in details. We highlight the analysis of the characteristics of the vascular tree, whose output contains the geometry information of each vascular branch in hierarchy.

3.3.1 The root of vascular tree

This is the start position for tracing the whole tree structure of vasculature. The root is identified automatically. We scan each slice from the first frame to the last one along six directions of the axes in the 3D image space Z³(±X, ±Y, ±Z). In the first frame every object voxel is searched, the largest and central connected region in the 2D frame is selected as the candidate root area. Therefore, there are six candidate root regions in terms of the six scanning directions. Then we calculate the largest inscribed sphere inside the vasculature to obtain the approximate location of the root. The candidate root area being the closest to the largest inscribed sphere is considered as the reference root. If a 3D image has more than one vascular tree, we could label all the connected components to differentiate all trees and identify trees’ roots based on the tree labels. Besides, the manually specified root can also be accepted if this is preferred.

3.3.2 Elements of the vasculature

The 3D vascular tree analysis method is developed from the 2D image analysis of neurite outgrowth ^[19]. This method outputs quantitative information of the vascular trees both on statistics and tree hierarchy images.

Fig. 4 illustrates the tree structure of vasculature. The r₁, s_1–8, b_1–3 and e_1–4 are the root, branches, bifurcations and extremities, respectively. Table 1 lists definitions of the elements in a vascular tree, and the serial properties concerning measurements and structure information.

Table 1Full-screen View

The elements in a vascular tree and the serial properties.

Elements	Serial properties
rm, vasculature root	RL(rm), root label
	RSL(rm), connected branches, as sub-trees
s i,vessel segments	VL (si), branch label
	VLen(si), branch length
	VRad(si),  mean radius of vessel branch
	VLV(si),branch level
	VR(si), associated root label
	VP(si), associated connection points’ labels
	VT(si), associated sub-tree label
b j,bifurcations	BL(bj), connection point label
	BV(bj), associated branches’ labels
e k  , extremities	EL(ek), extremities label
	ELR(ek), length from root to extremity
	EB(ek), associated branch

Show more

Fig. 4Full-screen View

Tree structure of vasculature.

The following values are given to r₁, s₂, b₃ and e₄ in Fig.4.

r₁: RL(r₁) = 1, RST(r₁) = s₁

s₂: VL(s₂)= 2; VLen(s₂)= length of s₂; VRad(s₂)= radius of s₂; VLV(s₂)= 2; VR(s₂)= r₁; VP(s₂)= b₁, b₂; VT(s₂)= sub-tree 1

b₃: BL(b₃)=3: BV(b_j)= s₃, s₆, s₇, s₈

e₄: EL(e₄)=4. ELR(e₄)= VLen(s₁) + VLen(s₃) + VLen(s₇), EB(e₄)=s₇

Then, the vasculature quantification is to traverse the vascular tree and assign property values to each element.

3.3.3 Traversing vascular tree and the related measurement

There are several pre-processing workflows before constructing the vasculature tree. Firstly, a label is given to each root of vascular tree. Next, the skeleton of vasculature is divided into unique vessel segment, by removing the intersection voxel from the skeleton. An intersection point refers to points having more than two 26-connected neighbours. Each segment is assigned with a unique label. Then, all the intersection points are isolated and morphologically dilated by a $\text{3×3×3}$ structure element. Each component is considered, and labelled, as a bifurcation of the vasculature. Finally, tree roots are dilated to overlap the vessel segment connected to the root. These segments are considered as root segment which are assigned the tree level VLV(s_i) = 1, the associated root label VR(s_i) and a sub-tree labels VT(s_i).

Finally, some properties listed in last section can be accumulated before the tree growing process. The branch length VLen(s_i) is assessed by the path from one extremity to the other. The distance between two neighbour voxels is Euclidean distance 1, $\sqrt{2}$ or $\sqrt{3}$ in terms of the connection type, including face-connected, edge-connected or vertex-connected respectively. The associated connection points VP(s_i) are the overlap bifurcations. Each voxel of individual segment is regards as the centre of the maximal inscribed-sphere of the vasculature. Therefore, the mean radius of a vessel VRad(s_i) is computed using the radius of all these maximal inscribed-sphere.

The graph-based watershed methodology ^[20] is used to associate isolated vessel segments to grow vascular tree. Typically, the watershed method is applied for image analysis based pixels which are named as nodes and seeds. In this case, the nodes are the individual vessel segments instead of pixels and seeds are the root segments. Root segments are initially put in a priority queue. Vessel segments are repeatedly taken from the top of the queue and its neighbouring vessel segments are added with priority of mean radius VRad(s_i). Radius is used to determine the priority as it is the reasonable indicator of the tree hierarchy. Once a segment is taken out of the queue, it passes the properties to its neighbours being never added into the queue. Accordingly, the segments inherit the associated root VR(s_i) and the sub-tree label VT(s_i). The same branch level VLV(s_i) is assigned to the segment with the highest priority, while a higher level to the others. Meanwhile, the length to the root is accumulated at this stage. The process that segments are passed in and out of the queue terminates until all the segments are removed out of the queue.

3.3.4 The quantification output

The quantification result is presented in both statistics and image patterns. The image-wide summary of vasculature includes the following statistics:

• Number of roots—number of vascular tree

• The total length of all the vascular branches (in pixels)

• Total number of vascular segments

• The longest length extended from the root

• Total number of extreme vessels

• Total number of branch points

For the tree structure-based quantification, the statistics are grouped in terms of sub-tree label. Each sub-tree contains following measurement:

• Tree label

• Total length of the tree

• Length of longest vessel from the root

• Max branch layer

• Mean branch layer

• Number of branch points

• Number of vessel segments

• Number of extreme vessels

The length and mean radius of the individual vessel segment are output in a list. Two images produced to present the vascular tree intuitively are shown in Fig. 5. One is vessel trees, in which the vessel branches are of grey values according to its sub-tree label. This image presents different trees in whole. The other one is the branch level image. It shows the hierarchy of the vasculature with different grey values in terms of the branch level.

Fig. 5Full-screen View

Visualization of vascular tree structure by output images. (a) Vessel trees; (b) Branch levels

4 Quantification results

A series of mouse liver images are used to evaluate the proposed quantification method of tree analysis. The images are grouped according to different stage of hepatic fibrosis. There are 8 set of images for mouse livers in total, including three normal livers, three sets at mild stage of hepatic fibrosis and two sets at advanced stage. In order to achieve statistical significance, we truncated the images of relatively intact samples to obtain more images, which simulate the incomplete or non-intact samples. First, the longest axis is found inside each sample’s convex volume. The orthogonal plane at one third of the axis is used to remove some part of the liver samples. One non-intact sample is generated from samples called “NA”, “NB”, “3A”, and “mouse2”, respectively. In advanced group, the sample “mouse3” generated two non-intact samples, and the existing non-intact sample “mouse5” removes one fourth of its tissue to generate worse non-intact samples. Therefore, there are 5 samples in each fibrotic group, which are enough to demonstrate the statistical significance to some extent. They are quantified according to the tree analysis workflows described in the section above. As an example, the images for the tree analysis of one normal mouse liver (Sample “NB”) are given in Fig. 6.

Fig. 6Full-screen View

Quantification for a normal mouse liver: (a) raw image, (b) the vasculature, (c) the vascular skeleton, (d) vascular trees of the liver, and (e) the branch levels of all vessel segments.

Fig. 6(b) is the segmentation result of the vasculature and (c) is its skeleton. They are the critical intermediate step of our quantification method. Figs.(d) and (e) are the final output images of proposed tree analysis, which intuitive display of the vasculatures in hierarchy. Corresponding statistics of the liver “NB” are listed in Table 2.

The data in Table 2 tend to be too abstract to study the hepatic fibrosis. To evaluate the level of angiogenesis for the investigation of hepatic fibrosis or other vasculature related diseases, the density of vascular segments, terminal vessel branches and bifurcations are usually used^{[21, 22]}. These parameters are employed in our analysis to all the three sample groups. These densities, described by the volume used for density evaluation, are computed by integrating the convex hull of the vasculature. The correspondent results are shown in Fig. 7.

Fig.7Full-screen View

The common use indices of evaluating angiogenesis for all specimens. (a) micro-vessel density (MVD), (b) ratio of junction point to volume (J/V), (c) ratio of extremities to volume, (d) fractional volume.

According to the density values given in, only the average value of segment density show obvious increasing trend (106.4, 121.0 and 142.4 for samples at normal, mild and advanced levels of fibrosis, respectively), while those for the extremities and branch points do not follow this regularity of fibrotic progression. Abnormal values could be found in the charts, such as the density for sample “NA” and “3A”. In addition, the large standard deviations in show obvious fluctuation of density value among samples in the same group, implying that a normalized parameter should be defined to evaluate the level of fibrosis more accurately.

The angiogenesis is a prominent characteristics of hepatic fibrosis ^[22]. It is a process that vasculature rapidly proliferates. The vessel segments outgrow from existing vessels^[23], and the newly generated vessels normally occupy higher levels of the vascular tree. Thus the mean level of a vascular tree increases as the angiogenesis develops. The indicator for branch level is normalized so that it can be used in statistics. To evaluate the level of angiogenesis, a new parameter called High-level Vascular Ratio is defined as HVR = Average level / Maximum level, where the average level denotes the mean level of a single vascular tree and the maximum level is the highest level of the same tree.

The HVR is extracted from the level quantification of the vascular tree. The parameter is defined based on a typical kind of vascular tree so that criterion for comparison among samples can be standardized. As well known, most of liver blood flow is through the hepatic portal vein, vasculature of which contains the thickest vessel segment in the liver tissues. The hepatic portal vein is characterized in terms of the vessel radius and the branch level in vasculature trees. Therefore, the hepatic portal vein of each specimen is identified from the quantitative parameters, including mean radius of vessel segment and branch level. For a compromise between accuracy and efficiency, branch level of 5 for the vessel segments is usually large enough to identify the hepatic portal vein. Fig. 8 shows the HVR value of each hepatic portal vein in the eight samples. The HVR correctly reflects the trend of angiogenesis when the degree of hepatic fibrosis aggravates, the larger the HVR value is, the higher the angiogenesis level. For samples in the same group in Fig. 8, the fluctuation of HVR value is greatly reduced compared to the conventional density parameters in Fig. 7. This means that the proposed HVR parameter can be used to evaluate the degree of fibrosis accurately. Among all the samples, most of the vasculatures used for the quantitative analysis are not intact, implying that the HVR parameter is robust even if a specimen misses a portion of vasculature.

Fig. 8Full-screen View

HVR of each hepatic portal vein inside the liver samples.

5 Discussion

The segment or bifurcation density is usually adopted to evaluate the degree of hepatic fibrosis in some cases. However, These parameters do not always work and may lead to incorrect assessment results, in cases that the three dimensional vasculature of the samples is not imperfect due to factors in sample preparation or image collection. The experimental results demonstrate that the proposed parameter HVR is robust and reliable in the evaluation of the angiogenesis, compared to those commonly used density measurements. In the supplementary material we provided, the videos for three dimensional CT images of all the liver specimens are given, some of which are not intact liver samples with missing tissues in the stages of the sample preparation or imaging. As the vessel segments are not uniformly distributed inside a liver, the density measurements relying on the segment number and the vascular convex volume may vary significantly in different liver regions. With the proposed measurement in this paper, the tree analysis of hepatic portal vein for each specimen shows that the new measurement can obtain a robust and reliable trend of the angiogenesis development which accurately reflects the degree of hepatic fibrosis. Therefore, we can conclude that the proposed method is practical for the evaluation of the angiogenesis, despite of all kinds of image data imperfections. The proposed method can also be extended to investigate other diseases in which vasculature change is a key indicator, such as tumor in the period of formation and growth, and cardiovascular diseases.