Vascular Tree Segmentation in Medical Images Using Hessian-Based Multiscale Filtering and Level Set Method

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.


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 tubularlike or line-like structure of vessels. Among existing vascular extraction methods, Hessian-based multiscale filtering has received much attention [1,[5][6][7][8][9][10][11][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.

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][16][17][18][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 : ] . (1) 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 ∇ 2 be 1 , 2 , and 3 (| 1 | ≤ | 2 | ≤ | 3 |), and their corresponding eigenvectors be 1 , 2 , and 3 , respectively. Table 1 summarizes the relations between and orientations of different structures in the image. Based on the eigenvalues of ∇ 2 , the dissimilarity measures and are defined as The first ratio accounts for the deviation for a bloblike 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 min and max 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. Medicine   3 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:

Computational and Mathematical Methods in
where = 1 / 2 is the blobness measure in 2D images.

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, regionbased 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 Ω ⊂ 2 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 higherdimensional 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 , < 0, and on the exterior of , > 0. 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.

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 cv ( 1 , 2 , ) defined by where 1 and 2 denote the average intensity of pixels inside and outside (i.e., 1 = 1 ( ), 2 = 2 ( )), respectively. The Heaviside function and the Dirac Delta function 0 will be used to partition the level set function. The objective function can be rewritten in terms of ( , ) as (1 − ( , )) .
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 0 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 max /2, in our paper). | ( , ) − 0 | is the evolution distance in a point ( , ) from the initial contour 0 .
By combining cv with , the new objective function is represented as and it will be rewritten in terms of ( , ) as follows: 4
If one regularizes the Heaviside function and the Dirac Delta function 0 by two suitable smooth functions and , the evolving equation will be obtained as with the natural boundary condition [26]. In addition In general, ≥ 0, ] ≥ 0, 1 > 0, and 2 > 0 are fixed parameters. As suggested in [24], these parameters are set as = 0, 1 = 2 = 1. Then the evolving equation can be simplified as We define regularizations of and (where ( ) = ( / ) ( )) as The values used in the simulations are = ℎ = 1 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-levelset 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 mm 2 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.

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. 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 counterclockwise 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. 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 vessellike structures.

Evaluation of Segmentation Accuracy and Robustness.
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.

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 Hessianbased multiscale filtering, Hessian-based multiscale filtering  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 Hessianbased 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.

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.