Abstract
Vascular segmentation plays an important role in medical image analysis. A novel technique for the automatic extraction of vascular trees from 2D medical images is presented, which combines Hessian-based multiscale filtering and a modified level set method. In the proposed algorithm, the morphological top-hat transformation is firstly adopted to attenuate background. Then Hessian-based multiscale filtering is used to enhance vascular structures by combining Hessian matrix with Gaussian convolution to tune the filtering response to the specific scales. Because Gaussian convolution tends to blur vessel boundaries, which makes scale selection inaccurate, an improved level set method is finally proposed to extract vascular structures by introducing an external constrained term related to the standard deviation of Gaussian function into the traditional level set. Our approach was tested on synthetic images with vascular-like structures and 2D slices extracted from real 3D abdomen magnetic resonance angiography (MRA) images along the coronal plane. The segmentation rates for synthetic images are above 95%. The results for MRA images demonstrate that the proposed method can extract most of the vascular structures successfully and accurately in visualization. Therefore, the proposed method is effective for the vascular tree extraction in medical images.
1. Introduction
Accurate segmentation and quantification of vascular structures in medical images is a critical task for clinical practices such as computer-aided diagnosis, treatment, surgical planning, and navigation. However, it is highly challenging to extract vascular structures in 2D and 3D medical images. The reasons lie in two aspects. On one hand, some vascular structures involve numerous vascular branches and complex patterns [1]. On the other hand, noise, variations in intensities, and low image contrast pose difficulties in vascular tree extraction [2].
Various extraction techniques have been proposed for vascular tree segmentation, that is, pattern recognition techniques, model-based approaches, mathematical morphology, multiscale filtering approaches, vessel tracking, and matched filtering (see Kirbas and Quek [3] and Lesage et al. [4] for comprehensive reviews). Almost all the vascular extraction techniques take advantage of the characteristics of tubular-like or line-like structure of vessels. Among existing vascular extraction methods, Hessian-based multiscale filtering has received much attention [1, 5–12]. These methods share a common idea that the images are convolved with 2D or 3D Gaussian filters at multiple scales, and the eigenvalues of the Hessian matrix at each pixel or voxel are analyzed in terms of a response function to determine the shape of the local structures in the images [13]. The response of the Hessian-based multiscale filtering will be strongest when the scale of the filter matches the size of the local structures, which means scale selection is keeping with the strongest response among multiple scales. Thus, local structures can be extracted using the local strongest response. However, the Hessian-based multiscale filtering is only based on geometry structures of the vessels, which will lead to the discontinuous, even fake vascular structures [1]. Furthermore, Gaussian filter convolution with the image tends to blur vessel boundaries and thus makes the scale selection inaccurate. To address these problems, we use morphological top-hat transformation and Hessian-based multiscale filtering to enhance vascular structures in medical images and an improved level set method involving an external constrained term related to the standard deviation of Gaussian function to extract vascular structures from the enhanced images since the blur level of vessel boundaries is associated with [14].
The paper is structured in five sections; Section 2 describes morphological top-hat transformation and Hessian-based multiscale filtering for vessel enhancement. Vessel segmentation with the improved level set method is presented in Section 3. Section 4 provides some experimental results for synthetic and clinical medical images, as well as evaluation of robustness and segmentation accuracy of the proposed method. The conclusions and future directions are given in Section 5.
2. Vessel Enhancement
The presence of numerous nonvascular structures in clinical medical images such as liver and kidney, will negatively affect the extraction of vascular structures. Considering that morphological top-hat transformation is a powerful technique for image enhancement, especially in extracting bright features from a dark background [15–19], it is adopted in our method to suppress nonvascular structures by using a structuring element larger than the maximum vessel scale in the medical images.
Following the morphological top-hat transformation, the Hessian-based multiscale filtering [5] is used for enhancing the medical image. The filter is based on eigenvalue analysis of the scale space of the Hessian matrix. The eigenvalues and eigenvectors of Hessian matrix are closely related to vascular intensity and direction. For a 3D input image , Hessian matrix is a 3 × 3 matrix composed of second-order partial derivatives of the input image :
In practice, the second-order partial derivatives of input image at a point are defined as a convolution with derivatives of Gaussian filter at scale : where denotes a Gaussian convolution kernel at scale . Let the eigenvalues of be , , and , and their corresponding eigenvectors be , , and , respectively. Table 1 summarizes the relations between and orientations of different structures in the image.
Based on the eigenvalues of , the dissimilarity measures and are defined as
The first ratio accounts for the deviation for a blob-like structure, but it cannot distinguish between a line-like and a plate-like pattern. The second ratio is applied for distinguishing between plate-like and line-like structures. In order to reduce the response of the background pixels, Frangi et al. used the Frobenius norm of the Hessian matrix to define the measure of “second-order structureness” [5] as follows:
Assuming a bright blood image, the vesselness function can be defined as follows: where , , and are thresholds which control the sensitivity of the filter to the measures , , and . The filter is applied at multiple scales, and the maximum response is selected to be a final estimate of vesselness. Consider where and are the minimum and maximum scales at which relevant structures are expected to be found. The choice of the two values must ensure that they will cover the range of vessel widths [5]. The maximum response output is the enhanced image which corresponds to line-like structures.
For 2D images, we use the following vesselness measure defined by Frangi et al. [5] which follows from the same reasoning as used for 3D: where is the blobness measure in 2D images.
3. Vessel Segmentation
3.1. Active Contour Methods
Active contour methods have been popular in a wide range of problems including visual tracking and image segmentation since they were first proposed in 1988 by Kass et al. [20]. The basic idea of active contour methods is to evolve a curve from a given initial state towards an object boundary. In this paper, we used region-based active contour methods under the level set framework to segment vascular trees. Region-based active contour methods assume intensities of image regions to be constants. Compared with edge-based methods, region-based approaches have many advantages such as robustness and insensitivity to image noise [21].
In region-based active contour methods, a curve is iteratively evolved by optimizing an objective function to find the boundary of an object. Let the bounded open subset represent the image domain. Each image is defined as , and is a spatial variable representing a single point within the image domain . In the level set framework, the boundary is embedded in a higher-dimensional function. For example, a simple curve on a 2D plane can be embedded in a 3D surface . By convention, is represented as the zero-level-set of such that the curve is located where crosses a plane at the 0th level. Thus, on the interior of , , and on the exterior of , . During the segmentation process, the function is evolved rather than explicitly evolving the boundary itself when using a parametric boundary representation. The level set evolution equation is given by where is speed function. In our implementation, is initially represented as a signed distance function of the boundary and is evolved via the optimization of an objective function representing the goal of segmentation.
It is well known that level set methods are the most widely used way to represent a contour because of their simple implementation. In addition, it allows very complex curve behavior and automatic topology adaptation [22]. But the primary drawback of level set methods is that they are slow to compute. In this work, we borrow the idea of Lankton’s work [22] to implement our approach with the sparse field method (SFM) proposed by Whitaker [23] which allows one to implement level set active contours efficiently. The objective function used in this work for vascular tree extraction will be described in Section 3.2.
3.2. Vascular Tree Segmentation
In this paper, we use Chan-Vese model [24], a region-based active contour model, to overcome the drawback of Hessian-based multiscale filtering. The object of Chan-Vese model is to minimize the objective function defined by where and denote the average intensity of pixels inside and outside (i.e., , ), respectively. The Heaviside function and the Dirac Delta function will be used to partition the level set function. The objective function can be rewritten in terms of as
Since Gaussian convolution tends to blur vessel boundaries and makes scale selection inaccurate, and the blur level of vessel boundaries is associated with standard deviation of Gaussian function, we define a new class of external constrained term related to . is the penalty on the evolution distance from the initial contour which is obtained by using Otsu’s thresholding [25] method on the enhanced image . Here, is defined as where is the maximum standard deviation of Gaussian function used in Hessian-based multiscale filtering and is a preset value for the largest contour evolution distance ( is set to , in our paper). is the evolution distance in a point from the initial contour .
By combining with , the new objective function is represented as and it will be rewritten in terms of as follows: The evolution follows that when evolution distance in point , , , and , where the contour at point can be evolved as the Chan-Vese model [24]. Meanwhile, when , , , and , where the contour at point will stop evolution. In other words, is the largest permissible evolution distance (i.e., ). Obviously, can control the contour evolution to avoid the segmentation leakage of nonvascular structures and blurred boundaries caused by Gaussian convolution.
If one regularizes the Heaviside function and the Dirac Delta function by two suitable smooth functions and , the evolving equation will be obtained as with the natural boundary condition [26]. In addition
In general, , , , and are fixed parameters. As suggested in [24], these parameters are set as , . Then the evolving equation can be simplified as
We define regularizations of and (where ) as The values used in the simulations are with denoting the space step.
In each step, the should be reinitialized to be the signed distance function [24, 26]. This procedure prevents the level set function from becoming too “flat” due to the use of the regularized Dirac Delta function [26]. The reinitialization process is expressed as The solution of this equation will have the same zero-level-set as , and will converge to 1 since it should be a distance function.
The process of the vascular extraction terminates when the evolution does not change within bounds 0.4 mm2 on successive iterations or the maximum number of iterations is reached. The improved active contour method converges to the boundary of vascular structures exactly in a few iterations.
3.3. Implementation of the Algorithm
(1)Vessel enhancement with morphological top-hat transformation and Hessian-based multiscale filtering; (2)Get the initial contour using Otsu’s thresholding method on the enhanced image ; (3)Initialize based on ; (4)Compute and ; (5)Solve (17); (6)Reinitialize as by using (19); (7)Check whether the solution is stationary or the stopping criteria is met. If not, go back to Step 4; Otherwise stop evolution.
4. Experiments and Discussions
4.1. Experiments on Synthetic Images
To evaluate the performance of the proposed method, we test it on synthetic images containing different vessel-like structures with different diameters and different directions. For quantification, we use the segmentation rate to measure the effectiveness of our method. The segmentation rate is used to estimate the completeness of a segmented vessel, and it is defined as the ratio of the number of segmented pixels to the number of gold standard pixels whose coordinates are known in synthetic images.
4.1.1. Evaluation of Segmentation Accuracy and Robustness
Figure 1 shows three types of synthetic images of size 512 × 512. In Figure 1(a), the diameters of the simulated vascular structures range from 1 pixel to 15 pixels. In Figure 1(b), directions of the vascular structures are simulated by counter-clockwise rotation with an interval of 30 degrees starting from the vertical direction. In Figure 1(c), the intensities of the vessels from left to right are set between 32 to 256 with an increment of 32, while the intensity of background is 0. Meanwhile, five vascular tree models with different complexities are presented in Figure 2.
(a)
(b)
(c)
(a)
(b)
(c)
(d)
(e)
Figure 3 shows the segmentation rate with different vessel diameters, vessel directions, vessel intensities, and vessel complexities. It can be seen that the proposed method can provide accurate segmentation results, and the segmentation rates for all the synthetic images are over 95%. Moreover, the segmentation accuracy is insensitive to vessel diameter or vessel direction and all the cases could be segmented completely. Unfortunately, it has a limit on vessel intensity. If a vessel structure is not obvious compared with the background in the image, the proposed method cannot obtain the expected outcome as shown in Figure 3(c). When the intensity of the vessel structures is lower than 64 with the background intensity equal to 0, the segmentation rate is close to 0. But this kind of vessel structure is rare in clinical practice. Simultaneously, with the increasing complexity of vessel structure, the segmentation rate presents a slight downward tendency as shown in Figure 3(d). Therefore, the proposed method is effective for the extraction in images with vessel-like structures.
(a)
(b)
(c)
(d)
To investigate the sensitivity of the proposed method to noise, we used the synthetic image of size 256 × 256 with a vessel-like structure of varying width and orientation in Figure 4(a), and added zero mean Gaussian noise of standard deviations ranging from 5 to 30 to this image. Figure 4(b) shows the segmentation rate for Gaussian noise of the various standard deviations. It is easy to see from Figure 4(b) that the segmentation rate decreases slightly but remains above 97% with increasing noise levels in the image. Figure 5 shows the segmentation results under different noise levels. Obviously, the proposed method is robust to noise in that it can extract the tree structure effectively at the various noise levels.
(a)
(b)
(a)
(b)
(c)
(d)
(e)
4.1.2. Comparison of Vessel Segmentation Methods
In this section, we compared the segmentation results of the proposed method with those of other two vessel segmentation techniques, Hessian-based multiscale filtering [5] and Hessian-based multiscale filtering combined with Chan-Vese model [24], on the synthetic image. The synthetic image of size 512 × 512 used in this part is given in Figure 6(a). Figure 6(b) shows that the Hessian-based multiscale filtering can locate vessel structures accurately but with inaccurate scales, which means that it is suitable to generate the initial contour. Hessian-based multiscale filtering combined with Chan-Vese model can converge to the boundary of the vascular structure exactly, but it will leak into neighboring nonvascular structures where the contrast is low as shown in Figure 6(c). The result of the proposed method presented in Figure 6(d) is of high accuracy and completely unaffected by nonvascular structures because of the introduction of a new class of external constrained term to penalize the evolution distance of the contour. The segmentation rate of Hessian-based multiscale filtering, Hessian-based multiscale filtering combined with Chan-Vese model, and the proposed method are 88.71%, 92.32%, and 98.14%, respectively.
(a)
(b)
(c)
(d)
4.2. Experiments on MRA Images
In this section, we applied the proposed method on 2D slices extracted from a 3D abdomen MRA image. The image size is 512 × 512 × 60 voxels with pixel spacing 0.51 mm × 0.51 mm × 1 mm which is acquired from syngo MR B15 by routine clinical scan. The 2D slices were generated by slicing through the 3D image in the direction of the coronal plane with 3D Quantify (a multiplanar visualization software) [27]. For the clinical images, there is no “ground truth” to prove presence or absence of the vessel structures or their sizes or positions. Thus, we evaluated segmentation results of the proposed method and the other two methods by visual inspection.
Figure 7(a) shows one of the 2D slices used in our experiments. We compared our proposed method with the Hessian-based multiscale filtering and Hessian-based multiscale filtering combined with Chan-Vese model. The results are presented in Figures 7(b), 7(c), and 7(d). Figure 8 shows a partially enlarged view of the marked green box in Figure 7. Obviously, Hessian-based multiscale filtering cannot estimate the scales of the vessel structures exactly inferred from the segmentation of the main renal artery on the left of Figure 8(b). Similar to Figure 6(c), the result of Hessian-based multiscale filtering combined with Chan-Vese model leaks into adjacent nonvascular structures due to the low contrast as shown on the right side of Figure 8(c). Figure 8(d) demonstrates that the proposed method successfully and accurately extracts most of the vascular structures.
(a)
(b)
(c)
(d)
(a)
(b)
(c)
(d)
5. Conclusions
This paper presented an automatic technique for extracting vascular tree in medical images. Distinctively, the proposed method introduces an external constrained term based on used in Hessian matrix with Gaussian function convolution into the level set to avoid the segmentation leakage of nonvascular structures. In the evaluation based on synthetic images, the segmentation rate of the proposed method is over 97% and it is robust to noise. The results for clinical datasets demonstrate that the proposed method is suitable and effective for the extraction of vascular tree in medical images. The main drawback of our method is that it cannot obtain expected results when the image contrast is very low. Future work will concentrate mainly on optimizing the performance on low contrast images and accuracy evaluation of our method for clinical datasets.
Acknowledgments
This work was partly supported by the National Natural Science Foundation of China (NSFC) (Grant no. 30911120497), the National 973 project (Grant no. 2011CB933103), the Project of the National 12th-Five Year Research Program of China (Grant no. 2012BAI13B02), and Graduates’ Innovation Fund of Huazhong University of Science and Technology (Grant no. HF-11-08-2013).