pp. X–XX LEVEL SET BASED BRAIN ANEURYSM CAPTURING IN 3D

Brain aneurysm rupture has been reported to be closely related to 
aneurysm size. The current method used to determine aneurysm size 
is to measure the dimension of the aneurysm dome and the width of 
the aneurysm neck. Since aneurysms usually have complicated 
shapes, using just the size of the aneurysm dome and neck may not 
be accurate and may overlook important geometrical information. In 
this paper we present a level set based surface capturing 
algorithm to first capture the aneurysms from the vascular tree. 
Since aneurysms are described by level set functions, volumes, 
curvatures and other geometric quantities of the aneurysm surface 
can easily be computed for medical studies. Experiments and 
comparisons with models used for capturing illusory contours in 2D 
images are performed. Applications to medical images are also 
presented to show the accuracy, consistency and robustness of our 
method in capturing brain aneurysms and volume quantification.

1. Introduction.Subarachnoid hemorrhage, primarily from brain aneurysm rupture, accounts for 5 to 10% of all stroke cases with a high fatality rate [1].Advancements in neuroimaging technology have helped these aneurysms to be more frequently found prior to rupture.A method to determine if aneurysms are at higher risk of rupturing would be extremely valuable.Brain aneurysm rupture has been reported to be related to the size of aneurysms [2].It is known that the risk of rupture greatly increases as the aneurysm becomes larger [3,4].Currently, methods to determine the aneurysm size are to simply manually measure the size of the neck and the dome of aneurysms.However, these methods may overlook important geometric information and are very hard to perform consistently across subjects.Our goal in this paper is to first segment the aneurysm from the entire blood vessel with minimal human interaction, then compute its volume and other geometric quantities.This problem can actually be realized as an illusory surface capturing problem by observing Fig. 1, which is an extension from illusory contours in 2D.The boundaries of the aneurysm that are not part of the blood vessel surface can be completed naturally by illusory surfaces.Illusory contours have been intensively studied in cognitive neuroscience, where people find that the human vision system is capable of combining nonexistent edges and making meaningful visual organization of both the real and imaginary contour segments [5,6,7] (e.g. the Kanizsa square and triangle [5] in Fig. 1).Various researchers have introduced mathematical models and techniques to mimic the human vision system in detecting and capturing perceptual contours in images [8,9,10,11,12,13].These mathematical models can be used to describe the process in medical evaluation when the location of an aneurysm needs to be identified.Given that our problem is to first capture and then calculate volumes and geometries of aneurysms, representing surfaces using level set functions and designing a proper surface evolution PDE is essential.Therefore, we introduce a level set and PDE based illusory surface model, inspired by the illusory contour models in [11], to capture aneurysms, and calculate their volumes and geometries.
The focus of this paper is to introduce a novel method to capture a specific part of a given pre-segmented surface obtained from 3D images.Therefore, we will not place emphasis on the techniques of surface reconstruction from 3D images.However, we note that different segmentations from a 3D image may result in rather different surfaces in terms of geometry.In fact, surface segmentation from 3D images is highly nontrivial, and is a very active research area.Interested readers can consult [14,15,16] and their references for detailed techniques of blood vessel segmentations.As an example, we applied a simple thresholding method (with carefully chosen thresholds) followed by fast sweeping method [19] and Gaussian smoothing to reconstruct the surface represented by a level set function [20], which takes positive values inside the vessel region and negative values outside.We note that one can replace the Gaussian smoothing by some more sophisticated smoothing techniques, e.g.nonlocal means [17,18], if there are some sharp and delicate features in the surfaces need to be well reserved.Fig. 2 illustrates the idea of reconstruction of vessel surfaces from 3D images and an example of brain aneurysm.2. Review of Level Set Based Illusory Contour Capturing.Many PDE based methods have been proposed to identify illusory contours as well as explain the phenomena [8,9,10,11,12,13].One of the most typical ones was introduced by Sarti et al [9,10].In their work, they first chose a fixation point inside the domain bounded by the ideal illusory contour, and constructed a surface on the whole domain on the basis of the point, then evolve the entire surface based on the image gradient.In fact, our user interactive initialization strategy (in Section 3.3) is very similar to their fixation point idea.More details can be found in [10].
In this section, we shall focus on reviewing the level set formulations of the illusory contour problems introduced by Zhu and Chan [11], because our PDE model is motivated by theirs.As a convention, all level set functions in this paper take negative values inside the domain of interest and positive outside.For example in Figure 1, denoting the regions inside the red curves as Ω, then the corresponding level set function φ satisfies and ∂Ω, which is the zero level set of φ, represents the illusory contours (red curves).
The first model considered in [11] is where ψ is the signed distance function obtained from a given image (see e.g.Fig. 1) whose zero level set represents the boundaries of the objects in the image, and d = |ψ| is the corresponding unsigned distance function.The symbol ∇ is the gradient operator, δ(φ) is the Dirac delta functional, and H(φ) is the Heaviside function.The energy term αH(ψ) ensures that the model will capture only the inner contour, instead of the outer one (see e.g. the Kanizsa square and triangle in Fig. 1).We note that in [11], the authors also had an additional term Ω κ 2 dx in the energy (1) to ensure the continuity of the curvature of φ.However, this term will result in a fourth order evolution PDE which significantly slows down the computations.Since what is important for us is the consistency of segmentation of aneurysms and their volumes, it seems unnecessary to have the curvature term in the energy.
From equation (1), the corresponding gradient flow can be written as Since the function δ(φ) is concentrated only on the zero level set of φ, the PDE (2) only describes a motion for the zero level set of φ.Similar to [21], to ensure all level sets of φ have similar motions and to be able to solve the PDE on the entire 3D rectangular domain, we replace δ(φ) in ( 2) by |∇φ| and obtain the following PDE Numerical results in [11] show that solving the PDE (3) gives fairly good results.However, the sharp corners, e.g. the corners in Kanizsa squares and triangles in Fig. 1, are not well captured.Therefore in that paper, the authors also considered the following improved model which enables the final curves to stick to sharp corners more closely where μ is some constant, c a,b is some restriction function defined in (12) and κ + (ψ) is the positive part of curvature.The corresponding evolution PDE to the energy ( 4) is Numerical experiments in [11] show that the model ( 5) does an excellent job in capturing illusory contours, especially at regions with sharp features.

Our Method and Model Comparisons.
In this section, we will first discuss the possibility of extending ( 1) and ( 4) directly from 2D to 3D, and what difficulties and problems one may encounter.Then, motivated by these discussions, we shall introduce our model in Section 3.2.

Direction Extension from 2D to 3D?.
As discussed in [11], the energies (1) and ( 4) are not convex, and hence proper initialization is required.They showed that if the initial contour contains all key features of the object, e.g.corners, one can obtain satisfactory results.However for our particular application here, it is not reasonable to assume that one can always start with some surface that includes the entire aneurysm within or as a subsurface.In addition, it is always desirable to have a method with less restrictive initialization constraints.In this subsection, we will observe that, generally speaking, the energies (1) and ( 4) require more restrictive initializations in 3D than in 2D.We also find that Gaussian curvature is a more suitable quantity to use than mean curvature.
Let us first consider the model ( 1) and call it E 1 .Notice that it can be extended to 3D trivially.Since it gives a fairly good result in 2D, one may expect it to do so in 3D.However, the following special examples says otherwise.Take α = β = 0 in (1), and take ψ represents a unit circle/sphere, and φ represents a circle/sphere with radius r ∈ [0, 1] (see Fig. 3).Then we can write E 1 explicitly as a function of r, The plot of E 1 (r) in Fig. 3 shows that one should initialize φ by a circle with radius r > 1 2 and r > 2 3 respectively in 2D and 3D in order to converge to the right solution.Thus, the initialization constraint in 3D is more restrictive than that in 2D.This is because the shrinking force from minimizing the surface area in 3D is stronger than that from minimizing the curve length in 2D.On the other hand, model (1) does not perform well in terms of capturing sharp corners in both 2D and 3D and the problem is even more severe in 3D than in 2D for similar reasons.Let us now consider (4), and denote it by E 2 .Take the same example as in Fig. 3, and let α = β = 0 and c a,b ≡ 1.If we choose κ to be mean and Gaussian curvature, we shall have the following two formulas for E 2 (r) in 2 and 3 dimensions, The plots of E 2 (r) in Fig. 4 show that if we choose mean curvature for κ, then in 2D, we have global convergence, while in 3D, we still need to initialize properly in order to get the correct solution.However if we choose Gaussian curvature, we have global convergence in both 2D and 3D.This is one of the motivations for using Gaussian curvature in 3D.Another, yet more important, motivation for choosing Gaussian curvature instead of mean curvature for our particular application is illustrated in the following Fig. 5.As is shown below, Gaussian curvature can discriminate between aneurysm and blood vessels naturally, while mean curvature cannot.This is essentially because cylindrical objects have small Gaussian curvatures in general, while the Gaussian curvatures for bulb like objects are relatively large.In contrast, mean curvatures for both cylindric and bulb like objects generally have comparable magnitudes.However, directly applying (4) in 3D using Gaussian curvature for κ does not give satisfactory results in general.We now consider another simple example in 2D to illustrate this.In Section 3.4 we shall give some comparisons in 3D.Let us now consider an example as given in Fig. 6, where the target object is no longer convex.We note that the aneurysms are usually not convex (see e.g.Fig. 6(d)).Therefore, the 2D example in Fig. 6 is a rather typical example for our applications.Fig. 6(a,b) shows that if we initialize using the blue curve, we will converge to the red curve, which is not a satisfactory result because we lost the sharp feature (the small bump).To further explain the reason we obtain such a result, we recall the evolution PDE corresponding to the energy ( 4) where As we can see from (9), the force field −∇(A(ψ)d) is indeed enhanced by A(ψ), which means around regions with large curvature we have a stronger force pushing the zero level set of φ towards the sharp tip (see Fig. 6(b)).However, because of the concavity change, the force vectors are almost tangential to the blue curve and hence most of the forces are wasted.Meanwhile, the shrinking force given by the second term of ( 9) is also enhanced by the factor A(ψ). Therefore, the blue curve will eventually shrink to the red one, instead of moving forward and capturing the entire small bleb, which is a very important feature for aneurysms (see Fig. 6(d)).To overcome this, one may choose the initial curve containing the entire small bump.From a practical point of view, however, it is not reasonable to assume that one can always start with a good initial surface, and it is always desirable to have a method with less restrictive initialization constraints.
where μ is a constant parameter and κ + (ψ) is the positive part of the Gaussian curvature of ψ.
1.The major difference of our model ( 10) from ( 5) is that, instead of putting the factor A(ψ) in the energy as in (4), we modify (3) directly and only enhance the force field −∇d.A 2D result using (10) is given in Fig. 6(a,c) showing improvement of our model over (4).More numerical experiments will be given in Section 3.4 which show significant improvement of the results using our model (10) for the blood vessel shown in Fig. 2.More details are discussed in Section 3.4.2. The choice of the positive component (instead of some other choices such as the absolute value) of the Gaussian curvature is to ensure that the resulting surface does not contain any part of the vessels.Indeed, assuming that the initial surface contains part of the blood vessels, and if the vessel locally looks like a cylinder, then its Gaussian curvatures are small, and the part of the surface on the vessel area will shrink and disappear eventually.More often than not, vessels are curved instead of straight as cylinders, as shown in the left figure of Figure 7. Since in (11) we do not enhance the force field at the region with negative Gaussian curvatures, then the part of the surface on those regions of the vessel will be peeled off quickly, and eventually all the surface parts within vessel regions will shrink and disappear.On the other hand, if the curved vessel is small in diameter, the mean curvature term β∇ • ∇φ |∇φ| dominates A(ψ)∇d • ∇φ |∇φ| , and the zero level set of φ in these regions also shrink.To illustrate the above observations visually, we show in Figure 7 the process of evolution when solving (10).The left figure of Figure 7 superimposes Gaussian curvature of the surface onto the surface itself, where the red curves specify the regions where the blood vessels are curved, with one side having positive curvature and the other side negative curvatures.The figures from the second to the last show the evolutions of φ.As one can see, all the vessels that were included in the initial guess of φ disappeared at the final stage.3. We also note that our model ( 10) can be used in other types of surface capturing problems.We just need to fashion the factor A(ψ) according to the type of surfaces and applications we have.4) is less restrictive in terms of initialization than the models (3) and ( 5), properly chosen initial guesses are still desirable.To obtain a reasonable initial surface, we adopt a user interactive strategy to initiate the computation.We let users select points around the area of interest and use the selected points to determine a sphere/ellipsoid with level set function φ s .Then φ(x, 0) is defined as the intersection of φ s with ψ, or mathematically φ(x, 0) = min{φ s (x), ψ(x)}.In our proposed method, the selection of the points for the region of interest is the only part that needs user interaction.Although automated computation is desirable, determining a pathologic region is a medical diagnosis which needs an expert's supervision.Therefore, it is reasonable to have experts' inputs and use them to initiate the computation.In the numerical section, we will show numerous clinical examples of allowing the user to select only a few (no more than six) points to capture the surfaces of brain aneurysms.
With the initial condition described above, we employ the local level set method and finite difference discretizations [22] to solve equation ( 3) and (10), as well as to minimize (4), in order to alleviate the time step restrictions and lower the complexity of our numerical computations.Generally speaking, we solve There are three parameters in our model (10), i.e. μ, α and β.The parameter μ controls the amount of force one wishes to apply near the regions with sharp features.The term αH(ψ) prevents the zero level set of φ from passing through that of ψ.Since we initialize our φ within ψ, this term only acts as a barricade and we shall fix its value through out our experiments.The parameter β controls the global smoothness of φ.The larger β is, the smoother our final results will be.In Section 4, we shall fix all the parameters μ, α and β, and the stopping tolerance.Our numerical results show that the set of parameters we have chosen gives consistently good results for all of the ten different subjects we tested.
3.4.Models Comparison.The algorithms (3), (5) and our method (10) are applied to a set of brain images acquired by 3D CT angiography.The images have 512 × 512 in-plane spatial resolutions with each voxel size approximately 0.125mm 3 .We then extract subimages of size 54 × 37 for the aneurysm from the entire brain images.The reconstruction of the surface is shown in Fig. 2. The initial surface, i.e. the zero level set of φ(x, 0), is visualized in Fig. 8 (right).
The numerical results of solving (3) are shown in Fig. 9, top row.Although this model has been reported with fairly good results for 2D images [11], direct application to capturing 3D surface is not satisfactory, as discussed in Section 3.1.
Here we also tested the model (5) which was developed to improve the illusory contours at corners [11].The results are shown in the second row of Fig. 9.This model provides improvement at the tip of the aneurysm in comparison with the model (3); however, it still can not capture the entire tip which is a very important medical feature.The reason is the concavity change near the sharp tip, which is consistent with our earlier discussion in Section 3.1.The results of our surface capturing model (10) are shown in the third row of Fig. 9. Using our method (10), we are able to capture the entire aneurysm.4) and ( 10) respectively.For the visually best results, parameters are β = 1 for (3), (μ, β)=(500, 0.05) for ( 4) and (μ, β)=(2700, 0.05) for (10).In each row, the five figures are results at iteration=0, 100, 500, 1000 and 2000 respectively.4. Validations of Our Method on Various Brain Aneurysm Data.In this section, we further validate our model (4) on ten different brain aneurysms.Brain aneurysms are classified as the narrow-neck aneurysms or wide-neck aneurysms by their dome/neck ratios.A narrow-neck aneurysm has a dome/neck ratio more than 1.5; otherwise, it is a wide-neck aneurysm [24].We shall test the consistency and robustness of our method on both types(five subjects for each category).Throughout the numerical experiment, parameters μ, α and β are taken as α = 0.5, β = 0.001, and μ = 2/mean(κ + (ψ)).Only the stopping criteria are different depending on whether the aneurysm is classified as narrow-neck or wide-neck.Numerical results presented in Fig. 10 to Fig. 13 show that this choice of parameters and stopping criteria gives consistently good results.
All the numerical experiments were performed using MATLAB on a Windows Laptop (Duo processor, 2.0GHz CPU and 2GB RAM).It took approximately one minute to capture an aneurysm with volume 100mm 3 , and an additional one minute for every 100mm 3 increase in size.
4.1.Narrow-Necked Aneurysms.We test our model (10) on five narrow-neck aneurysms data.The reconstructed surfaces from 3D images are given in the top row of Fig. 10.We initialize our computation and perform the calculation by the algorithm in Section 3.2 and 3.3, and employ the following stopping criteria for narrow-necked aneurysms where n is the iteration number which comes from the discretization of time variable t.Fig. 4.2.Wide-Necked Aneurysms.For wide-neck aneurysms, using (13) as stopping criteria may cause the zero level set of φ shrink to zero.This is because in general there does not exit a stable state for the zero level set of φ due to the fact that the necks of the aneurysms are usually wide open.Thus, we adopt the following stopping criteria based on the special geometry of wide-neck aneurysms, The above equation means that the computation stops whenever the change of φ n picks up some constant pace.We test our model (10) with the above stopping criteria ( 14) on five wide-neck aneurysms data.The reconstructed surfaces from 3D images are given in Fig. 12 top row.Fig. 12 bottom row shows the numerical results of aneurysm capturing.To test the robustness, we also randomly choose 6 different sets of initial points on one of five aneurysms in Fig. 12, which generates 6 different initial surfaces.The final results from the 6 different initializations are also nearly identical as one can see in Fig. 13 bottom row.The volumes captured by different initial points are 122.42 ± 5.37mm 3 (mean ±standard deviation).As a result, we expect the deviation of the volume computation which can be caused by different users is approximately 4.4% of the total aneurysm volume.Note that, although we do not have a theoretical guarantee that ( 14) works for all wide-neck aneurysms, we believe it should work for typical wide-necks and it certainly works well for the 5 subjects we tested on in Section 4.2.

Conclusion.
A method to quantify the volume and other geometries of brain aneurysms is needed to better study how they associate with aneurysmal growth and rupture.In this paper, we introduced a level set based PDE model to capture brain aneurysms.We also introduced a supervised strategy where users only need to select no more than six initial points on the surface to initialize the algorithm.The numerical results showed that the final surface captured the entire target region and we were able to compute the volume and curvatures of the aneurysms for clinical studies.There is huge variation among brain aneurysms and being able to quantify the geometry of irregular shapes is especially important to the study of the associations of shape with rupture.Our future work will involve applying this algorithm to diverse aneurysm shapes and adjusting the algorithm for different clinical purposes.

Figure 2 .
Figure 2. Reconstruction from 3D images.Left figure shows a few slices of images corresponding to the reconstructed surface; right figure shows the reconstructed surface with the red curve circling the aneurysm.

Figure 3 .
Figure 3. Left figure illustrates the special example consider above, while the right one shows the plot of E 1 (r).

Figure 6 .
Figure 6.Figure (a): special example with blue curve being the initial curve, red one being the result obtained by (4) and green one being the result obtained by our proposed model (10); Figure (b): one zoom-in with the blue arrows specifying the vector field −∇(A(ψ)d); Figure (c): one zoom-in with the blue arrows specifying the vector field −A(ψ)∇d; Figure (d): the blood vessel (same object as in Fig. 2 with another view) with concavity change on the aneurysm region (within the red circle).

Figure 7 .
Figure 7. Figures from upper left to lower right are: Gaussian curvature on the surface; initial surface (red) superimposed with the reference surface; initial surface; iteration 200; iteration 500; iteration 1000; and final result.
the normal velocity depending on φ (e.g.V n = −κ for mean curvature motion), and the restriction function c a,b introduced to confine all effective calculations within a narrow band of zero level set of φ.The restriction function c a,b is defined as

Figure 8 .
Figure 8. Left: selected points on the target vessel; right: initial surface (red).

Figure 10 .
Figure 10.Top row shows the surfaces of narrow-necked aneurysms.Second row shows the sets of points given by users.Third row is the corresponding initial surfaces.Bottom row is the corresponding final captured surfaces.The surfaces in row 2-4 are shown with close-up views.The volumes of the aneurysms captured are 213.527mm 3 , 520.196mm 3 , 602.7mm 3 , 319.296mm 3 and 516.399mm 3 respectively from left to right.

Figure 11 .
Figure 11.Top row is the set of points given by users.Middle row is the corresponding initial surfaces.Bottom is the corresponding final captured surfaces.The resulting volumes for the different points chosen by users from left to right are: 319.296mm 3 , 317.275mm 3 , 307.781mm 3 , 302.881mm 3 , 315.499mm 3 and 310.905mm 3 .

Figure 12 .
Figure 12.Top row shows the surfaces of narrow-necked.Second row shows the sets of points given by users.Third row is the corresponding initial surfaces.Bottom row is the corresponding final captured surfaces.The surfaces in row 2-4 are shown with close-up views.The volumes of the aneurysms captured are 78.767mm 3 , 95.823mm 3 , 117.355mm 3 , 300.493mm 3 and 748.23mm 3 respectively from left to right.

Figure 13 .
Figure 13.Top row is the set of points given by users.Middle row is the corresponding initial surfaces.Bottom is the corresponding final captured surfaces.The resulting volumes for the different points chosen by users from left to right are: 117.355mm 3 , 122.133mm 3 , 131.136mm 3 , 122.5mm 3 , 124.95mm 3 and 116.436mm 3 .
10 (bottom row) shows the numerical results of aneurysm capturing for the five subjects.The robustness of the numerical solutions is also tested by randomly choosing 6 different sets of initial points (Fig 11 top row) on one of five aneurysms in Fig. 10, which generates 6 different initial surfaces (Fig 11 middle row).The final results from the 6 different initializations are nearly identical to each other as shown in Fig. 11 bottom row.The volumes captured by different initial points are 312.273±6.245mm 3 (mean ±standard deviation).As a result, we expect the deviation of the volume computation which can be caused by different users is approximately 2% of the total aneurysm volume, which can be considered well acceptable in practice.