An improved variational level set method for MR image segmentation and bias field correction

An improved variational level set method for MR image segmentation and bias field correction

Magnetic Resonance Imaging 31 (2013) 439–447 Contents lists available at SciVerse ScienceDirect Magnetic Resonance Imaging journal homepage: www.mri...

2MB Sizes 7 Downloads 226 Views

Magnetic Resonance Imaging 31 (2013) 439–447

Contents lists available at SciVerse ScienceDirect

Magnetic Resonance Imaging journal homepage: www.mrijournal.com

An improved variational level set method for MR image segmentation and bias field correction Tianming Zhana, Jun Zhangb, Liang Xiaoa, Yunjie Chenc, Zhihui Weia, b,⁎ a b c

School of Computer Science and Technology, Nanjing University of Science and Technology, Nanjing, 210094, China School of Science, Nanjing University of Science and Technology, Nanjing, 210094, China School of Math and Phy, Nanjing University of Information Science and Technology, Nanjing, 210044, China

a r t i c l e

i n f o

Article history: Received 16 November 2010 Revised 11 June 2012 Accepted 22 August 2012 Keywords: Medical image segmentation Bias field correction Three-phase level set framework

a b s t r a c t In this paper, we propose an improved variational level set approach to correct the bias and to segment the magnetic resonance (MR) images with inhomogeneous intensity. First, we use a Gaussian distribution with bias field as a local region descriptor in two-phase level set formulation for segmentation and bias field correction of the images with inhomogeneous intensities. By using the information of the local variance in this descriptor, our method is able to obtain accurate segmentation results. Furthermore, we extend this method to three-phase level set formulation for brain MR image segmentation and bias field correction. By using this three-phase level set function to replace the four-phase level set function, we can reduce the number of convolution operations in each iteration and improve the efficiency. Compared with other approaches, this algorithm demonstrates a superior performance. © 2013 Elsevier Inc. All rights reserved.

1. Introduction Automatic segmentation plays a vital role in numerous biomedical imaging applications. Segmentation is a method to delineate the anatomical structures and other regions of interest. It has many applications in the medical field such as the quantifying of tissue volumes [1], diagnosing [2], localization of pathology [3], studying anatomical structure [4], and planning the treatment [5]. However, due to the limitations in imaging instruments and other external effects, intensity inhomogeneities often occur in the image which can be ascribed to a spatially varying field (generally called bias field). The traditional intensity-based models [5–11] which assume that the intensities are homogeneous in each subregion cannot obtain accurate segmentation results. Therefore, intensity correction is often a necessary step before MR image segmentation. In the past two decades, bias field correction has been extensively studied. The most popular methods for bias field correction are segmentation based approaches [12–15]. In these methods, segmentation and bias field estimation are interleaved in an iterative process to benefit each other, thereby allowing their optimal solutions to be simultaneously achieved. In this process, the segmentation is usually achieved by using maximum likelihood or maximum a posteriori based methods [12,13] and Fuzzy-C-Means based methods [14]. In [15], Li et al. proposed a variational level set method to estimate the bias field and to segment the images with ⁎ Corresponding author. E-mail address: [email protected] (Z. Wei). 0730-725X/$ – see front matter © 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.mri.2012.08.002

inhomogeneous intensity. A unique feature of this method is that the smoothness of the computed bias field is intrinsically ensured by the data term in the variational formulation without any additional effort for maintaining the smoothness of the bias field. However, only using intensity means as a criterion of classification, this method may not be able to precisely extract the object boundaries, especially when the object texture distributions have identical means but different variances. In this paper, we use Gaussian distribution with bias field as a local region descriptor to improve the method in [15] for simultaneous tissue classification and bias field estimation of MR images. The ideal of this paper is similar to that in [17] which uses local Gaussian distribution to improve the Euclidean distance of the energy function in [16]. Compared with Li et al.'s method, the local variance in our method enables the proposed algorithm to distinguish regions with different intensity variances. Moreover, this method is extended to three-phase level set formulation for brain MR image segmentation and bias correction. By using this three-phase formulation, it is able to reduce the number of convolution operations in each iteration and improve the efficiency. 2. Li et al.'s method Li et al. proposed a variational level set method to estimate the bias field and to segment the images with intensity inhomogeneities based on the following model to describe images in [15]: I ¼ b⋅J þ n:

ð1Þ

440

T. Zhan et al. / Magnetic Resonance Imaging 31 (2013) 439–447

where I is the measured image intensity, b is the unknown bias field, J is the true signal to be restored, n is additive noise. Our goal is to estimate the unknown bias field b and true signal J from the measured image I. In [15], it is assumed that the true image J and the bias field b have the following properties: (P1) The bias field b is slowly varying in the entire image domain. (P2) The true image intensities J are approximately a constant within each class of issue, i.e. J(x)≈ ci for x ∈ Ωi, with {Ωi}iN=1 being a partition of Ω. Let Ox = {y : ‖x − y‖ ≤ ρ} be the circular neighborhood with a small radius ρ. Then Ωi ∩ Ox is a partition of the neighborhood Ox produced by Ωi. According to the above assumption (P1), the values b(y) for all y in the circular neighborhood Ox can be well approximated by b(x). Therefore, the intensities b(y)J(y) in each subregion Ωi ∩ Ox are approximately the constant b(x)ci. Considering the task of classifying the intensities I(y) in the neighborhood Ox into N classes, Li et al. proposed an energy function as follows: Ex ¼

N X

2

∫Ox ∩Ωi K ðx−y ÞjI ðy Þ−bðx Þci j dy;

ð2Þ

where ˜p i;x ðI ðy ÞÞ is the shortened form of ˜p ðI ðy Þjy∈i ∩Ox Þ, which is the probability density in region Ωi ∩ Ox. b(x)ci and σi(x) are the local intensity mean and standard deviation respectively. Compared to [17], the local mean ui(x) is replaced by b(x)ci in the distribution. Based on this new distribution, our descriptor can recognize the intersection of bias correction and segmentation for images with intensity inhomogeneity. According to the Bayesian formulation and maximum a posteriori probability (MAP) method [17,18], we define the following objective function for each point x in the image domain Ω: Ex ¼

N X

∫Ωi −K ðx−y Þlog˜p i;x ðI ðy ÞÞdy;

To minimize Ex for all the points x, the integral of Ex is minimized over Ω, and the level set functions are employed to represent each partition {Ωi}iN= 1. The entire objective energy function for image segmentation and bias correction can be written as follows: F ðϕ; b; u 1 ; u 2 ; 1 ; ϕ2 Þ

i¼1

where b(x)ci are the cluster centers to be optimized, and K(x − y) is a non-negative weighting function so that K(x − y) = 0 for |x − y| N ρ and ∫Ox K ðx−y Þdy ¼ 1. By minimizing εx for all the points x in the whole image domain Ω and using level set functions to represent a partition {Ωi}iN= 1, Li et al. defined the unified variational level set formulation for bias correction and segmentation. However, only using intensity means as a criterion of classification, this method is hard to correctly segment the image when the mean values of intensities in background and object are rather close to each other. In this paper, we use a Gaussian distribution with the information of a bias field as a local region descriptor and propose an improved method for simultaneous bias correction and segmentation. With more statistical information in our new region descriptor, our method can obtain more accurate results than Li et al.'s method. In Section 4, our energy function is extended to three-phase level set formulation for brain MR image segmentation and bias field correction. By using this three-phase level set formulation, the number of convolution operations in each iteration is reduced to improve the efficiency of segmentation and bias field correction.

  ¼ −∫ ∫K ðx−y Þlog˜p 1;x ðI ðy ÞÞH ε ððy ÞÞdy dx   −∫ ∫K ðx−y Þlog˜p 2;x ðI ðy ÞÞð1−H ε ðϕðy ÞÞÞdy dx þv∫j∇H ε ðϕðx ÞÞjdx þ μ∫pðj∇ϕjÞdx; ð6Þ

where the third term is the length term to smooth the contour and the last term is the level set regularization term to penalize the deviation of the level set function ϕ from a signed distance function [19]. H ðx Þ ¼ 12½1 þ 2arctan ðx Þ is an approximate smoothing function of p the ffiffiffi Heaviside function. Obviously, if we set σi(x) in Eq. (6) to 1= 2, our method and Li et al.'s method in [15] are similar. However, containing the information of local variance in the local region descriptor, the proposed method is able to obtain more accurate results. The minimizer of F in each variable is obtained by fixing other variables. By fixing ci, σi, i = 1, 2, and b, minimization of the energy functional F in Eq. (6) with respect to ϕ can be achieved by solving the gradient descent flow equation:       ∂ϕ ∇ þ μ div d p j∇ϕj ∇ϕ ; ¼ −δε ðϕÞðe1 −e2 Þ þ δε ðϕÞdiv j∇ϕj ∂t



u i ðx Þ ¼

y∈Ωi ∩O x

jΩj

I ðy Þ





y∈Ωi ∩O x

bðy ÞJ ðy Þ

jΩj





y∈Ωi ∩O x

jΩj

ei ¼ ∫Ω ωðx−y Þ logðσ i ðx ÞÞ þ

ð3Þ

# ðI ðy Þ−bðx Þci Þ2 dy: 2σ i ðx Þ2

ð8Þ

By fixing ϕ, σ and c, we find an optimal bias field b that minimizes F: 2 X



bðx Þci ¼ bðx Þci ;

2

"

Gaussian distribution is employed as the descriptor of local intensities in the proposed method to extend Li et al's model in [15]. Let ui(x) be the local intensity means in the partition Ωi ∩ Ox. From (P2) in the above section, we can prove that ui(x) can be approximated by b(x)ci:

ð7Þ

ε where δε ðx Þ ¼ H ε ðx Þ ¼ π1 ε þx , div(⋅) is the divergence operator, dp is a function defined by d p ðsÞ ¼ p′ ðsÞs and e1 and e2 are the functions as follows: 2

3. Improved variational level set method based on Local Gaussian distribution

ð5Þ

i¼1

ðI ci M i ðϕÞÞ  K :=σ i

i¼1 2  X

 c2i M i ðϕÞ  K :=σ i

ð9Þ

i¼1

where |Ω| is the number of pixels in Ωi ∩ Ox. Then, similar to other statistical approaches for medical Image segmentation [12,13], the I(y) is assumed to obey the following Gaussian distribution in subregion Ωi ∩ Ox:

It is worth noting that, similar to Li et al.'s method [15], the smoothness of b in Eq. (9) is intrinsically ensured by the first two terms of our model. Using a similar optimal process, the global means ci and local variance σi2(x) can be written as:

! 1 ðbðx Þci −I ðy ÞÞ2 ˜p i;x ðI ðy ÞÞ ¼ pffiffiffiffiffiffi ; i ¼ 1; 2: exp − 2σ i ðx Þ2 2πσ i ðx Þ

∫ðb  K ÞI M i ðÞdx ; ci ¼   ∫ b2  K M i ðϕÞdx

ð4Þ

ð10Þ

T. Zhan et al. / Magnetic Resonance Imaging 31 (2013) 439–447

The curve evolution equations are given as:

and 2

2 i ðx Þ

441

¼

∫K ðy−x ÞðI ðy Þ−bðx Þci Þ M i ðϕðy ÞÞdy ∫K ðy−x ÞM i ðϕðy ÞÞdy

:

ð11Þ

To minimize the energy, the level set ϕ, bias field b, global means and local standard deviations are updated by using Eqs. (7)–(11) alternately in each iteration. After the ultimate iteration, we can compute the corrected image J with J ¼ I :=˜b, where ˜b is the optimal bias field.

For brain MR image segmentation, two level set functions are always used to create four phases which represent three tissues and background in brain the MR image respectively [7,15,20–22]. However, if using the same way to represent different brain tissues, the efficiency of our energy minimization will be much lower because of employing convolution calculation. To reduce the unnecessary convolution computations and improve the efficiency of energy minimization, a three-phase level set formulation is created in this section. In brain MR images, the intensity values in the background are close to zero. Based on this characteristic, a symbolic matrix is used to separate the background from brain tissues as follows: ðx Þ ¼

1 x∉background 0 x∈background

  ∂ϕ2 ∇ϕ2 ¼ −δðϕÞððe1 −e2 ÞH ðϕ1 ÞÞ þ δðϕÞdiv j∇2 j ∂t     þμ div d p j∇ϕ2 j ∇ϕ2 :

ð14Þ

ð15Þ

Global mean ci and local standard deviation σi can be calculated by using the following equations:

4. Three-phase level set formulation for brain MR image segmentation



  ∂ϕ1 ∇ϕ1 ht ¼ −δðϕÞððe1 −e2 ÞH ðϕ2 Þ−ðe2 −e3 ÞÞ þ δðϕÞdiv j∇ϕ1 j ∂t     þ div d p j∇ϕ1 j ∇ϕ1 ;

ð12Þ

∫ðb  K ÞI M i ðΦÞdx ; ci ¼   ∫ b2  K M i ðΦÞdx 2

σ i ðx Þ ¼

ð16Þ

∫ωðy−x ÞðI ðy Þ−bðx Þci Þ2 M i ðΦðy ÞÞðy Þdy ∫ωðy−x ÞM i ðΦðy ÞÞðy Þdy

:

ð17Þ

Similar to Eq. (10), the bias field is given as: 3 X



ðI ci M i ðϕÞÞ  K :=σ i

i¼1 3  X

 ci 2 M i ðϕÞ  K :=i

:

ð18Þ

i¼1

Then, a two level set functions and symbolic matrix χ are used to construct a three-phase level set formulation to detect the remaining three tissues (see Fig. 1). The three regions are represented by {ϕ 1 b 0, χ = 1}, {ϕ 1 N 0, ϕ 2 b 0, χ = 1}, and {ϕ 1 N 0, ϕ 2 N 0, χ = 1} respectively. Our entire energy function in three-phase level set formulation can be written as:

Unlike the energy minimization in four-phase level set formulation, it is not necessary to calculate the mean and the standard deviation of background region. Moreover, in Eqs. (14) and (15), we only calculate ei − ej twice rather than four times in four-phase framework, where ei − ej is a linear combination of many convolutions. By reducing these unnecessary convolution operations in every iteration step, this method significantly improves the efficiency of energy minimization.

F ðϕ1 ; ϕ2 ; b; c; σ Þ ¼ 3 h i X − ∫ ∫K ðx−y Þlog˜p i;x ðI ðy ÞÞM i ððy ÞÞðy Þdy dx;

5. Experimental results and Performance Analysis

ð13Þ

i¼1

5.1. Two-phase segmentation results and bias field corrections

where c = {c1,c2,c3}, σ = {σ1,σ2,σ3}, ˜p i;x ðI ðy ÞÞ are given as Eq. (4), M1(Φ) = Hε(ϕ1)Hε(ϕ2), M2(Φ) = Hε(ϕ1)(1 − Hε(ϕ2)), and M3(Φ) = 1 − Hε(ϕ1).

Fig. 1. The graphic explanation of the three-phase framework.

In this subsection, the two-phase level set formulation of our method will be tested with synthetic and real images from different modalities. Unless otherwise specified, we use the following parameters in this subsection: σ = 3.0, time step: Δt = 0.1, μ = 1.0, and ν = 0.001 × 255 × 255. We first show the experimental results on a synthetic image in Fig. 2. Fig. 2(A) shows a synthetic image with intensity inhomogeneity. Using the initial curve in Fig. 2(B), the segmentation results of Chan–Vese method, Wang et al.'s method (LGD) and our method are shown from Fig. 2(C) to Fig. 2(E), respectively. From these segmentation results, we can see that the Chan–Vese model fails to extract the object boundary, because it assumed the image consists of statistically homogeneous regions. Taking local intensity into account, both the LGD method and our method are able to detect boundary correctly. Compared to LGD method, our method also can obtain the bias corrected image shown in Fig. 2(F) beside the accurate segmentation. The improvement of the image quality can be demonstrated by comparing the histograms of the original image and the bias corrected image shown in Fig. 2(G) and Fig. 2(H). In the next experiment, we test Li et al.'s method and our method on a synthetic image with a star shaped object shown in Fig. 3(A). The object and background in this image have same mean value but different variance. The segmentation result and bias field obtained

442

T. Zhan et al. / Magnetic Resonance Imaging 31 (2013) 439–447

A

B

C

D

E

F

Histogram of original image

Histogram of bias corrected image

3000

3000

2500

2500

2000

2000

1500

1500

1000

1000

500

500

0

0

50

100

150

200

250

0

0

50

G

100

150

200

250

H

Fig. 2. Synthetic image. (A) The original image; (B) initial curve; (C) segmentation result of Chan-Vese model; (D) segmentation result of LGD method; (E) segmentation result of proposed method; (F) bias corrected image of proposed method; (G) histogram of original image; (H) histogram of corrected image.

by Li et al.'s method are shown in Fig. 3(B) and Fig. 3(C). We can see that Li et al.'s method which only use the Euclidean distance as a criterion of classification cannot detect the object boundary correctly. Fig. 3(D) and Fig. 3(E) show the segmentation result of our method. Taking more statistical characteristics of local region into account, our model is able to accurately distinguish the background and foreground and estimate a smooth bias field shown in Fig. 3(E). 5.2. Three-phase Segmentation results and bias field corrections for Brain MR Image In this subsection, we will demonstrate the effectiveness of our method in the segmentation and bias field correction of a brain MR image. In our experiments, we use the following parameters: σ = 4.0, time step: Δt = 0.1, μ = 1.0, and ν = 0.0001 × 255 × 255. To construct a three-phase framework, the symbolic matrix is set as follows: ðx Þ ¼

1 0

if if

I ðx Þ≥T ; I ðx ÞbT

ð19Þ

where the threshold T is estimated by using the Ostu method [23]. Fig. 4 shows the comparison results of the segmentation accuracy and efficiency of energy minimization. Fig. 4(A) is a synthetic brain

MR image corrupted with intensity inhomogeneity. The first row and the second rows show the curve evolution procedure and bias corrected images of three-phase and four-phase level set formulation respectively. From these results, we can see that the segmentation results and the bias corrected images are similar. The CPU times at the same iterations of these two methods are listed in Table 1. The CPU times are recorded with Matlab code run on a Lenovo computer, with a Pentium 4 processor, 2.80 GHz, 512M RAM, and Matlab R2007 on Windows XP. From the table, we can see that our model is superior in terms of computational efficiency. To quantitatively compare our method with other methods, we use Jaccard Similarity (JS) [24] as an indicator of the segmentation accuracy, which is defined by J ðS 1 ; S 2 Þ ¼

jS 1 ∩S 2 j jS 1 ∪S 2 j

ð20Þ

where | ⋅ | represents the area of a region. S1 is the segmented region by the method, and S2 is the corresponding region given by groundtruth. The closer the JS value to 1, the better segmentation results. We tested our method and the methods of Wells et al., Leemput et al. and Li et al. on 100 simulated images. These images are obtained from McGill Brain Web (http://www.bic.mni.mcgill. ca/brainweb/) with 3% noise level, and 40%, 60%, 80%, 100%

T. Zhan et al. / Magnetic Resonance Imaging 31 (2013) 439–447

A

443

B

C

D

E

Fig. 3. Artificial image. (A) The original image; (B) segmentation result of Li et al.'s model; (C) the bias field of Li et al's model; (D) segmentation result of proposed method; (E) the bias field calculated using Eq. (15).

A

B

C

D

E

F

G

H

I

J

K

Fig. 4. Comparison results between proposed method and four-phase framework on a synthetic image.

444

T. Zhan et al. / Magnetic Resonance Imaging 31 (2013) 439–447

Table 1 Comparison of CPU time for the results in Fig. 5

Four-phase framework Our method

iter = 2

iter = 4

iter = 8

iter= 15

iter = 30

iter = 40

1.159 0.817

2.347 1.626

4.630 3.263

8.844 6.134

17.769 12.223

24.133 16.272

1

0.9

0.9

0.8

0.8

0.7

JS for GM of Axial

JS for GM of Axial

0.7 0.6 0.5 0.4 0.3

0.5 0.4 0.3 0.2

0.2 Wells et al Leamput et al Li et al our method

0.1 0

0.6

40% IUN

60% IUN

80% IUN

Wells et al Leamput et al Li et al our method

0.1 0

100% IUN

40% IUN

60% IUN

1

0.9

0.9

0.8

0.8

0.7

0.7

0.6 0.5 0.4 0.3

0.6 0.5 0.4 0.3 0.2

0.2 Wells et al Leamput et al Li et al our method

0.1 0

40% IUN

60% IUN

80% IUN

Wells et al Leamput et al Li et al our method

0.1 0

100% IUN

40% IUN

60% IUN

C 0.9

0.9

0.8

0.8

100% IUN

0.7

0.7

JS for GM of sagittal

JS for WM of sagittal

80% IUN

D

1

0.6 0.5 0.4 0.3

0.6 0.5 0.4 0.3 0.2

0.2 Wells et al Leamput et al Li et al our method

0.1 0

100% IUN

B

1

JS for GM of conoral

JS for WM of conoral

A

80% IUN

40% IUN

60% IUN

80% IUN

E

100% IUN

Wells et al Leamput et al Li et al our method

0.1 0

40% IUN

60% IUN

80% IUN

100% IUN

F

Fig. 5. Comparison of Jaccard similarity between our method and the other methods based on McGill Brain Web. The axis x denotes the INU values, INU = 40%,⋯, 100%. The axis y denotes the values of JS for WM (left column) and GM (right column).

T. Zhan et al. / Magnetic Resonance Imaging 31 (2013) 439–447

445

A

B

C

D

E

F

G

H

I

J

K

L

M

N

O

P

Q

R

Fig. 6. Comparison of segmentation results for simulated images from McGill Brain Web. Red, green and blue regions represent WM, GM, and CSF respectively. Column 1: original images; Column 2: ground truth segmentations; Column 3–Column 6: the results of Wells et al., Leemput et al., Li et al., and the proposed method respectively.

A

B

C

D

E

F

G

H

I

J

K

L

M

N

O

Fig. 7. Comparison of the results of bias correction for simulated images from McGill Brain Web Column 1: original images; Column 2–Column 5: the corrected images of Wells et al., Leemput et al., Li et al. and the proposed method respectively.

446

T. Zhan et al. / Magnetic Resonance Imaging 31 (2013) 439–447

6. Conclusion

Table 2 Coefficient of variation Experiments

Tissue

Original

Axial

WM GM WM GM WM GM

12.63% 15.29% 17.92% 24.27% 17.34% 26.16%

Coronal Sagittal

Bias correction Wells

Leamput

Li

Our method

6.35% 10.39% 6.70% 12.84% 6.16% 13.82%

6.44% 10.62% 7.01% 11.61% 6.49% 12.92%

5.99% 10.24% 6.81% 11.83% 5.91% 11.73%

5.85% 10.20% 6.66% 11.36% 5.68% 11.39%

intensity inhomogeneity respectively. From Fig. 5, we can see that JS values of our method are better than those of the other three methods. In Fig. 6, we show the segmentation results of brain simulated image with 3% noise level and 60% intensity nonuniformity (INU) as an example. Fig. 7 shows the bias field corrected images estimated by our method and the other three methods. The images in the first row are obtained from McGill Brain Web with 3% noise level and 100% intensity inhomogeneity. The bias corrected images of the method by Wells et al., Leemput et al., Li et al., and the proposed method are shown in columns 2, 3, 4, and 5, respectively. To evaluate the performance of the algorithms for bias correction and segmentation, we use the Coefficient of Variation (CV) as a metric [25]. Coefficient of variance is defined as a quotient between the standard deviation and the mean value of selected tissue class. A good algorithm for bias correction and segmentation should give low CV values for the bias corrected intensities within each segmented region. From Table 2, we can see that the CV values of our method are lower than those of the other methods. It indicates that the bias corrected images of our method are more uniform than those of the other methods. Fig. 8 shows the results of our method on true 3D brain data. The data size is 181 × 217 × 181. The first and second rows show the surfaces of WM and GM, respectively. It can be observed that the good result can be obtained by our method.

In this paper, we use a Gaussian distribution with the information of bias field as the local region descriptor, and propose an improved method for simultaneous MR image segmentation and bias field correction. Benefiting from the local variance, our method can obtain more accurate segmentation results. Moreover, for brain MR image segmentation, we use three-phase level set formulation in our energy function to improve the efficiency of energy minimization. Comparison results and analysis demonstrate the good performance of our method.

Acknowledgments This work is supported by National Nature Science Foundation of China under Grant No. 61071146, No. 61003209 and No. 61171165, China Postdoctoral Science Fund No. 2012M511281, Natural Science Foundation of Jiangsu under Grant No. SBK201022367 and No. BK2012800, Graduate Student Research and Innovation Program of Jiangsu province No. CXZZ1102pt]0.2cm0.5pt259 and Qing Lan Project of Jiangsu Province. The authors would like to thank Dr. Li Wang for suggestions on the manuscript, and would also like to thank the anonymous reviewers for their comments which are valuable in improving this paper.

References [1] Larie S, Abukmeil S. Brain abnormality in schizophrenia: a systematic and quantitative review of volumetric magnetic resonance imaging studies. J Psych 1998;172:110-20. [2] Akbari H, Josygu Y, Kojima K, Tanaka N. Detection and analysis of the intestinal ischemia using visible and invisible hyperspectral imaging. IEEE Tran Biomed Eng 2010;57(8):2011-7. [3] Zijdenbos A, Dawant B. Brain segmentation and white matter lesion detection in MR images. Crit Rev Biomed Eng 1994;22:401-65. [4] Akabari H, Kosugi Y, Kojima K. Segmentation of arteries in minimally invasive surgery using change detection. IEICE Trans Inf Syst 2009;E92-D(3):498-505.

A

B

C

D

Fig. 8. 3D segmentation results of our method for simulated images from McGill Brain Web. The first and second columns show 3D results of WM and GM respectively.

T. Zhan et al. / Magnetic Resonance Imaging 31 (2013) 439–447 [5] Khoo V, Dearnaley D, Finnigan D, Padhani A, Tanner S, Leach M. Magnetic resonance imaging (MRI): considerations and applications in radiotherapy treatment planning. Radiother Oncol 1997;42:1–15. [6] Chan T, Vese L. Active contours without edges. IEEE Trans Imag 2001;10:266-77. [7] Vese L, Chan T. A multiphase level set framework for image segmentation using the Mumford and Shah model. Int J Comp Vis 2002;50:271-93. [8] Rousson M, Deriche R. A variational framework for active and adaptative segmentation of vector valued images. Workshop on Motion and Video, Computing. IEEE; 2002. p. 52-62. [9] Lie J, Lysaker M, Tai X. A variant of the level set method and applications to image segmentation. Mathematics for Computation 2006;75(255):1155-74. [10] Tai X, Li H. A piecewise constant level set methods for elliptic inverse problems. Applied Numerical Mathematics 2007;57(5–7):686-96. [11] Lie J, Lysaker M, Tai XC. A binary level set model and some applications for Mumford–Shah image segmentation. IEEE Trans Image Process 2006;15(5): 1171-81. [12] Wells W, Grimson E, Kikinis R, Jolesz F. Adaptive segmentation of MRI data. IEEE Trans Med Imag 1996;15(4):429-42. [13] Leemput V, Maes K, Vandermeulen D, Suetens P. Automated model-based bias field correction of MR images of the brain. IEEE Trans Med Imag 1999;18(10):885-96. [14] Pham D, Prince J. Adaptive fuzzy segmentation of magnetic resonance images. IEEE Trans Med Imag 1999;18(9):737-52. [15] Li C, Huang R, Ding Z, Gatenby C, Metaxas D, Gore J. A level set method for image segmentation in the presence of intensity inhomogeneities with application to MRI. IEEE Trans Image Process 2011;20(7):2007-16.

447

[16] Li C, Kao C, Gore JC, Ding Z. Minimization of region-scalable fitting energy for image segmentation. IEEE Trans Imag Proc 2008;17(10):1940–19498. [17] Wang L, He L, Mishra A, Li C. Active contours driven by local gaussian distribution fitting energy. Signal Processing 2009;89(12):2435-47. [18] Li C, Xu C, Gui C, Fox M. Distance regularized level set evolution and its application to image segmentation. IEEE Trans Image Process 2010;19(12):3243-54. [19] Paragios N, Deriche R. Geodesic active regions: a new framework to deal with frame partition problems in computer vision. J Vis Commun Image Representation 2002;13:249-68. [20] Wang L, Li C, Sun Q, Xia D, Kao C. Active contours driven by local and global intensity fitting energy with application to brain MR image segmentation. Comput Med Imaging Graph 2009;33(7):520-31. [21] Chen Y, Zhang J, Macione J. An improved level set method for brain MR images segmentation and bias correction. Comput Med Imaging Graph 2009;33(7): 510-9. [22] Jeon M, Alexander M, Pedrycz W, Pizzi N. Unsupervised hierarchical image segmentation with level set and additive operator splitting. Pattern Recognit Lett 2005;26(10):1461-9. [23] Ostu N. A threshold selection filter from grey level histogram. IEEE Trans System Man Cybernet 1972;9:62-6. [24] Vovk U, Pernus F, Likar B. A review of methods for correction of intensity inhomogeneity in MRI. IEEE Trans Med Imag 2007;26(3):405-21. [25] Dawant B, Zijdenbos A, Margolin R. Correction of intensity variations in MR images for computer-aided tissue classification. IEEE Trans Med Imaging 1993;12(4):770-81.