Probabilistic Tracking and Model-based Segmentation of 3D Tubular Structures Stefan Wörz, William J. Godinez, Karl Rohr University of Heidelberg, BIOQUANT, IPMB, and DKFZ Heidelberg, Dept. Bioinformatics and Functional Genomics, Biomedical Computer Vision Group s.woerz@dkfz.de Abstract. We introduce a new approach for tracking-based segmentation of 3D tubular structures. The approach is based on a novel combination of a 3D cylindrical intensity model and particle filter tracking. In comparison to earlier work we utilize a 3D intensity model as the measurement model of the particle filter, thus a more realistic 3D appearance model is used that directly represents the image intensities of 3D tubular structures within semi-global regions-of-interest. We have successfully applied our approach using 3D synthetic images and real 3D MRA image data of the human pelvis. 1 Introduction Accurate segmentation and quantification of 3D tubular structures is crucial in many biomedical applications. Examples of such structures are blood vessels, airways, or the spinal cord. In particular, the accurate quantification of human blood vessels is indispensable for diagnosis, treatment, and surgical planning. In clinical practice, the human vascular system is typically imaged by 3D magnetic resonance angiography (MRA) or computed tomography angiography (CTA). Previous work on the automatic segmentation of tubular structures from 3D images can be divided into two main classes of approaches, one based on two-step centerline segmentation schemes and the other based on incremental tracking schemes. In two-step centerline segmentation schemes, first a local vesselness measure or a binary segmentation is determined. In the second step, the centerline of a tubular structure is constructed based on the results of the first step. Such approaches are, e.g., based on differential measures [1], minimal cost paths [2], or region growing with subsequent skeletonization [3]. A disadvantage of two-step approaches is that often only the centerline is computed and, thus, an additional step is required to estimate the local vessel shape (e.g., radius) [3]. In contrast, incremental tracking schemes segment vessels by incrementally proceeding along a vessel. In each increment, the image information of a vessel is analyzed locally, and a prediction of the centerline position of the next increment is generated. Incremental tracking schemes differ w.r.t. the used measurement model and the prediction method. Concerning the prediction, three different methods are typically used. Often, a simple forward propagation method is utilized where the predicted position of the next increment is directly based on 42 Wörz, Godinez & Rohr the segmentation result of the current increment, usually proceeding a short distance along the local direction of a vessel (e.g., [4]). A more robust prediction method is the Kalman filter, which assumes a linear measurement model (e.g., [5, 6]). A promising alternative method are particle filters, which are more general than the Kalman filter (e.g., inclusion of a nonlinear measurement model) and exploit more effectively the image information (e.g., [7, 8, 9]). Regarding the measurement model, for tracking schemes based on forward propagation or a Kalman filter, typically contour information (e.g., [6]) or intensity information (e.g., [4, 5]) is used. For approaches employing particle filters, different measurement models have been used, e.g., based on circular shortest path search [7], gradient-based shape detection [9], or maximization of the image contrast [8]. However, these approaches utilize only relatively coarse appearance models based on mean intensity levels [7, 8, 9]. Also, these approaches often analyze tubular structures in 2D planes orthogonal to the centerline, i.e., the full 3D information is not exploited. In this contribution, we introduce a new approach for the segmentation of 3D tubular structures, which is based on a combination of 3D parametric intensity models and particle filter tracking. In contrast to previous approaches based on particle filters (e.g., [7, 8, 9]), our approach relies on a 3D cylindrical intensity model as the measurement model in conjunction with a model fitting scheme. This 3D appearance model directly represents the image intensities of 3D tubular structures within semi-global regions-of-interest (ROIs). Thus, we utilize an advanced and more realistic appearance model. Also, since we directly exploit the 3D intensity information, a segmentation step or computation of image gradients is not required. Moreover, we directly quantify the tubular structure (e.g., centerline and local shape) based on the integrated model fitting scheme. 2 Materials and methods In our approach, we use a 3D parametric intensity model which represents the shape and the image intensities of a tubular structure, and which serves as the appearance model of the particle filter. The model consists of an ideal sharp 3D cylinder convolved with a 3D Gaussian and incorporates, for example, the blurring effect of the image formation process. This cylindrical model comprises parameters for the width of the tubular structure (radius R) and the image blur σ, and is well-suited to represent tubular structures of different widths. Since the exact solution of a Gaussian smoothed cylinder cannot be expressed in closed form, we use an accurate approximation gCylinder , which is defined as the weighted superposition of two approximations for thin and thick cylinders [5]. Moreover, we incorporate the intensity levels a0 (surrounding tissue) and a1 (vessel) as well as a 3D rigid transform R with rotation parameters α = (α, β, γ) and translation parameters x0 = (x0 , y0 , z0 ). This results in the parametric model gM,Cylinder (x, p) = a0 + (a1 − a0 ) gCylinder (R (x, α, x0 ) , R, σ) where x = (x, y, z) and p = (R, a0 , a1 , σ, α, β, γ, x0 , y0 , z0 ). To segment a certain vessel segment, we use a model fitting approach based on least-squares fitting of the Model-based segmentation of 3D tubes 3D model to the image intensities g(x) within a ROI: X 2 (gM,Cylinder (x, p) − g (x)) → min. x∈ROI 43 (1) For 3D vessel segmentation, we use a particle filter approach where tracking is formulated as a Bayesian sequential estimation problem. We represent the configuration of a vessel by a state vector θ k and assume that the measurement zk reflects the true state of θ k . At the increment k, the aim is to estimate the state θ k given a sequence of measurements z1:k . By modeling the state’s evolution via a dynamical model p(θ k |θ k−1 ) and incorporating measurements derived from the images via a measurement model p(zk |θ k ), a Bayesian filter estimates the posterior distribution p(θ k |z1:k ) via stochastic propagation and Bayes’ theorem: R p(θ k |z1:k ) ∝ p(zk |θ k ) p(θ k |θ k−1 ) p(θ k−1 |z1:k−1 ) dθ k−1 (2) An estimate of θ k can be obtained from the posterior p(θ k |z1:k ), which, in our case, is determined using a particle filter. The main idea of this algorithm is to i s employ a set {θ ik ; wki }N i=1 of Ns weighted random samples θ k (the ‘particles’) to approximate the posterior distribution. Based on the particle filter framework we have developed two different tracking schemes. In the first scheme, the state vector is defined by the position x0 of the 3D tubular model gM,Cylinder and by its velocity vector, thus θ k = (x0,k , x′0,k , y0,k , y ′0,k , z0,k , z ′0,k ). For the dynamical model p(θ k |θ k−1 ) we adopt a constant velocity model. Measurements for p(zk |θ k ) are computed based on the difference between the predicted position and the measured position zk = x̂0,k , which is obtained by fitting the 3D tubular model based on (1). In our second tracking scheme, the parameters p including the position x0 of the 3D tubular model gM,Cylinder constitute the state θ. For the dynamical model p(θ k |θ k−1 ) we assume that the radius R, intensity levels a0 and a1 , the image blur σ, and the rotation parameters α follow independent Gaussian random walk dynamics. For the position x0 , the dynamical model is defined by ¢ ¡ xi0,k = xi0,k−1 + R−1 (0, 0, l)T , αik−1 , 0 + nik , (3) where l denotes the magnitude of the displacement vector and n is a noise vector comprising independent zero-mean Gaussian noise statistics. Our measurement model p(zk |θ k ) quantifies the probability that the measured image intensities zk of the image g conform to the predicted state θ k within a 3D ROI. 3 Results We have applied our new approach to different 3D synthetic images as well as 3D MRA image data. In the first part of the experiments, we have used 3D synthetic images containing different tubular structures (straight and curved structures) with a spectrum of different widths, curvatures, and noise levels. 44 Wörz, Godinez & Rohr From the experiments we found that our two tracking schemes generally yield accurate results for the centerline and shape of the tubular structures. As an example, we have used tubular spirals which comprise several short low contrast parts along the spiral to simulate vessels with poor local contrast (see Fig. 1, left). Both tracking schemes performed well for a range of different noise levels. In the case of a relatively high noise level (normal contrast a = 50, low contrast alow = 10, and added Gaussian noise with standard deviation σn = 30), it turned out that the first tracking scheme terminates in the third winding. In comparison, the second tracking scheme, which directly employs the cylindrical appearance model, is able to segment the full spiral (seven windings) despite the poor signal-to-noise ratio. Fig. 1 shows the segmentation result represented by the centerline (center) and shape (right). It can be seen that the spiral including the highly curved inner part has been fully tracked and generally well segmented. In contrast, a previous approach based on a Kalman filter [5] failed to segment the full spiral, and terminated even earlier than the first tracking scheme based on a particle filter (see the marked termination points in Fig. 1, left). To quantify the results, we have computed the maximal error eR,max of the estimated radius. For the first 1000 voxels along the centerline (about two windings), which have Particle filter (scheme 2) Particle filter (scheme 1) Kalman filter Start Fig. 1. 3D spiral with low contrast parts and indicated end points for the tracking schemes based on the Kalman filter and both particle filter schemes (left) as well as segmentation result of our new approach using the second tracking scheme represented by the centerline (center) and shape (right). Fig. 2. Maximum intensity projection of a 3D MRA of a human pelvis (left) and segmentation result for two arteries using the second particle filter tracking scheme (right). Model-based segmentation of 3D tubes 45 been segmented by all three approaches, we found that the new approach yields a maximal error of eR,max = 1.12 and eR,max = 0.87 voxels for the first and second tracking scheme, respectively, and the Kalman filter yields eR,max = 1.28 voxels. Moreover, we have also applied the new approach using both tracking schemes to 3D MRA images. For example, Fig. 2 shows the maximum intensity projection of the MRA of a human pelvis (left) as well as the 3D segmentation result using the second tracking scheme (right). It can be seen that arteries of varying sizes and high curvatures have been well segmented. 4 Discussion We introduced a new approach for tracking-based segmentation of 3D tubular structures. The approach is based on a combination of a 3D cylindrical intensity model and particle filter tracking. In comparison to earlier work we utilize a 3D intensity model as the measurement model of the particle filter, thus a more realistic appearance model is used that directly represents the image intensities of 3D tubular structures within semi-global regions-of-interest. We have successfully applied our approach using 3D synthetic images and real 3D MRA image data. Acknowledgement. This work has been funded by the Deutsche Forschungsgemeinschaft (DFG) within the project QuantVessel (RO 2471/6-1). References 1. Frangi AF, Niessen WJ, Vincken KL, et al. Multiscale vessel enhancement filtering. Proc MICCAI. 1998; p. 130–137. 2. Wink O, Niessen WJ, Viergever MA. Multiscale vessel tracking. IEEE Trans Med Imaging. 2004;23(1):130–133. 3. Volkau I, Ng TT, Marchenko Y, et al. On geometric modeling of the human intracranial venous system. IEEE Trans Med Imaging. 2008;27(6):745–51. 4. Noordmans HJ, Smeulders AWM. High accuracy tracking of 2D/3D curved line structures by consecutive cross-section matching. Pattern Recogn Lett. 1998;19(1):97–111. 5. Wörz S, Rohr K. Segmentation and quantification of human vessels using a 3-D cylindrical intensity model. IEEE Trans Image Process. 2007;16(8):1994–2004. 6. Guerrero J, Salcudean SE, McEwen JA, et al. Real-time vessel segmentation and tracking for ultrasound imaging applications. IEEE Trans Medical Imaging. 2007;26(8):1079–1090. 7. Florin C, Paragios N, Williams J. Globally optimal active contours, sequential Monte Carlo and on-line learning for vessel segmentation. Proc ECCV. 2006; p. 476–489. 8. Schaap M, Manniesing R, Smal I, et al. Bayesian tracking of tubular structures and its application to carotid arteries in CTA. Proc MICCAI. 2007; p. 562–570. 9. Lesage D, Angelini ED, Bloch I, et al. Medial-based bayesian tracking for vascular segmentation: Application to coronary arteries in 3D CT angiography. Proc IEEE ISBI. 2008; p. 268–271.