A robust method for high-precision qu­antification of the complex three-dimensional vasculatures acquired by X-ray microtomography

The quantification of micro-vasculatures is important for the analysis of angiogenesis on which the detection of tumor growth or hepatic fibrosis depends. Synchrotron-based X-ray computed micro-tomography (SR-µCT) allows rapid acquisition of micro-vasculature images at micrometer-scale spatial resolution. Through skeletonization, the statistical features of the micro-vasculature can be extracted from the skeleton of the micro-vasculatures. Thinning is a widely used algorithm to produce the vascular skeleton in medical research. Existing three-dimensional thinning methods normally emphasize the preservation of topological structure rather than geometrical features in generating the skeleton of a volumetric object. This results in three problems and limits the accuracy of the quantitative results related to the geometrical structure of the vasculature. The problems include the excessively shortened length of elongated objects, eliminated branches of blood vessel tree structure, and numerous noisy spurious branches. The inaccuracy of the skeleton directly introduces errors in the quantitative analysis, especially on the parameters concerning the vascular length and the counts of vessel segments and branching points. In this paper, a robust method using a consolidated end-point constraint for thinning, which generates geometry-preserving skeletons in addition to maintaining the topology of the vasculature, is presented. The improved skeleton can be used to produce more accurate quantitative results. Experimental results from high-resolution SR-µCT images show that the end-point constraint produced by the proposed method can significantly improve the accuracy of the skeleton obtained using the existing ITK three-dimensional thinning filter. The produced skeleton has laid the groundwork for accurate quantification of the angiogenesis. This is critical for the early detection of tumors and assessing anti-angiogenesis treatments.

1. Introduction

The analysis of micro-vasculature is important for the understanding, diagnosis and treatment of diseases involving vasculopathy, such as cardiovascular disease, cerebral infarction and tumors. The development of micro-computed tomography (micro-CT), particularly synchrotron-radiation-based micro-tomography (SR-µCT) (Chen et al., 2014), provides means to examine the minute detail of vasculature within an organ (Wan et al., 2000; Heinzer et al., 2008' Folarin et al., 2010; Vasquez et al., 2011; Lang et al., 2012). X-ray propagation-based phase-contrast imaging allows researchers to inspect the microstructures of vasculatures without contrast agents. Therefore, many micro-vasculatures wherein the particles of the contrast agent are too large to penetrate can be observed by such non-invasive imaging. Images with elaborate vessel structures introduce challenges to the structural analysis and the extraction of quantitative information. Generally, skeletonization is an effective way to analyze the vasculature because the structural parameters are readily conveyed in terms of this `compact' representation. Thinning is a frequently used skeletonization method for the analysis of vasculature because the generated skeleton possesses various properties (Cornea et al., 2007). Some open source software packages and libraries provide implementations of thinning algorithms, including ImageJ (Arganda-Carreras, 2008), OpenCV (Bradski, 2008) and the Insight Segmentation and Registration Toolkit (ITK) (Homann, 2007). However, most of the thinning algorithms focus on representing the topological structure rather than the geometrical features of a volumetric object. Consequently, the generated skeletons have some limitations which would significantly affect the accuracy of the quantitative results of the vasculature analysis (Vallotton et al., 2007; Marks et al., 2013) concerning the geometrical features, such as the number of vascular segments, number of branching points, length and number of terminal vessel branches. These are normally considered as critical parameters to evaluate the angiogenesis in the detection of tumors, hepatic fibrosis (Paternostro et al., 2010; Bocca et al., 2015) or other related diseases.

A side effect of the boundary peeling process of a typical thinning method is the shortening of the length of an elongated object as shown in Fig. 1(a). For a high-resolution image, these shortened lengths can accumulate to become a significant error which is critical for the detection of minor changes of the vasculature. Moreover, some thinning algorithms can result in the complete elimination of some vessel branches in a vascular tree. Fig. 1(b) demonstrates that the frequently used ITK three-dimensional (3D) thinning method with weak end-point constraint fails to preserve the geometry by eliminating branches, although it does not change the overall topology of the original object. Preservation of topology can be stated simply as follows: two objects have the same topology if they have the same number of connected components, tunnels and cavities (Kong & Rosenfeld, 1989; Cornea et al., 2007). The preservation of geometry is another critical requirement for the quantification of vasculature besides topology. For example, an object such as the character `b' cannot be transformed into a skeleton similar to the character `o' (Palágyi & Kuba, 1999). Another problem produced by the existing thinning algorithms is that the resulting skeleton may contain many spurious branches as shown in Fig. 1(a) because the skeleton is intrinsically sensitive to small changes on the object's boundary. These spurious branches are the erroneous representation of the original object.

Figure 1 (a) Skeleton showing the reduced length of an elongated structure and a spurious branch. (b) Skeleton produced from the ITK 3D thinning filter fails to preserve the geometrical structure and has missing branches.

This article investigates the addition of the requirement for preserving geometrical information as accurate as possible for the thinning algorithm in order to obtain precise quantitative information of the vasculature. In §2 we review the existing thinning methods, focusing on the issue of the geometrical properties of the obtained skeleton. Then, in §3, in order to overcome the limitations of the existing thinning algorithms to generate a desired skeleton, we propose a new solution to generate a robust end-point constraint in the thinning process. The method is implemented by detecting the end-points of terminal branches in a vessel tree before thinning starts. We have applied the new method in a thinning algorithm to extract a skeleton from several synthetic and real images to evaluate the results in §4. The evaluation results show that the proposed method can produce more accurate quantitative results from the geometry-preserving skeletons with preserved length and shape information and much less spurious branches. §5 discusses why some vascular tree branches are eliminated during thinning if no end-point constraint or a weak end-point constraint is applied. This is followed by conclusions in §6.

2. Geometrical properties of thinning algorithms

As a dimensionality-reduced representation, the skeleton provides an efficient way to analyze 3D volumetric objects. The first and the most important feature of a skeleton is the preservation of topology for the object. Many thinning algorithms have been developed in the past two decades to generate skeletons which preserve topology rather than geometry. A comprehensive review of digital topology can be found by Kong & Rosenfeld (1989), and more information on topology preserving based thinning has been detailed (Lee et al., 1994; Saito & Toriwaki, 1995; Ma & Sonka, 1996; Saha et al., 1997; Palágyi & Kuba, 1999; Palágyi et al., 2001; Xie et al., 2003; Wang & Basu, 2007; Palágyi et al., 2012).

2.1. Geometry features

Topology preservation is necessary but not sufficient for a skeleton to represent the structure of a vascular tree, which is the prominent structure of vasculatures. The preservation of geometrical features is crucial for the quantitative analysis of the vasculatures. As the prominent information about the geometric shape of the object is concentrated at the boundary of an object (Attneave, 1959), it will be helpful to add some boundary constraints to preserve the geometrical features in the thinning process. In fact, there are reports in the literature that end-point constraints are used in the thinning of linear structures to maintain their geometrical properties (Lee et al., 1994; Palágyi & Kuba, 1999; Xie et al., 2003; Wang & Basu, 2007; Palágyi et al., 2012). In image analysis, pixel connectivity is the way in which pixels (i.e. voxels in a 3D image) relate to their neighbors. Twenty-six connected voxels are neighbors to every voxel that touches one of their faces, edges or corners. As the curve skeleton for quantitative analysis is one pixel wide, the commonly used end-point constraint is defined below:

Definition 1. An object point P is a curve end-point if the set of its 26 neighbors N₂₆(P) contains exactly one object point.

The geometrical features of a vascular structure (Vallotton et al., 2007), such as shape and length (Wang et al., 2010), are directly related to the counting of the vascular segments and the length measurement of the vascular vessels. Some work has been reported using the end-point constraints to preserve the geometry of an object (Lee et al., 1994); however, the end-point constraint in the literature is not suitable for the generation of a desired skeleton that maintains the adequate geometrical features of vasculature. Cornea et al. (2007) have pointed out that an additional condition needs to be used in cases where removing end-points results in eliminating branches; however, they did not provide a solution.

Saha et al. (1997) described a shape-preserving algorithm, to generate a medial surface for object recognition and description. Their algorithm implemented shape preservation by identifying specific `shape points' which indicate the corners of the object, and keeping them non-removable during thinning. The algorithm can be used to maintain the crucial structure in addition to the main topology of the object.

The skeletons generated by traditional thinning algorithms are often shorter than the actual length of an object (Cornea et al., 2007; Choi & Seong, 2008). This reduction of the length is much more remarkable in high-resolution images, which would introduce considerable errors to the length measurement. In order to mitigate the length error, Choi & Seong (2008) proposed a length-preserving thinning method for extracting river skeletons from two-dimensional (2D) land cover image data. This algorithm maintains the length by keeping the boundary of the rivers' ends invariant during thinning and using a following-up skeleton growth procedure to compensate the over-reduced length. Their method is designed for 2D images, and the out-growth procedure following the skeleton extraction never fully resolves the issue associated with spurious branches and the eliminated branches in vascular skeletonization.

Pruning (Shaked & Bruckstein, 1998) is a common solution for removing the spurious branches. It is implemented by removing the short terminal branches whose lengths are smaller than a threshold. However, there is no strict criterion to evaluate a spurious branch in a skeleton. Besides pruning, Palágyi (2014) proposed a 3D curve-thinning algorithm based on isthmuses, which results in a skeleton with less spurious branches of the 3D vasculature. All these post-processing workflows and algorithms focus on removing spurious branches but do not consider the issues of branch elimination and length shortening.

In some medical applications, the detection of subtle changes of blood vessels is critical. For example, when using micro-tomography to analyze the angiogenesis during tumor formation and growth, or the evaluation of the tumor angiogenesis inhibitors, the detection and quantification of the minor variations of the micro-vasculature is important for the early detection of the tumor formation and growth. The blood vessels never actually end but the inspection devices lead to terminal blood vessels of the vasculature. To some degree, the density of the terminal vessel branch is often a useful parameter for evaluating the change of vasculature caused by angiogenesis (Bocca et al., 2015). There is no report so far that a thinning algorithm can address all three problems mentioned in the abstract, namely, the excessively shortened length of elongated objects, eliminated branches of tree structures and numerous spurious branches.

2.2. End-point constraint

End-points are critical in the extraction of a skeleton. They correspond to the terminal points of the vascular branches. They are located on the boundary of the volumetric object and can be potentially used as constraints to preserve geometrical features such as the lengths of the terminal vessels. The end-point of an elongated object is the voxel which has exactly one object neighbor in its 26 neighbors N₂₆(P). However, properly and promptly identifying the end-points of the volumetric vasculature before skeletonization is not an easy task.

Zhou & Toga (1999) proposed a voxel-coding method based on distance transformation to detect the end-points of the vasculature from volumetric objects. Methods given by Shahrokni et al. (2001) and Maddah et al. (2003) also demonstrated such a voxel-coding approach used to extract the end-points. They all used these pre-detected end-points to extract a skeleton by means of a path planning process. However, their methods were only tested using a simple vessel tree with limited number of branches. When the structure of the vascular trees becomes complex and the number of end-points increases enormously, their algorithms will soon encounter problems in identifying the end-points effectively and efficiently.

In this paper, we propose a new method to detect the end-points of vascular and microvascular branches in high-resolution synchrotron images. The end-points are identified using a voxel-coding and a lookup point list. Then the detected end-points that are located at the boundary of the object are applied as the constraints in the extraction of the skeleton. The extracted skeleton can preserve the geometry properties of the vasculature and therefore leads to much more accurate quantification of the vasculature.

3. Proposed method

Thinning is used to extract the skeleton of the vasculature. Before the thinning procedure, all of the branch end-points are detected using our proposed method. Fig. 2(a) illustrates a flow chart of the proposed method. Each step in the flow chart is explained below.

Figure 2 (a) Flow chart of the proposed algorithm for extracting the end-points. (b) Extracted end-points in a 9 × 9 × 9 cube. (c) End-points after refinement.

3.1. Input image

The 3D input image is in binary format, which can be depicted as (Z ³,m,n,B), where Z ³ is the 3D image space; (m, n) means the connected relationship of voxels, which has m-adjacency for the object voxels and n-adjacency for the background voxels in the digital picture; and B is the set of object voxels. We use the frequently adopted relationship (m,n) = (26, 6) in this paper. The binary image is the segmentation result of the blood vessel network from a high-resolution CT image, and some pre-processing algorithms have been applied before identifying end-points, such as filling up small holes inside a blood vessel and joining the disconnected vessels. Therefore the input 3D vascular image is one connected component without any tunnel or cavity.

3.2. Detecting root area

The reference point is fundamental in the distance transform. It is also the start point for tracing the structure of a vascular tree. An arbitrary point near the root of the vasculature could be a reasonable reference point to start with. Instead of using a single arbitrary point (Zhou & Toga, 1999), we use the entire root area within a slice of the image as a reference to guarantee the shortest path from each end-point to the unique reference. We scan each slice from the first frame to the last one along the six directions of the axes of the 3D image space Z ³(X,Y,Z,-X,-Y,-Z). In the first frame, if any object voxel is detected, the largest and central 2D connected regions are selected as the candidate root area. Then we calculate the largest inscribed sphere inside the vasculature to obtain the approximate location of the root. One of the six candidate root areas, which is the closest to the inscribed sphere, is considered as the reference root for voxel-coding. A manually specified point can also be accepted as the reference if this is preferred.

3.3. Distance transformation

We use a distance map to identify the end-points of the vascular branches. The distance map is an image showing the distance values of all the object points to a set of voxels. In order to generate the distance map for identifying end-points, we calculate the distance of each object voxel representing the nearest connected path from the voxel to the reference root area.

The voxel-coding scheme we employed is similar to the one proposed by Zhou & Toga (1999), but with the following improvements: (i) we use the automatically detected root area rather than the manually designated point as the reference; and (ii) we employ a point list procedure to automatically refine the end-points as there are false ones extracted only by voxel-coding in a high-resolution image. We will explain the refinement processing in §3.5.

The voxel-coding is a sequential procedure of accessing points and calculating distance values. When a point is accessed, a distance value is assigned to all of its 26 neighbors based on the distance value of the current accessing point by the Chamfer distance transform (CDT) 〈3, 4, 5〉 metric (Zhou & Toga, 1999). We initialize all the object voxels' distances with a large value, and the reference area is initialized with 1. For instance, when we visit a voxel P (with a distance value of V), we assign distance values to its 26 neighbors. Its face-connected neighbors, edge-connected neighbors and vertex-connected neighbors are assigned with the distance values (V + 3), (V + 4) and (V + 5), respectively, if the assigned value is smaller than the current distance value of the assigned point. In order to ensure that each point we access contains a distance value before assigning values to its neighbors, we use the mechanism of region growing to traverse all object points (Hauke & Groher, 2007). The root area is considered as the seeds for region growing. This procedure continues until all object points are assigned a distance value indicating how far away it is from the reference area.

3.4. End-points extraction

On the basis of the Chamfer distance map generated in the previous step, the end-points are identified as the local maxima. Nevertheless, the lump or distortion on the surface of the blood vessels is frequently seen in vascular images, particularly in the high-resolution 3D images; therefore, some points near the lump or protuberance of the coarse vessel surface could be detected as end-points. This could potentially cause the spurious branches in the skeleton. The spurious branches are not desirable because they could lead to false structure in the skeleton and introduce errors in the quantification of the vasculature.

There are several factors that contribute to the coarse boundary of the blood vessel, and subsequently result in spurious branches in the skeleton of the blood vessels. The first is that the blood vessels are not in the ideal tubular shape because of certain vasculopathy. The second concerns some artificial intervention at the imaging stage. For instance, the SR-µCT requires that samples are to be dehydrated. The dehydration often introduces distortion on the vascular wall. Another factor causing the coarse boundary is the use of threshold-based segmentation when producing the binary image. Using a single threshold for the entire image does not guarantee preservation of the continuous and smooth boundary for the segmented blood vessel.

In identifying the end-points, a small region size is used to find the local maximum. False end-points may be introduced because of the lump or protuberance on the lateral surface of vessels, such as the point `A' illustrated in Fig. 2(b). On the other hand, a large region size may result in higher computational cost. Therefore, the size of the region used for the end-points detection shall be decided carefully in order to avoid false end-points while not significantly increasing the computation time.

In our study, a 9×9×9 cube is used to find local maximum values and eliminate false end-points. As the computing cost of traversing a cubic region is O(n³) in the 3D space, we devise a strategy which splits the process of finding local maximum values into two steps. Firstly, we traverse the entire image and find the local maxima using a 3×3×3 cube; then check the end-points detected in the previous step by using a 9×9×9 cube. In comparison with directly traversing the 9×9×9 cube for each object point, this strategy significantly reduces the traversal number for accessing image points.

3.5. End-point refinement

There is still another issue that needs to be resolved in determining the final end-points. In the high-resolution images we used, the terminal of a branch is normally a flat surface rather than a sharp point. Therefore we find more than one end-point for some branches, such as the points of type `B' illustrated in Fig. 2(b). Increasing the size of the neighborhood region to find local maximum values can potentially reduce the number of end-points detected from the terminal surface of a vessel branch; however, it is not efficient because the computation increases tremendously when traversing a larger cubic region. We propose to use a lookup point list to refine the end-points because we need only one end-point for each terminal vessel branch.

Duplicated end-points are normally close to each other, likely located on the same terminal surface of a vessel branch. The connecting line of these points is limited within a length threshold and contains no background voxels. These two conditions can be used to remove duplicated end-points for a vessel branch. This refinement process can be depicted as follows:

Assume P(x₁,y₁,z₁) and Q(x₂,y₂,z₂) are the end-points extracted using the method described in §3.4, the algorithms used for the refinement can be formulated as

At first, we use the city block distance (Ren et al., 1998) as the distance metric for P(x₁,y₁,z₁) and Q(x₂,y₂,z₂) in the refinement. Based on the definition of the relationship for two distinct points, we add all the end-points in a lookup list and rank them in an ascending order.

As the duplicated end-points are normally located on the same terminal surface of a vessel branch, the diameter R of a maximal inscribed sphere inside the vasculature can be used as a threshold of distance to check two duplicated end-points of the same branch, i.e. < R. As mentioned earlier, if two duplicated end-points are from the same branch, their connected line contains no background-point. After finding duplicated end-points from a branch, we consider the middle point of these duplicated end-points as the only end-point of the branch. Some branches may contain more than two end-points; they are processed using the same principle. Finding the middle point for each duplicated end-point pair in the list can be processed in parallel. Therefore the processing of the lookup point list can be performed efficiently.

The refinement ensures that all the end-points of the terminal branches of the vascular tree are consolidated and ready to be used as the constraints for the following thinning process.

3.6. Thinning

The extracted end-points are then used as the end-point constraint for the thinning algorithm (Lee et al., 1994). The new end-point constraint can ensure that the resulting skeleton not only maintains the topology but also preserves accurate geometry properties such as shape and length. As each end-point corresponds to only one terminal branch of the vasculature, the resulting skeleton is clean and has no spurious branches.

4. Experimental results

The proposed method has been tested using both synthetic and actual images. The method is implemented in C++ and parallelized using OpenMP (The OpenMP API specification for parallel programming; https://openmp.org/). In order to verify the proposed method, we employed the proposed end-point constraint to the 3D binary image thinning filter in the ITK library (Homann, 2007), which is widely used in skeletonization (Lee et al., 1994). The proposed method could be used by other thinning algorithms as well. The quantitative features of both synthetic and real images are extracted to evaluate and compare the skeleton generated with and without the proposed end-point constraint. The vascular image used in this experiment is constructed from high-resolution micro-CT images of a typical vascular tree of a mouse liver.

4.1. Experiment with synthetic images

The synthetic image we created is composed of five connected cuboids (20 × 20 × 198 voxels) as shown in Fig. 3(a). Fig. 3(e) shows the skeleton extracted from a synthetic image using the ITK 3D binary image thinning filter with the proposed end-point constraint; Fig. 3(b) illustrates the skeleton produced using the original ITK 3D binary image thinning filter. As demonstrated in Fig. 3, skeleton (e) maintained the shape of the original object while skeleton (b) with end-point constraint Definition 1 failed to preserve the shape of the object because of missing branches (segments 10, 11, 12, 13 and 14).

Figure 3 (a) Binary image. (b) Skeleton directly generated using ITK 3D thinning filter with weak end-point constraint. (c) Extracted end-points using a 9 × 9 × 9 cube. (d) Refined end-points. (e) Skeleton generated by our algorithm.

The length of the skeleton is also assessed to evaluate the geometrical preservation for the volumetric object. The measurements of the individual branches of the two skeletons are shown in Fig. 4.

Figure 4 Measurements of the individual branch of the two skeletons.

The lengths of the segments 6–9 are almost identical in both skeletons. However, the terminal branches 1–5 in the skeleton showed in Fig. 3(b), which are generated without applying our end-point constraints, are obviously shortened; and the branches 10–14 are completely eliminated. The accumulated errors caused by the inaccurate skeleton will result in false measurement for the quantitative analysis of the vasculature, especially for the analysis of tumors or fibrosis livers, of which the number of terminal vessel branches is a critical parameter in the evaluation of angiogenesis.

4.2. Experiment with real image data

Hepatic vasculature images of mice are introduced to assess the accuracy of the proposed method. Phased-contrast SR-µCT was employed to capture images of the specimen at beamline BL13W1 (SSRF, Shanghai, China). The energy of the X-ray beam was set at 15 keV. The image voxel size was 3.7 µm × 3.7 µm × 3.7 µm. The detector was positioned at 70 cm downstream of the sample. The CT scanned 720 projections at angles evenly distributed between 0° and 180°. Exposure time per projection was 2 s. The total time of X-ray exposure was approximately 35 min. After image acquisition, the dataset was reconstructed using the software PITRE (Chen et al., 2012). Phase retrieval was carried out for each projection image, and then the stack image was reconstructed using the filtered back-projection algorithm. The stack image was denoised by a 3D Gaussian filter and a ring artifacts removal algorithm (Zhang et al., 2012). Then the liver was isolated from the image background so that the threshold segmentation could be used to extract more accurate vasculature within the liver region. As a result, the binary vasculature image was prepared for skeletonization. For the purpose of efficiency, thinning skeletonization was accelerated by a parallel algorithm (Hai et al., 2015), by which the computational time decreased by an order of magnitude. Finally, the quantitative parameters of the vasculature were extracted to evaluate the performance of the proposed method.

Fig. 5(a) shows the entire vascular network of the rat's hepatic portal vein; Figs. 5(b) and 5(d) show the skeletons extracted using the ITK thinning filter and the same thinning filter with our proposed end-point constraint, respectively. They both represent the overall shape of the vessel network. However, if we look into the micro-structures, we can identify a lot of differences between the two skeletons. For example, the skeleton shown in Fig. 5(d) looks much cleaner than the one shown in Fig. 5(b). There are many spurious branches in Fig. 5(b) compared with the skeleton shown in Fig. 5(d). This is evident inside the marked circles.

Figure 5 (a) Hepatic vein of a rat. (b) Skeleton generated using the ITK 3D thinning filter. (c) End-points extracted from the terminal vascular branches. (d) Skeleton generated based on our method.

In order to examine the vessel segments in detail, we truncate one piece of the skeleton to zoom in and visualize the difference between the two skeletons produced using different end-point constraints. As showed in Fig. 6, the skeleton in Fig. 6(b) is produced using the proposed end-point constraints, and it preserves the length of the vascular branches. On the contrary, the skeleton shown in Fig. 6(a) apparently shortens the vascular terminal branches. Figs. 6(d) and 6(e) show a comparison of the skeletons generated using the ITK thinning filter and our method, respectively. They are extracted using the image shown in Fig. 6(c) as the input. It clearly shows that our proposed method can be used to generate a much cleaner skeleton as shown in Fig. 6(e); however, the original ITK 3D thinning filter failed to represent the correct structure of the vessel tree. The skeleton shown in Fig. 6(d) contains many spurious branches. These figures show that the proposed method can significantly improve the skeleton in terms of preserving the length of the vascular branches and producing a clean skeleton with much less spurious branches.

Figure 6 (a) Skeleton showing shortened branches. (b) Skeleton with preserved length using our method. (c) Sub-tree of the rat's hepatic vessel network. (d) Skeleton with many spurious branches generated using the ITK thinning filter. (e) Skeleton generated based on our proposed end-point constraint.

To verify our method, we have also used manual measurement results as a standard to evaluate the skeletons generated by thinning algorithms with different end-point constraints. This type of manual verification is frequently used and considered as a reliable approach to assess the results produced from an automated method (Weidner, 1995; Haisan et al., 2013). As it is impossible to observe all the vessel branches by manual counting, in the experiment, we truncated 20 small vascular sub-trees from different regions of the hepatic vasculature for the evaluation of the performance of the proposed method. Fig. 7 shows skeletons of the two sub-trees generated using our method and the original ITK thinning filter. Five observers experienced in medical image analysis were asked to manually count the number of terminal branches for each sub-tree. To avoid human error which may be introduced by bias, we have carried out cross-checking and considered that a sub-tree has a valid terminal branch number only if three or more people arrive at the same count. We had 18 valid sub-trees out of 20 as shown in Fig. 8.

Figure 7 (a) Sub-tree 6. (b) Skeleton produced using our method. (c) Skeleton produced using the original ITK filter. (d) Sub-tree 18. (e) Skeleton produced using our method. (f) Skeleton produced using original ITK filter.

Figure 8 Comparison of the manual counts of the branches of the sub-trees and their skeletons.

We use a t-test to evaluate the difference between an automatically generated skeleton and the manual counting one in Fig. 8. At the level of 5% significance, the P value of Group 2 and Group 1 in Fig. 8 is smaller than 0.05, which indicates the obvious difference of the two group of samples. However, the P values of Group 3 and Group 1 are larger than 0.05, which shows no difference of the terminal branch number. The experimental results demonstrate that the skeleton produced using the proposed method is more accurate than that generated using the original ITK thinning method with the weak end-point constraint.

The noise in a vascular skeleton is presented as the spurious branches which result in an error count of the vessel branches. The mean error count between Group 2 and Group 1 is 2.61 while the value between Group 3 and Group 1 is 0.33. This means that the skeleton produced by the proposed method represents the vasculature well and is close to the human recognized skeleton. The results listed in Fig. 8 also show that more noise occurred in the skeleton generated by ITK 3D thinning with weak end-point constraint than that extracted using the proposed method.

4.3. Quantitative analysis of hepatic fibrosis vasculatures

Angiogenesis is one of the characteristics of hepatic fibrosis. The parameter, density of terminal vessel branches, is often used to evaluate the degree of hepatic fibrosis. In principle, this parameter has an increase trend when the fibrosis becomes more serious. Fig. 9 demonstrates analysis of the normal liver, mild stage and advanced stage of hepatic fibrosis. The parameters are extracted from the two skeletons generated using different end-point constraints, i.e. the proposed end-point constraint and the traditional one (Definition 1).

Figure 9 Quantitative result of hepatic fibrosis samples.

In the chart at the bottom right of Fig. 9, the blue columns show a correct trend of hepatic fibrosis when using the proposed method while the red ones do not demonstrate a reasonable trend. More liver samples are being collected now in order to present the statistical analysis results using our method in the future work.

5. Discussion

5.1. Missing branches in the skeleton generated by the ITK 3D thinning method

Thinning is widely used for quantitative assessment in medical research. It is unacceptable for the vasculature to lose branches after skeletonization. We have analyzed the cause of the missing branches in Figs. 1(b) and 3(b), and this turns out to be the problem with the weak end-point constraint (Definition 1) used in the ITK 3D thinning method. The missing branches are removed at the stage of thinning when a 2 voxel-wide object is peeled into a 1 voxel-wide skeleton. At this stage, none of these points is considered as an end-point by Definition 1. All points of the branch are added to the list to go through the sequential (by the order of the raster scan in Fig. 10) re-checking for deletion in thinning (Lee et al., 1994). Therefore, all of these points are sequentially deleted until only one voxel is left. In contrast, the proposed rigorous end-point constraint prevents this from happening by ensuring that there is always a connected path from the fixed end-point `1' to the last point `12'.

Figure 10 Procedure of eliminating a branch by ITK 3D thinning.

In fact, the skeleton Fig. 3(b) extracted by the ITK 3D thinning with the weak end-point constraint is a figure that contains one node point with five branches. The branches 1 and 6 are considered as one single branch since there is no voxel containing more than two voxels in its neighbors along them. This shape of the skeleton is a wrong representation of the original object so that the quantitative result from the skeleton is unreliable. The proposed method can be used to produce accurate skeletons by preserving the geometrical features of the vasculature without missing branches and apparently reducing spurious branches. This will lead to more accurate quantitative analysis.

5.2. Comparison with other image analysis workflows

More refined thinning algorithms have been published in the last two decades. Various useful post-processing workflows have been proposed to refine the skeleton. To compensate for the over-shortened length of the skeleton, the consideration of a maximal inscribed sphere at the end-point or the outgrowth from the end-point along the medial axis (Choi & Seong, 2008) should give an accurate result of the branch length. The compensation length needs a precondition of removing the false spurious branches of the skeleton and avoiding missing branches. Other workflows, such as pruning, are introduced, which is commonly used to carry out this based on a defined threshold of branch length. Any branches whose lengths are smaller than the threshold will be definitely removed, including the bifurcated branches close to the vessel terminal. This irregular structure is commonly seen in the analysis of angiogenesis. Therefore, the traditional pruning may accidently remove bifurcated branches beside the eliminated branches discussed in the previous section, which may introduce further errors and produce an incorrect skeleton of the original object. The experiments show that the proposed end-point constraint can be used in thinning to produce a robust skeleton effectively, not only preserving all branches and their lengths but also eliminating the spurious branches. Finally, the precise quantitative analysis of vasculature is achieved using the geometry-preserving skeleton produced by the proposed method.

5.3. Stubs in the skeleton produced by our method

Although the proposed end-point constraint can be used in the thinning to generate a much cleaner skeleton than the commonly used one (Definition 1), there is still a slight difference in terminal branch numbers between Group 1 and Group 3 as shown in Fig. 8. This is caused by the presence of some tiny stubs which do not correspond to any vessel branches; refer to the marked regions in Fig. 11. We have identified that these stubs are a false representation of the vasculature and they occur at the location where a coarse blood vessel bends sharply. Although these stubs do not frequently appear in our hepatic vasculature, they should be removed in order to obtain a more precise skeleton for the subsequent quantitative analysis.

Figure 11 (a) Sub-branch 19. (b) Skeleton with fake terminal branches caused by sharply bent vessels.

6. Conclusions

As the widely used ITK 3D thinning fails to generate a skeleton that preserves geometry, we have developed a novel approach for producing a robust end-point constraint which can be used in the skeletonization of vasculature. In order to identify the end-points, we use a distance map generated from a reference root area with the Chamfer distance metric 〈3, 4, 5〉 at first. Then the local maximal values in the distance map are examined to determine the end-points. The end-points are further consolidated by removing duplicated and false ones. Finally the consolidated end-points are employed as the end-point constraint in the thinning algorithm. The experimental results show that the ITK 3D thinning can be significantly improved by using the end-points identified using our proposed method. The improved thinning filter can be used to extract more accurate skeletons from high-resolution vasculature images.

By using the proposed method, the extracted skeleton of the vasculature can not only precisely reflect the topology of the vasculature but also preserve its geometry properties. This is especially critical for quantifying micro-vasculature in detecting subtle changes of the micro-vasculature with very fine structures for the early detection of tumor formation, growth and hepatic fibrosis, and for the assessment of angiogenesis treatment. The proposed method can also be potentially used in the skeletonization of any elongated 3D objects besides the 3D vasculature.