Piecewise Constant Level Set Method for 3D image segmentation
Piecewise Constant Level Set Method for 3D image segmentation Are Losneg˚ ard, Oddvar Christiansen, Xue-Cheng Tai Department of Mathematics, University of Bergen, Johannes Brunsgate 12, 5008 Bergen, Norway
Abstract. Level set methods have been proven to be efficient tools for tracing interface problems. Recently, some variants of the Osher- Sethian level set methods, which are called the Piecewise Constant Level Set Methods (PCLSM), have been proposed for some interface problems. The methods need to minimize a smooth cost functional under some special constraints. In this paper a PCLSM for 3D image segmentation is tested. The algorithm uses the gradient descent method combined with a Quasi-Newton method to minimize an augmented Lagrangian functional. Experiments for medical image segmentation are shown on synthetic three dimensional MR brain images. The efficiency of the algorithm and the quality of the obtained images are demonstrated.
1
Piecewise constant level set methods
In many applications one wants to divide an image into subsections based on intensity values, e.g. extract the gray matter, white matter and cerebrospinal fluid from a brain MRI, recognize the letters on a car plate or isolate any interesting objects in an image. There any many ways to achieve this segmentation. In the last decade, methods based on level sets have been quite popular [1] [2] [3] [4] [5] In the traditional level set methods [6] [7] the idea is to define a function φ(x), whose zero level set represents an interface Γ . This will divide a domain Ω into 2 subdomains, i.e. φ(x) > 0, φ(x) = 0, φ(x) < 0,
x inside Γ , x on Γ , x outside Γ .
Extensions of these methods results in the possibilities of separating the domain Ω into 2N disjunct regions, using N level set functions. Recently, methods based on Piecewise Constant Level Set Methods (PCLSM) have been introduced for multiphase segmentation [8] [9] [10]. Instead of using using N level set functions one uses a single piecewise constant level set function φ, defined as φ = i in Ωi , i = 1, 2, . . . , N . (1)
The discontinuities of φ give the curves that separate the regions. This extend the traditonal multiphase level set method of [11,12]. In [13], the layers for the φ functions are used to identify the phases. Another approach for multiphase segmentation is related to the phase-field models and the so-called binary level set methods, see [14] [15], [16], [17]. In this work we will extend the approach introduced in [9], where a novel method for multiphase image segmentation is proposed and tested on 2d images. In this work we extend the code to 3d images and test the method on synthetic MRI data. We shall note that the level set idea has been applied [18] to the geodesic active surface model for gray matter segmentation. In [19] they combine the Chan-Vese [7] model with a variant of the Haralick/Canny edge detector and the geodesic active surface model for thin structures segmentation in volumetric data. This work is organized in the following way. We start by introducing the minimization problem in section 2. In section 3 we show how the minimization problem can be solved using the gradient descent method and a quasi-Newton method, and in section 4 we show numerical results for the method.
2
The minimization problem
To segment an image u0 based on intensity values, Mumford and Shah [20] proposed to approximate the function u0 by a piecewise constant function u by solving: ( ) XZ MS 2 inf F (u, Γ ) = |ci − u0 | dx + ν|Γ | , (2) u,Γ
i
Ωi
where Ω = ∪i Ωi ∪Γ , i.e. they tried to find a decomposition Ωi of Ω, where u = ci (constant) inside each connected component Ωi . The length of the curve Γ is controlled by a positive parameter ν. For a fixed Γ , we see that (2) is minimized when ci = mean(u0 ) in Ωi . One challenge when solving (2) is to find a unique representation of the parametrized curve Γ . In [8] they propose to solve the segmentation problem by using a PCLSM, thus solving the following constrained minimization problem:
min
c,φ K(φ)=0
Z Z N Z o n X 1 |u − u0 |2 dx + β |∇ψi | dx + ν |∇φ| dx , (3) F (c, φ) = 2 Ω Ω i=1 Ω
in order to segment an image u0 into N phases. Above, The functions ψi are defined as N N Y 1 Y ψi = (φ − j), αi = (i − k) (4) αi j=1,j6=i
k=1,k6=i
and the constraint K is defined as K(φ) = (φ − 1)(φ − 2) · · · (φ − N ) =
N Y φ=i
(φ − i).
(5)
The piecewise constant image u is a linear combination of the characteristic functions N X u= ci ψi . (6) i=0
large approximation errors will be penalized by the fidelity term R We see that 2 |u − u | and the last two terms surpress oscillation, whereas the regular0 Ω ization parameters β > 0, ν > 0 control the effect of the two terms. 1 2
3
Steepest descent and Quasi-Newton
In [8] the augmented Lagrangian method was used to solve the constrained minimization problem (3), and they defined this as Z Z r |K(φ)|2 dx, (7) L(c, φ, λ) = F (c, φ) + λK(φ) dx + 2 Ω Ω where λ ∈ L2 (Ω) is the Langrange-multiplier and r > 0 is a penalty parameter. We normally choose the penalization parameter r large, compared to the other parameters. To minimize (3), we have to find the saddle points for L. The saddle points are found by minimizing L with respect to φ and c, and maximizing with respect to λ. By minimizing with respect to φ and c, we ensure that F (c, φ) is minimized, and by maximizing with respect to λ, the constraint must be fulfilled at convergence, otherwise the Lagrangian term of (7) will not vanish. The result is the the following algorithm: Algorithm 1. Choose initial values for φ0 , c0 and λ0 . For k = 1, 2, . . . , do: 1. Find ck from L(ck , φk−1 , λk−1 ) = min L(c, φk−1 , λk−1 ). c
PN 2. Use (6) to update u = i=1 cki ψi (φk−1 ). 3. Find φk from L(ck , φk , λk−1 ) = min L(ck , φ, λk−1 ). φ
(8)
(9)
PN 4. Use (6) to update u = i=1 cki ψi (φk ). 5. Update the Langrange-multiplier by λk = λk−1 + rK(φk ).
(10)
The minimizer for (9) is solved here by the gradient descent method [8]. We introduce an artificial time variable and look for a steady-state solution to the PDE ∂L φt = − . (11) ∂φ
We use the following numerical approximation to solve this problem: φnew = φold − ∆t
∂L (c, φold , λk−1 ). ∂φ
(12)
In order to compute this, we need the Gateux derivative: N X ∂L ∂u ∇φ ∇ψi ∂ψi = (u−u0 ) −β −ν∇· +λK 0 (φ)+rK(φ)K 0 (φ) ∇· ∂φ ∂φ |∇ψ | ∂φ |∇φ| i i=1 (13) and the step size ∆t is fixed during the whole iterative procedure. As always, the derivative of L with respect to λ recovers the constraint ∂L = K(φ). ∂λ
(14)
Because u is linear with respect to the ci values, we see that L is quadratic with respect to ci . Thus the minimization problem (8) can be solved exactly. We see that Z Z ∂L ∂L ∂u = = (u − u0 )ψi dx for i = 1, 2, . . . , N. (15) ∂ci Ω ∂u ∂ci Ω Therefore, the minimizer of (8) satisfies a linear system of equations Ack = b: n Z X j=1
Ω
(ψj ψi )ckj dx
Z =
u0 ψi dx,
for i = 1, 2, . . . , N.
(16)
Ω
Tests have shown [9] that this algorithm alone converges very slowly, i.e. it would require thousands of iterations for a brain MRI to achieve convergence. The strategy is to terminate it after a certain number of iterations and then use a quasi-Newton method to make it converge fast. The reason for choosing a quasi-Newton method, is to avoid inverting a huge linear algebraic system due to the regularization terms in (3). So we define Z Z 1 Q(c, φ, λ) = |u − u0 |2 dx + λK(φ) dx. (17) 2 Ω Ω Numerical experiments also reveals that it is not necessary to use the penalization term with the quasi-Newton updatings. Thus, we also define Z L0 (c, φ, λ) = F (c, φ) + λK(φ) dx. (18) Ω
We see that L0 is equal to L if we take r = 0 in (7). Also, the functional L0 reduces to Q if we take β = ν = 0. Thus the Hessian matrix for Q is a good approximation for the Hessian matrix of L0 , using the fact that β and ν are normally very small. So we arrive at the next algorithm: Algorithm 2. Choose initial values φ0 , c0 and λ0 . For k = 1, 2, . . ., do:
1. Find ck from L(ck , φk−1 , λk−1 ) = min L(c, φk−1 , λk−1 ). c
PN 2. Update u = i=1 cki ψi (φk−1 ). 3. Find φk , λk from ! 2 2 ∂ Q ∂ Q ∂φ2 ∂φ∂λ ∂2Q 0 ∂φ∂λ
φk − φk−1 λk − λk−1
(19)
∂L0
=−
∂φ ∂L0 ∂λ
(20)
PN 4. Update u = j=1 ckj ψj φk 5. If converged, end the loop. Otherwise go to step 1. The updating of c in (19) can be done in the same way as for algorithm 1. To solve (20), we need ∂2Q = ∂φ2
∂u ∂φ
2 + (u − u0 )
∂2u + λK 00 (φ), ∂φ2
∂2Q ∂2Q = = K 0 (φ), ∂φ∂λ ∂λ∂φ
where λ = λk−1 , φ = φk−1 and u = u(ck , φk−1 ). In addition we use that 0 K(φk−1 ) and ∂L ∂φ can be found from (13) by setting r = 0.
4
(21) ∂L0 ∂λ
=
Numerical results
We have applied our methods to some synthetic 3D brain MRI, and in this section we present some of the results. First we discuss how to choose the initial values, before looking into the tested images. To find φ0 we have scaled the image, which normally takes the values between 0 and 255, to take the values between 1 and n, where n is the number of phases for the image, i.e. φ0 (x) = 1 +
u0 (x) − minx∈Ω u0 · (n − 1). maxx∈Ω u0 − minx∈Ω u0
(22)
To find initial c0 , we have used the function kmeans in MatLab and numerical tests have shown that the results are improved by not updating c in (8) and (19). Initial value λ0 is set to 0. The tests are run on a computer with a 2 Opteron 270 dualcore processor and 8 GB memory. In the following, we will look at some results of a synthetic brain image obtained from [21]. This is an image of size 217×181×181 with 20% inhomogenity and 7% noise. In the numerical experiments for this image we have achieved the best results with 1000 iterations of algorithm 1 followed by 30 iterations of algorithm 2. The parameters are r = 500, ν = 500, β = 0, dt = 1e−6 and we have used 3 phases.
Fig. 1. Comparing the synthetic image, axial view. From left to right we display: The original image, the result from our method and the result from brainweb [21]. From top to bottom we show slice nr.: 61, 76, 91 and 106.
Fig. 2. Comparing the synthetic image, sagittal view. From left to right we display: The original image, the result from our method and the result from brainweb [21]. From top to bottom we show slice nr.: 45, 75, 105 and 135.
Fig. 3. A closer look at slice 91 from figure 1 for, from top to bottom, the original image and the result from our method.
In figure 1 and 2 we try to compare the original image, the image segmented by our method and the image from brainweb [21] indicating the 3 phases. This is done by looking at some of the slices from the brain. In figure 3 we take a closer look at slice 91, comparing the original image and the result from our method. Finally, in figure 4 we present the white matter extracted from our segmented image, and in figure 5 we also show the white matter, but now sliced and with slice 91 displayed in gray. By looking at figure 1 and 2, we can see that our results are very good compared to the images in the 3rd column. By looking at the slices in the 1st row in figure 1, we can observe that there is a difference in the gray matter in the lower middle part of the images. In the 2nd row, we have the same tendency and there are also minor differences with respect to the white matter in the same part of the image. In figure 2 we can observe that there is a difference in the lower right parts of the images with respect to the white matter. Also small differences with respect to the gray matter are observable in this rather complex part of the image, probably due to noise and inhomogeneity.
5
Summary
In this paper we have extended the Piecewise Constant Level Set Methods (PCLSM) and tested it on complex 3D medical MR images. The segmentation results are very promising and we believe that it has a great potential for these kinds of applications. Manual segmentation is very time consuming and we believe that our method has a potential of becoming quicker than it is today, but the convergence naturally depends on the size and the complexity of the image. In future work we will test our method on 3D clinical MRI data which is an even harder task.
References 1. Osher, S., Fedkiw, R.: An overview and some recent results. J.Comput.Phys, 169 No. 2:463-502 (2001) 2. Tai, X.C., Chan, T.F.: A survey on multiple level set methods with applications for identifying piecewise constant functions. Int. J. Numer. Anal. Model, 1(1):25-47 (2004) 3. Osher, S., Burger, M.: A survey on level set methods for inverse problems and optimal design. Cam-report-04-02, UCLA, Applied Mathematics (2004) 4. Vese, L.A., Chan, T.F.: A new multiphase level set framework for image segmentation via the Mumford and Shah model. International Journal of Computer Vision 50 (2002.) 271–293 5. Chan, T.F., Vese, L.A.: Image segmentation using level sets and the piecewise constant mumford-shah model. Technical report, CAM Report 00-14, UCLA, Math. Depart. (April 2000) revised December 2000. 6. Osher, S., Sethian, J.: Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys., vol. 79, no. 1 (1988)
7. Chan, T., Vese, L.: Active contours without edges. IEEE Image Proc., 10, pp. 266-277 (2001) 8. Lie, J., Lysaker, M., Tai, X.C.: A variant of the level set method and applications to image segmentation. Math. Comp. 75 no. 255 (2006) 9. Tai, X.c., Yao, C.h.: Image segmentation by piecewise constant Mumford-Shah model without estimating the constants. J. Comput. Math. 24(3) (2006) 435–443 10. Christiansen, O., Tai, X.C.: Fast implementation of piecewise constant level set methods. In: Image processing based on partial differential equations, Springer Verlag (2006) 253–272 11. Zhao, H.K., Chan, T., Merriman, B., Osher, S.: A variational level set approach to multiphase motion. J. Comput. Phys. 127(1) (1996) 179–195 12. Vese, L., Chan, T.: A new multiphase level set framework for image segmentation via the Mumford and Shah model. International Journal of Computer Vision, 50, pp.271-293 (2002) 13. Chung, J., Vese, A.: Energy minimization based segmentation and denoising using a multilayer level set approach. Lecture notes in computer science 3757 (2005) 439–455 14. Song, B., Chan, T.F.: Fast algorithm for level set segmentation. UCLA CAM report 02-68 (2002) 15. Gibou, F., Fedkiw, R.: A fast hybrid k-means level set algorithm for segmentation. 4th Annual Hawaii International Conference on Statistics and Mathematics, pp. 281-291 (2005) 16. Esedo¯ glu, S., Tsai, Y.H.R.: Threshold dynamics for the piecewise constant Mumford-Shah functional. J. Comput. Phys. 211(1) (2006) 367–384 17. Shi, Y., Karl, W.C.: A fast level set method without solving pdes. In: ICASSP’05. (2005) 18. Goldenberg, R., Kimmel, R., Rivlin, E., Rudzsky, M.: Cortex Segmentation: A Fast Variational Geometric Approach. IEEE Transactions on medical imaging, Vol. 21, No. 2 (2002) 19. Holtzman-Gazit, M., Kimmel, R., Peled, N., Goldsher, D.: Segmentation of Thin Structures in Volumetric Medical Images. IEEE Transactions on image processing, Vol. 15, No. 2 (2006) 20. Mumford, D., Shah, J.: Optimal approximation by piecewise smooth functions and associated variational problems. Comm. Pure Appl. Math, 42 (1989) 21. McConnell Brain Imaging Centre, Montr´eal Neurological Institute, M.U.: Brainweb. http://www.bic.mni.mcgill.ca/brainweb/ (12.6.2006)
Fig. 4. We have extracted the white matter from our test image and used the built-in MatLab function bwareaopen to strip the skull together with a little gaussian smoothing.
Fig. 5. Lower part of the white matter, where slice 91 is displayed in gray.
Abstract. Level set methods have been proven to be efficient tools for tracing interface problems. Recently, some variants of the Osher- Sethian level set methods, which are called the Piecewise Constant Level Set Methods (PCLSM), have been proposed for some interface problems. The methods need to minimize a smooth cost functional under some special constraints. In this paper a PCLSM for 3D image segmentation is tested. The algorithm uses the gradient descent method combined with a Quasi-Newton method to minimize an augmented Lagrangian functional. Experiments for medical image segmentation are shown on synthetic three dimensional MR brain images. The efficiency of the algorithm and the quality of the obtained images are demonstrated.
1
Piecewise constant level set methods
In many applications one wants to divide an image into subsections based on intensity values, e.g. extract the gray matter, white matter and cerebrospinal fluid from a brain MRI, recognize the letters on a car plate or isolate any interesting objects in an image. There any many ways to achieve this segmentation. In the last decade, methods based on level sets have been quite popular [1] [2] [3] [4] [5] In the traditional level set methods [6] [7] the idea is to define a function φ(x), whose zero level set represents an interface Γ . This will divide a domain Ω into 2 subdomains, i.e. φ(x) > 0, φ(x) = 0, φ(x) < 0,
x inside Γ , x on Γ , x outside Γ .
Extensions of these methods results in the possibilities of separating the domain Ω into 2N disjunct regions, using N level set functions. Recently, methods based on Piecewise Constant Level Set Methods (PCLSM) have been introduced for multiphase segmentation [8] [9] [10]. Instead of using using N level set functions one uses a single piecewise constant level set function φ, defined as φ = i in Ωi , i = 1, 2, . . . , N . (1)
The discontinuities of φ give the curves that separate the regions. This extend the traditonal multiphase level set method of [11,12]. In [13], the layers for the φ functions are used to identify the phases. Another approach for multiphase segmentation is related to the phase-field models and the so-called binary level set methods, see [14] [15], [16], [17]. In this work we will extend the approach introduced in [9], where a novel method for multiphase image segmentation is proposed and tested on 2d images. In this work we extend the code to 3d images and test the method on synthetic MRI data. We shall note that the level set idea has been applied [18] to the geodesic active surface model for gray matter segmentation. In [19] they combine the Chan-Vese [7] model with a variant of the Haralick/Canny edge detector and the geodesic active surface model for thin structures segmentation in volumetric data. This work is organized in the following way. We start by introducing the minimization problem in section 2. In section 3 we show how the minimization problem can be solved using the gradient descent method and a quasi-Newton method, and in section 4 we show numerical results for the method.
2
The minimization problem
To segment an image u0 based on intensity values, Mumford and Shah [20] proposed to approximate the function u0 by a piecewise constant function u by solving: ( ) XZ MS 2 inf F (u, Γ ) = |ci − u0 | dx + ν|Γ | , (2) u,Γ
i
Ωi
where Ω = ∪i Ωi ∪Γ , i.e. they tried to find a decomposition Ωi of Ω, where u = ci (constant) inside each connected component Ωi . The length of the curve Γ is controlled by a positive parameter ν. For a fixed Γ , we see that (2) is minimized when ci = mean(u0 ) in Ωi . One challenge when solving (2) is to find a unique representation of the parametrized curve Γ . In [8] they propose to solve the segmentation problem by using a PCLSM, thus solving the following constrained minimization problem:
min
c,φ K(φ)=0
Z Z N Z o n X 1 |u − u0 |2 dx + β |∇ψi | dx + ν |∇φ| dx , (3) F (c, φ) = 2 Ω Ω i=1 Ω
in order to segment an image u0 into N phases. Above, The functions ψi are defined as N N Y 1 Y ψi = (φ − j), αi = (i − k) (4) αi j=1,j6=i
k=1,k6=i
and the constraint K is defined as K(φ) = (φ − 1)(φ − 2) · · · (φ − N ) =
N Y φ=i
(φ − i).
(5)
The piecewise constant image u is a linear combination of the characteristic functions N X u= ci ψi . (6) i=0
large approximation errors will be penalized by the fidelity term R We see that 2 |u − u | and the last two terms surpress oscillation, whereas the regular0 Ω ization parameters β > 0, ν > 0 control the effect of the two terms. 1 2
3
Steepest descent and Quasi-Newton
In [8] the augmented Lagrangian method was used to solve the constrained minimization problem (3), and they defined this as Z Z r |K(φ)|2 dx, (7) L(c, φ, λ) = F (c, φ) + λK(φ) dx + 2 Ω Ω where λ ∈ L2 (Ω) is the Langrange-multiplier and r > 0 is a penalty parameter. We normally choose the penalization parameter r large, compared to the other parameters. To minimize (3), we have to find the saddle points for L. The saddle points are found by minimizing L with respect to φ and c, and maximizing with respect to λ. By minimizing with respect to φ and c, we ensure that F (c, φ) is minimized, and by maximizing with respect to λ, the constraint must be fulfilled at convergence, otherwise the Lagrangian term of (7) will not vanish. The result is the the following algorithm: Algorithm 1. Choose initial values for φ0 , c0 and λ0 . For k = 1, 2, . . . , do: 1. Find ck from L(ck , φk−1 , λk−1 ) = min L(c, φk−1 , λk−1 ). c
PN 2. Use (6) to update u = i=1 cki ψi (φk−1 ). 3. Find φk from L(ck , φk , λk−1 ) = min L(ck , φ, λk−1 ). φ
(8)
(9)
PN 4. Use (6) to update u = i=1 cki ψi (φk ). 5. Update the Langrange-multiplier by λk = λk−1 + rK(φk ).
(10)
The minimizer for (9) is solved here by the gradient descent method [8]. We introduce an artificial time variable and look for a steady-state solution to the PDE ∂L φt = − . (11) ∂φ
We use the following numerical approximation to solve this problem: φnew = φold − ∆t
∂L (c, φold , λk−1 ). ∂φ
(12)
In order to compute this, we need the Gateux derivative: N X ∂L ∂u ∇φ ∇ψi ∂ψi = (u−u0 ) −β −ν∇· +λK 0 (φ)+rK(φ)K 0 (φ) ∇· ∂φ ∂φ |∇ψ | ∂φ |∇φ| i i=1 (13) and the step size ∆t is fixed during the whole iterative procedure. As always, the derivative of L with respect to λ recovers the constraint ∂L = K(φ). ∂λ
(14)
Because u is linear with respect to the ci values, we see that L is quadratic with respect to ci . Thus the minimization problem (8) can be solved exactly. We see that Z Z ∂L ∂L ∂u = = (u − u0 )ψi dx for i = 1, 2, . . . , N. (15) ∂ci Ω ∂u ∂ci Ω Therefore, the minimizer of (8) satisfies a linear system of equations Ack = b: n Z X j=1
Ω
(ψj ψi )ckj dx
Z =
u0 ψi dx,
for i = 1, 2, . . . , N.
(16)
Ω
Tests have shown [9] that this algorithm alone converges very slowly, i.e. it would require thousands of iterations for a brain MRI to achieve convergence. The strategy is to terminate it after a certain number of iterations and then use a quasi-Newton method to make it converge fast. The reason for choosing a quasi-Newton method, is to avoid inverting a huge linear algebraic system due to the regularization terms in (3). So we define Z Z 1 Q(c, φ, λ) = |u − u0 |2 dx + λK(φ) dx. (17) 2 Ω Ω Numerical experiments also reveals that it is not necessary to use the penalization term with the quasi-Newton updatings. Thus, we also define Z L0 (c, φ, λ) = F (c, φ) + λK(φ) dx. (18) Ω
We see that L0 is equal to L if we take r = 0 in (7). Also, the functional L0 reduces to Q if we take β = ν = 0. Thus the Hessian matrix for Q is a good approximation for the Hessian matrix of L0 , using the fact that β and ν are normally very small. So we arrive at the next algorithm: Algorithm 2. Choose initial values φ0 , c0 and λ0 . For k = 1, 2, . . ., do:
1. Find ck from L(ck , φk−1 , λk−1 ) = min L(c, φk−1 , λk−1 ). c
PN 2. Update u = i=1 cki ψi (φk−1 ). 3. Find φk , λk from ! 2 2 ∂ Q ∂ Q ∂φ2 ∂φ∂λ ∂2Q 0 ∂φ∂λ
φk − φk−1 λk − λk−1
(19)
∂L0
=−
∂φ ∂L0 ∂λ
(20)
PN 4. Update u = j=1 ckj ψj φk 5. If converged, end the loop. Otherwise go to step 1. The updating of c in (19) can be done in the same way as for algorithm 1. To solve (20), we need ∂2Q = ∂φ2
∂u ∂φ
2 + (u − u0 )
∂2u + λK 00 (φ), ∂φ2
∂2Q ∂2Q = = K 0 (φ), ∂φ∂λ ∂λ∂φ
where λ = λk−1 , φ = φk−1 and u = u(ck , φk−1 ). In addition we use that 0 K(φk−1 ) and ∂L ∂φ can be found from (13) by setting r = 0.
4
(21) ∂L0 ∂λ
=
Numerical results
We have applied our methods to some synthetic 3D brain MRI, and in this section we present some of the results. First we discuss how to choose the initial values, before looking into the tested images. To find φ0 we have scaled the image, which normally takes the values between 0 and 255, to take the values between 1 and n, where n is the number of phases for the image, i.e. φ0 (x) = 1 +
u0 (x) − minx∈Ω u0 · (n − 1). maxx∈Ω u0 − minx∈Ω u0
(22)
To find initial c0 , we have used the function kmeans in MatLab and numerical tests have shown that the results are improved by not updating c in (8) and (19). Initial value λ0 is set to 0. The tests are run on a computer with a 2 Opteron 270 dualcore processor and 8 GB memory. In the following, we will look at some results of a synthetic brain image obtained from [21]. This is an image of size 217×181×181 with 20% inhomogenity and 7% noise. In the numerical experiments for this image we have achieved the best results with 1000 iterations of algorithm 1 followed by 30 iterations of algorithm 2. The parameters are r = 500, ν = 500, β = 0, dt = 1e−6 and we have used 3 phases.
Fig. 1. Comparing the synthetic image, axial view. From left to right we display: The original image, the result from our method and the result from brainweb [21]. From top to bottom we show slice nr.: 61, 76, 91 and 106.
Fig. 2. Comparing the synthetic image, sagittal view. From left to right we display: The original image, the result from our method and the result from brainweb [21]. From top to bottom we show slice nr.: 45, 75, 105 and 135.
Fig. 3. A closer look at slice 91 from figure 1 for, from top to bottom, the original image and the result from our method.
In figure 1 and 2 we try to compare the original image, the image segmented by our method and the image from brainweb [21] indicating the 3 phases. This is done by looking at some of the slices from the brain. In figure 3 we take a closer look at slice 91, comparing the original image and the result from our method. Finally, in figure 4 we present the white matter extracted from our segmented image, and in figure 5 we also show the white matter, but now sliced and with slice 91 displayed in gray. By looking at figure 1 and 2, we can see that our results are very good compared to the images in the 3rd column. By looking at the slices in the 1st row in figure 1, we can observe that there is a difference in the gray matter in the lower middle part of the images. In the 2nd row, we have the same tendency and there are also minor differences with respect to the white matter in the same part of the image. In figure 2 we can observe that there is a difference in the lower right parts of the images with respect to the white matter. Also small differences with respect to the gray matter are observable in this rather complex part of the image, probably due to noise and inhomogeneity.
5
Summary
In this paper we have extended the Piecewise Constant Level Set Methods (PCLSM) and tested it on complex 3D medical MR images. The segmentation results are very promising and we believe that it has a great potential for these kinds of applications. Manual segmentation is very time consuming and we believe that our method has a potential of becoming quicker than it is today, but the convergence naturally depends on the size and the complexity of the image. In future work we will test our method on 3D clinical MRI data which is an even harder task.
References 1. Osher, S., Fedkiw, R.: An overview and some recent results. J.Comput.Phys, 169 No. 2:463-502 (2001) 2. Tai, X.C., Chan, T.F.: A survey on multiple level set methods with applications for identifying piecewise constant functions. Int. J. Numer. Anal. Model, 1(1):25-47 (2004) 3. Osher, S., Burger, M.: A survey on level set methods for inverse problems and optimal design. Cam-report-04-02, UCLA, Applied Mathematics (2004) 4. Vese, L.A., Chan, T.F.: A new multiphase level set framework for image segmentation via the Mumford and Shah model. International Journal of Computer Vision 50 (2002.) 271–293 5. Chan, T.F., Vese, L.A.: Image segmentation using level sets and the piecewise constant mumford-shah model. Technical report, CAM Report 00-14, UCLA, Math. Depart. (April 2000) revised December 2000. 6. Osher, S., Sethian, J.: Fronts propagating with curvature dependent speed: algorithms based on Hamilton-Jacobi formulations. J. Comput. Phys., vol. 79, no. 1 (1988)
7. Chan, T., Vese, L.: Active contours without edges. IEEE Image Proc., 10, pp. 266-277 (2001) 8. Lie, J., Lysaker, M., Tai, X.C.: A variant of the level set method and applications to image segmentation. Math. Comp. 75 no. 255 (2006) 9. Tai, X.c., Yao, C.h.: Image segmentation by piecewise constant Mumford-Shah model without estimating the constants. J. Comput. Math. 24(3) (2006) 435–443 10. Christiansen, O., Tai, X.C.: Fast implementation of piecewise constant level set methods. In: Image processing based on partial differential equations, Springer Verlag (2006) 253–272 11. Zhao, H.K., Chan, T., Merriman, B., Osher, S.: A variational level set approach to multiphase motion. J. Comput. Phys. 127(1) (1996) 179–195 12. Vese, L., Chan, T.: A new multiphase level set framework for image segmentation via the Mumford and Shah model. International Journal of Computer Vision, 50, pp.271-293 (2002) 13. Chung, J., Vese, A.: Energy minimization based segmentation and denoising using a multilayer level set approach. Lecture notes in computer science 3757 (2005) 439–455 14. Song, B., Chan, T.F.: Fast algorithm for level set segmentation. UCLA CAM report 02-68 (2002) 15. Gibou, F., Fedkiw, R.: A fast hybrid k-means level set algorithm for segmentation. 4th Annual Hawaii International Conference on Statistics and Mathematics, pp. 281-291 (2005) 16. Esedo¯ glu, S., Tsai, Y.H.R.: Threshold dynamics for the piecewise constant Mumford-Shah functional. J. Comput. Phys. 211(1) (2006) 367–384 17. Shi, Y., Karl, W.C.: A fast level set method without solving pdes. In: ICASSP’05. (2005) 18. Goldenberg, R., Kimmel, R., Rivlin, E., Rudzsky, M.: Cortex Segmentation: A Fast Variational Geometric Approach. IEEE Transactions on medical imaging, Vol. 21, No. 2 (2002) 19. Holtzman-Gazit, M., Kimmel, R., Peled, N., Goldsher, D.: Segmentation of Thin Structures in Volumetric Medical Images. IEEE Transactions on image processing, Vol. 15, No. 2 (2006) 20. Mumford, D., Shah, J.: Optimal approximation by piecewise smooth functions and associated variational problems. Comm. Pure Appl. Math, 42 (1989) 21. McConnell Brain Imaging Centre, Montr´eal Neurological Institute, M.U.: Brainweb. http://www.bic.mni.mcgill.ca/brainweb/ (12.6.2006)
Fig. 4. We have extracted the white matter from our test image and used the built-in MatLab function bwareaopen to strip the skull together with a little gaussian smoothing.
Fig. 5. Lower part of the white matter, where slice 91 is displayed in gray.
