Distributed under Creative Commons Attribution License Contents 1 Intensity enhancement 2

We present a semi-automatic algorithm for Carotid lumen segmentation on CTA images. Our method is based on a variant of the minimal path method that models the vessel as a centerline and boundary. This is done by adding one dimension for the local radius around the centerline. The crucial step of our method is the definition of the local metrics to minimize. We have chosen to exploit the tubular structure of the vessels one wants to extract to built an anisotropic metric giving higher speed on the center of the vessels and also when the minimal path tangent is coherent with the vessel’s direction. Due to carotid stenosis or occlusions on the provided data, segmentation is refined using a region-based level sets.

tubular anisotropy segmentation method [1] on the enhanced image, in order to obtain a first segmentation, as described in section 2. The tubular anisotropy method represents vessels with circular cross sections. This may yield various inaccuracies in the provided datasets, which contain vessels exhibiting non circular lumen and strong lumen occlusions. Therefore, we propose to refine the segmentation using a region based level set model, as described in section 3. The method presented in this paper was developed for the MICCAI'09 workshop 3D Segmentation in the Clinic: A Grand Challenge -Carotid Bifurcation Lumen Segmentation and Stenosis Grading [3]. Results are discussed in section 4.

Intensity enhancement
Due to the injection of contrast agent before imaging, the vessel lumen is brighter than neighboring structures like muscles and fat. It remains however darker than bone. As a preprocessing step before computing the metric, we enhance intensities around a reference intensity I ref , making the reasonable assumption than intensity within the lumen follows a gaussian distribution. Working with 256-grayscale images, the transformation is implemented through a look-up table computed with For the reference intensity, we take advantage of the provided initial points. The intensity at the point in the Common Carotid Artery turns out to be a relevant choice. The variance was fixed to σ = 20.  2 Tubular anisotropy Tubular Anisotropy Segmentation [1] unable us to extract centerline and boundary of a vessel from an image using source points and destinations. A vessel is modeled by a higher dimensional path [5]. This is done by incorporating an extra non-spatial dimension into the search space, that is the radius of a moving sphere along vessel's centerline. Each point of the 4D path (after adding the extra dimension for the 3D image) consists of three spatial coordinates plus a fourth coordinate which describes the vessel thickness at that corresponding 3D point. Thus, each 4D point represents a sphere in 3D space, and the vessel is obtained by taking the envelope of these spheres as we move along the 4D curve. Then a minimal path approach is employed to obtain the segmentation. A minimal path is the shortest path between two points according to a given metric. A crucial step of this method is to build an adequate metric that drives the propagation. We choose to design an anisotropic metric from the image that takes into account vessel orientation.
In section 2.1, we will give some details on the chosen minimal path model, then in section 2.2, we will show how the metric is constructed and we will give details on the chosen parameters.

Short background on minimal path
A minimal path is a pathway, parametrized along its length (i.e γ = 1), minimizing the energy functional, (1) In this paper, we are interested only on the particular case of Riemannian manifold, that is the potential in under the form P (γ(.), γ (.)) = γ (.) T M (γ(.))γ (.) describing an infinitesimal distance along a pathway γ relative to a metric tensor M (symmetric definite positive). Thus, we are considering only the case of an elliptic medium. A curve connecting p 1 to p 2 that globally minimizes the above energy (1) is a minimal path between p 1 and p 2 , noted C p 1 ,p 2 and also called a geodesic.
The solution of this minimization problem is obtained through the computation of the minimal action map U : Ω → R + associated to p 1 on the domain Ω which is a 4D domain in our case. The minimal action is the minimal energy integrated along a path between p 1 and any point x of the domain Ω : where A p 1 ,x is the set of paths linking x to p 1 . The values of U may be regarded as the arrival times of a front propagating from the source p 1 with oriented velocity related to the metric tensor M −1 . The map U has only one local minimum, the point p 1 , and its flow lines satisfy the Euler-Lagrange equation of functional (1).
Thus, the minimal path C p 1 ,p 2 can be retrieved with a simple gradient descent on U from p 2 to p 1 , solving the following ordinary differential equation with standard numerical methods like Heun's or Runge-Kutta's : Details on the computation of the minimal action map U can be found in [1].
For vessel segmentation, it is not natural to consider orientations on the 4 th dimension, i.e the radii dimension. Thus one propose to decompose by block the metric M as follows: whereM (x, r) is a 3 × 3 symmetric definite positive matrix giving the spatial anisotropy and P radius (x, r) is the radius potential (also strictly positive).

How to construct the metric
The used metric is based on the Optimally Oriented Flux [4] (OOF). At the position x on an image I, the amount of the image gradient projected along the axis v flowing out from a 3D sphere S r is measured as in where G is a Gaussian function with a scale factor of 1 voxel, r is the sphere radius, h is the position vector along ∂S r and da is the infinitesimal area on ∂S r . To detect vessels having higher intensity than the background region, one would be interested in finding the vessel direction which minimizes f (x, v; r), i.e. we are looking for: Using the divergence theorem, it can be shown that f (x, v; r) is a quadratic form on v and its associated matrix, called oriented flux matrix, can be calculated using a simple convolution, where (∂ i, j G) is the Hessian matrix of function G and 1 S r is the indicator function inside the sphere (or circle) S r . Q is called oriented flux matrix. By differentiating the above equation with respect to v, minimization of function f is in turn acquired as solving a generalized eigenvalue decomposition problem. Solving the aforementioned generalized eigen decomposition problem gives 3 eigenvalues, To handle the vessels having various radii, a multi-scale approach should be used along with the OOF method. Then, we normalize the OOF's eigenvalues by the sphere surface area when the OOF method is incorporated in a multi-scale approach for 3D image volumes. The eigenvalues are normalized by the sphere area 4πr 2 .
In practice we compute the oriented flux matrix Q on the downsampled image at different scales, using the enhanced image I output . Here, the minimal and maximal values of the radius, r min and r max , are fixed. We took r min equal to the smallest spacing of the downsampled image and r max = 3.5mm. And the spacing for the radius is equal to half the smallest spatial spacing.
In order to make the minimal path coincides with the vessel centerline, the metric is crucial. The spatial metricM has to be well oriented along the vessel centerline. And the radius potential P radius has to be small for the adequate scale for any point of the image. √ P radius corresponds to the inverse speed for the radius dimension. SinceM is symmetric definite positive, one can decompose it as follows: where 0 < m 1 ≤ m 2 ≤ m 3 are the eigenvalues and u i are the associated eigenvectors. The velocity of the propagating front along direction u i is equal to 1/ √ m i . Therefore, we used the OOF descriptor and the enhanced image I output to construct the metric as follows: The chosen metric for the challenge is slightly different than the original one in [1]. We added a contrast enhancer, (I output + ε) −1 , based on the image itself, where ε = 0.1. The constant α is controlled by an intuitive parameter, which is the maximal spatial anisotropy ratio, noted µ. It corresponds to the maximal spatial speed ratio one wants to impose: By choosing the maximal spatial anisotropy ratio µ, the constant α is fixed. And by doing so, the anisotropy descriptor M becomes affine contrast invariant because the OOF is linear on the image. The parameter β controls the radius speed. For our experiments, parameters α and β where fixed such that µ = 6 and such that the propagation along the radius dimension, 1/ √ P radius is at least four times the best spatial propagation speed. This yields to The tubular anisotropy method yields continuous representations of paths, since point coordinates and corresponding radii are stored as real numbers. As a post-processing step, in order to make output data suitable for evaluation, generated paths are embedded into the discrete integer space of images. At each point located on the path, an ellipsoid is voxelized at the nearest integer point. Radii of the ellipsoid are determined thanks to the real radii of the continuous path and voxel spacings provided in the header data.

Refining segmentation with region-based level sets
The minimal path is computed on a downsampled image and is voxelized into the full image afterwards. In addition, it can only represent vessels with circular cross sections. This may yield various inaccuracies in the provided datasets, which contain vessels exhibiting high curvature and non circular lumen cross sections. Hence, we perform a final refinement step using level-set based segmentation, where the previously discretized path is used as the initial surface. We consider the level set function φ : D ⊂ R 3 → R. The surface is the zero level set of ψ. We define the region enclosed by the surface by R in = {x|ψ(x) ≤ 0}. Function ψ deforms according to the evolution equation: where speed function F is a weighted sum of smoothness and data terms: where coefficient ω controls the significance of the smoothness term and was fixed to 0.5. In the level set framework, regularization is usually performed with a curvature-dependent term. With this technique, the effect is limited to the direct neighborhood of pixels. In order to achieve a more diffuse regularization, we replace the usual curvature term with a gaussian convolution, as in [6][7]: where W is a circular window of a given radius around x. Vessels are brighter than neighboring structures like muscles and fat. With respect to intensity, they may be segmented using a straightforward region criterion. We consider the data term of the Chan-Vese model [2], which we use on the enhanced image.
where the Heaviside step function H selects image points inside or outside the surface. The data speed term is determined from the variational derivative of energy E data with respect to ψ. Thanks to gradient descent, variables k in and k out are assigned to average intensities inside and outside the current region, respectively. One may note that for implementation purpose, we use regularized versions of the Heaviside H and Dirac δ functions, as in [2]. The level set function is updated according to the narrow band technique with a single pixel-wide band. Unlike the minimal path, which is computed on a downsampled image, level set-based refinement is performed in the initial resolution to achieve accurate segmentation.

Results
Comparing coarse segmentations depicted in figure 2 on the one hand and refined segmentations, it turns out that refining the volume obtained with the minimal path approach significantly improves fidelity to actual vessel boundaries. As a drawback of using unconstrained level sets, we should point out the fact that a tradeoff has to be made between boundary fitting and leakage preventing. On a few datasets, we found difficult to segment the carotids accurately on all slices and prevent the surface from propagating into neighboring or branching vessels simultaneously. This phenomenon is depicted in fig. 3 bottom, where the refined surface exhibits bumps corresponding to partial propagations into small branching vessels. Since level sets handle topological changes in a natural way, they allow the carotid volume to split during evolution. As an example, fig. 3   The 31 competition data sets were segmented using the method described above. Intermediate and final segmentation results are shown on figures 2 and 3. The obtained segmentations were evaluated using four performance measures: dice similarity index, mean surface distance, root mean square surface distance and the maximum surface distance, see [3] for more details. A summary of lumen segmentation scores is presented in Table 1. Averages lumen scores are presented in Table 2. The average processing time for a single data set was approximately 2 minutes.

Conclusion and discussion
A semi-automatic method for Carotid lumen segmentation was presented. We tested our method on 31 CTA data sets within the scope of the CLS'2009 MICCAI Challenge. Our dice mean score is 83.6%, making our method relatively robust. The tubular anisotropy method provides us a good first estimate of Carotid lumen. However, as discussed previously, it detects vessels with circular cross sections yielding to inaccuracies on strong lumen occlusions. Therefore, we used a region-based level sets approach in order to refine our segmentation. For future work, we would like to improve the metric for the tubular anisotropy method.