Cerebrovascular segmentation for MRA data using level sets

International Congress Series 1256 (2003) 246 – 252

Cerebrovascular segmentation for MRA data using level sets Hossam Hassan *, Aly A. Farag Computer Vision and Image Processing Laboratory, University of Louisville, Lutz Hall, Louisville, KY 40292, USA Received 14 March 2003; received in revised form 14 March 2003; accepted 18 March 2003

Abstract In this paper, we use a level set based segmentation algorithm to extract the vascular tree from Phase Contrast Magnetic Resonance Angiography, ‘‘PCMRA’’. Classification model finds an optimal partition of homogeneous classes with regular interfaces. Regions and their interfaces are represented by level set functions. The algorithm initializes level sets in each image slice using automatic seed initialization and then iteratively, each level set approaches the steady state and contains the vessel or non-vessel area. The results are validated using a phantom that simulates the ‘‘PCMRA’’. The approach is fast and accurate. Results on various cases demonstrate the accuracy of the approach. D 2003 Published by Elsevier Science B.V. Keywords: Level set segmentation; Vascular tree; Seed initialization; PCMRA; PDE

1. Introduction The human cerebrovascular system is a complex 3-D anatomical structure. Serious types of vascular diseases such as carotid stenosis, aneurysm and vascular malformation may lead to brain stroke, which are the third leading cause of death and the number one cause of disability. An accurate model of the vascular system from MRA data volume is needed to detect these diseases at early stages and hence may prevent invasive treatments. A variety of methods have been developed for segmenting vessels within MRA. One class of methods is based on a statistical model, which classifies voxels within the image volume into either vascular or nonvascular class for time-of-flight MRA [1]. Another class of segmentation is * Corresponding author. Tel.: +1-502-852-2789; fax: +1-502-852-1580. E-mail addresses: [email protected] (H. Hassan), [email protected] (A.A. Farag). URL: http://www.cvip.uofl.edu. 0531-5131/03 D 2003 Published by Elsevier Science B.V. doi:10.1016/S0531-5131(03)00445-X

H. Hassan, A.A. Farag / International Congress Series 1256 (2003) 246–252

247

based on intensity threshold where points are classified as either greater or less than a given intensity. This is the basis of the iso-intensity surface reconstruction method [2– 4]. This method suffers from errors due to image inhomogenities. In addition, the choice of the threshold level is subjective. An alternative to segmentation is axis detection known as skeletonization process, where the central line of the tree vessels is extracted based on the tubular shape of vessels [5]. Other approaches for MRA vessel segmentation are the manually defined seed locations for segmentation [6]. In this paper, we use level set method for image segmentation to improve the accuracy of the vascular segmentation. This work is a supervised classification, which means that the number of classes, and the class distribution are assumed to be known. Usually, the class distribution is assumed to be Gaussian. The classification model taken from Ref. [7] and assumes that the classes are phases separated by interfaces boundaries. A set of functionals are developed with properties of regularity. The level set function representation depends on these functionals. Each class occupies certain areas (regions) in the image. The classes have no common areas, i.e., the intersection between classes is not allowed. The sum of lengths of the interfaces between the areas is taken in consideration. Herein, the functionals are dependent mainly on these properties and they expect to have a minimum, which is the segmented image. The change of each level set is guided by two forces, the minimal length of interfaces which is the internal force, and the homogeneous class distribution which is the external one. A Partial differential equation (PDE) guides the motion of each level set. In this paper, we use level sets for segmentation of phase contrast MR (PCMRA). A related study to segment MRI data has been represented in Ref. [7]. The vessels represent a critical ratio in the PCMRA image with low contrast which represents a challenge in this case. After segmenting the volume, a connectivity filter is used to exploit the fact that the vascular system is a tree-like structure and makes use of the 3-D computer graphics region-filling algorithm to extract the vascular tree. The used algorithm with PCMRA data volumes is evaluated using a phantom showing a good accuracy. We finally show various experimental results. The results are shown using a visualization toolkit ‘‘VTK’’, which can be viewed on a stereo graphics workstation or in a virtual reality environment. Level set [8] evolution combining global smoothness with the flexibility of topology changes offer significant advantages over conventional statistical classification. Sabry et al. [9] used a statistical approach with Expectation Maximization (EM) algorithm for mixture of densities to extract the vascular tree from PCMRA. This approach depends on the grey level. In this paper, we use only the grey level and the interfaces between the class regions which can be extended to include geometric properties. This leads to detecting more details in the vascular tree.

2. Review of related work The level set function is a function defined over the image space. It is a surface that is positive inside the region, negative outside and zero on the interfaces between regions. Ui ðx; y; tÞ > 0 If xaXi ; Ui ðx; y; tÞ ¼ 0 If xaCi Ui ðx; y; tÞ < 0 otherwise: Xi is the region and Ci is the interface between two regions [7].

ð1Þ

248

H. Hassan, A.A. Farag / International Congress Series 1256 (2003) 246–252

The level set method is an analytic framework for working with evolving geometries. The level set method embeds geometries into a scalar field function. Its variation with time can be described by a PDE: BUi ¼ FArUi A ð2Þ Bt F is an extremely imposed boundary propagation function or velocity function. Also, F is application dependent. Eq. (2) will take a more detailed form. By defining F such that it incorporates both an image dependence as well as interface boundaries dependence, we can cause boundaries to be attached to image features subject to some constraints. The steady-state solution for each level set represents the segmented area in the positive part. The classification model satisfies the following three conditions: .

.

.

The partition condition: Z K k X F1 ¼ ðHa ð/i Þ  1Þ2 dx; kaRþ and K is the number of classes 2 X i¼1

ð3Þ

The data term condition (Gaussian property of the class): Z K X ðuo  ui Þ2 ei Ha ð/i Þ dx; F2 ¼ B2i X i¼1 ei aRbi ; ui and B2i are the mean and variance The length shortening of interfaces between regions: Z K X F3 ¼ ci da ð/i ÞArUi Adx; ci aRbi i¼1

ð4Þ

ð5Þ

X

H function represents the Heaviside distribution function and d represents the Dirac distribution function in Fig. 1. The first represents the entire region and the second represents the boundary of the region. As a ! 0+, the sum F = F1 + F2 + F3 is minimized and taking the derivative of F with respect to Ui and set to 0, we get the solution represented by Eq. (6). This solution represents the level set function variation with time. When the function approaches the steady state, it is then not changed. It has positive parts,

Fig. 1. Plot of H and d functions for a = 0.5.

H. Hassan, A.A. Farag / International Congress Series 1256 (2003) 246–252

249

Fig. 2. Average histogram for the first patient data set.

negative and zero parts. We are interested only in the positive parts. Each pixel in the positive parts belongs to the associated class of its function. " !#   K X ðuo  ui Þ2 jUi tþ1 tþ1 Ui ¼ Ui  Dtda ðUi Þ ei  ci div Ha ðUi Þ  1 ð6Þ þk AjUi A B2i i¼1 Dt is the step in time. By this representation, jUi the level set function formulation allows breaking and merging fronts where div AjU represents the local curvature of the level iA set function. Signed distance function is commonly used for level set schemes where the representation defined by Eq. (6) does not maintain this condition. jjUij = 1 is the condition for the level set function to be signed distance. This condition is maintained through time for many practical reasons. From Ref. [10], we find that during image segmentation, it has been found that the level sets function can be changed into a nondistance function with the initial level set function defined as a distance function. This

Fig. 3. Automatic seed initialization of slice #60.

250

H. Hassan, A.A. Farag / International Congress Series 1256 (2003) 246–252

Fig. 4. Slices #1, 30, 60 and 117.

causes some applications to fail, such as Coupled Surfaces Propagation. Yu et al. [10] proved that the signed distance function could be preserved to the level set function and they developed a representation for the level set function to maintain this property. But the analysis introduced in Ref. [10] is not suitable to use since our application is different. Another solution found in Ref. [11] is to recompute the level set function periodically based on Eq. (7). In our implementation, for each n iterations, we apply this equation to regularize the level set function to remain signed distance function. BUi Ui ¼ ð1  AjUi AÞ Bt AUi A

ð7Þ

3. Volume segmentation algorithm 3.1. Algorithm The algorithm in Ref. [7] is modified for segmenting the volume, and also, we added a smoothing layer to enhance the segmentation areas:

Step Step Step Step Step Step Step

0: 1: 2: 3: 4: 5: 6:

Initialize each level set function. Initialize Counter to 0. Counter = Counter + 1. Update each level set function using Eq. (6). Update each level set function using Eq. (7) for each n iterations. Smooth each level set function and remove noise. If steady state is not reached, then go to Step 2, else go to next slice.

The steady state means that the difference between two iterations is not significant or a certain number of iterations is reached. Step 0 uses automatic seed initialization to speed up

Fig. 5. Segmentation results for slices #1, 30, 60 and 117.

H. Hassan, A.A. Farag / International Congress Series 1256 (2003) 246–252

251

Fig. 6. Visualization of data set #1.

the process and it is also less sensitive to noise. Automatic seed initialization [7] is used do divide the image into nonoverlapped windows of predefined size. Then, the average grey level is calculated and compared to the mean of each class to specify the nearest class it belongs to. A signed distance function is initialized to each window. After segmentation, the connectivity filter is applied. 3.2. Segmentation quality measurement A 2-D phantom is designed to simulate the PCMRA. This phantom image contains many circles with decreasing diameters like the cerebrovascular tree shape, which is a coneshaped. Then, using the level set segmentation algorithm with this image we obtain a resultant image. By counting the number of white pixels (segmented vessels) in the resultant image and dividing it by the number of pixels in the phantom image (in the circular areas) we get the evaluation factor. A value of 0.94 was computed, indicating a high-quality segmentation of the phantom image.

4. Results MRA data set consists of 117 slices. Each slice has a size of 256 256. Automatic seed initialization is used in each slice and each slice is divided into windows of size 5 5. Our

Fig. 7. Visualization of the data set #2.

252

H. Hassan, A.A. Farag / International Congress Series 1256 (2003) 246–252

experiments are made for two data sets. An average mean is estimated for each class from the average histogram of the volume. Fig. 2 shows the average histogram for the first data set. Fig. 3 shows the automatic seed initialization for slice #60 where the white color represents the vessels and the dark gray represents the non-vessel tissues. Fig. 4 shows slices #1, 30, 60 and 117 (for the patient #1) and the segmentation results are shown in Fig. 5 before applying the connectivity filter. Applying the connectivity filter and using VTK the two data sets are visualized and shown in Figs. 6 and 7. It is clear that patient #2 suffers from a malformation, which is the accurate diagnosis as produced by a radiologist. The connectivity filter removed the non-vessels tissues that appeared in the segmentation.

5. Conclusion and future work We have presented a fast and accurate supervised level set based segmentation algorithm to extract the cerebrovascular system from phase contrast MR angiography. The results are shown using VTK and a connectivity filter was used to remove tissues that form small islands in the segmented volume. A phantom is used to evaluate the algorithm and shows a good accuracy. The approach provides comparable results to that reported in Ref. [9]. The use of the phantom will help provide a quantitative evaluation of the algorithm. The evaluation criteria will include registration. The 2-D phantom can be modified to be a 3-D phantom simulating the whole volume leading to more accuracy.

References [1] D.L. Wilson, J.A. Noble, An adaptive segmentation algorithm for time-of-flight MRA data, IEEE Trans. Med. Imag. 18 (10) (1999) 938 – 945. [2] H.E. Cline, W.E. Lorensen, R. Kikinis, R. Jolesz, Three-dimensional segmentation of MR images of the head using probability and connectivity, Neurosurgery 14 (1990) 1037 – 1045. [3] S. Nakajima, H. Atsumi, A.H. Bhalerao, et al., Computer-assisted surgical planning for cerebrovascular neurosurgery, Neurosurgery 41 (1997) 403 – 409. [4] H.E. Cline, W.E. Lorensen, S.P. Souza, F.A. Jolesz, R. Kikinis, G. Gerig, T.E. Kennedy, 3D surface rendered MR images of the brain and its vasculature, JCAT 15 (1991) 344 – 351. [5] P.J. Yim, P.L. Choyke, R.M. Summers, Gray-scale skeletonization of small vessels in magnetic resonance angiography, IEEE Trans. Med. Imag. 19 (6) (2000) 568 – 576. [6] E. Bullitt, et al., Symbolic description of intracerebral vessels segmented from magnetic resonance angiograms and evaluation by comparison with X-ray angiograms, Med. Imag. Anal. 5 (2001) 157 – 169. [7] C. Samson, L.B. Feraud, G. Aubert, J. Zerubia, Multiphase Evolution and Variational Image Classification, INRIA Research Report, April 1999. [8] J.A. Sethian, Level Set Methods and Fast Marching Methods, Cambridge Univ. Press, USA, 1996, 1999. [9] M. Sabry, C.B. Sites, A.A. Farag, S. Hushek, T. Moriarty, Statistical cerebrovascular segmentation for phase-contrast MRA data, Proc. of the 1st International Conf. on Biomedical Engineering, Cairo, Egypt, December, 2002. [10] H. Yu, D. Wang, Z. Tang, Level set methods and image segmentation, IEEE, International Workshop on Medical Imaging and Augmented Reality (MIAR’01) 0-7695-1113-9/01, June, 2001. [11] M. Sussman, P. Smereka, S. Osher, A level set approach for computing solutions to incompressible twophase flow, J. Comput. Phys. 114 (1994) 146 – 159.