segmentation of three - dimensional anatomical brain images into 
  tissue classes has 
 applications in both clinical and research settings .
although 
 numerous methods to segment brain mri for extraction of cortical 
 white matter , gray matter , and cerebrospinal fluid ( csf ) have been 
 proposed for the past two decades , little work has been done to 
 evaluate and compare the performance of different segmentation 
 methods on real clinical data sets , especially for csf . segmenting three - dimensional anatomical brain images into tissue 
 classes
, segmentation can provide volumetric 
 quantification of cortical atrophy and thus aid in the diagnosis 
 of degenerative diseases .
these volumetric measurements apply to 
 the research setting as well , where segmentation can also be used 
 to define regions of interest for quantifying the physiological 
 responses measured with fmri or pet acquired on the same patients 
 and coregistered with the mri data .
clinical studies based on 
 quantitative measurements of cortical brain structures include 
 alzheimer 's disease [ 13 ] , epilepsy , schizophrenia 
 , cerebrovascular deficiency , and multiple 
 sclerosis .
a good review of these methods can be found in 
 and we can distinguish two general families .
( 1 ) statistical methods the first family of methods is 
 based on classification of brain tissues into different 
 classes , based on intensity values ( direct values of features 
 computed from these values ) .
one common difficulty 
 with this method is the selection of the threshold level . 

many selection approaches have been proposed based on histograms 
 [ 10 , 11 ] , combination of morphological operators and region 
 growing , and so forth .
derived from a statistical 
 framework , bayesian analysis [ 13 , 14 ]   is a popular 
 classification method for brain tissue where automatic 
 segmentation is performed with expectation maximization ( em ) 
 .
the freely distributed statistical parametric mapping 
 ( spm ) software tool   is widely used by neuroradiologists 
 for clinical research .
additional classification methods include 
 clustering [ 16 , 17 ] , fuzzy classification , neural 
 networks , deterministic annealing .
hybrid 
  neuro - fuzzy  methods , see [ 2 , 20 , 21 ] , combining fuzzy 
 logic and neural networks perform an unsupervised learning 
 process . another class of statistical methods is based on markov 
 random field ( mrf ) [ 2225 ] .
mrf model encodes spatial 
 information through the mutual influences of neighboring voxels 
 for class assignments .
a major issue with mrf - based classification 
 methods is the requirements for training the model and setting the 
 mrf parameters which typically require supervised learning and 
 a priori information from manual labeling or from an 
 atlas . in this context , classification and nonuniform registration 

finally , these statistical methods can be applied to multispectral 
 mri ( i.e. , mri data of the same patient acquired with different 
 protocols such as t1-weighted , t2-weighted , proton diffusion ) 
 using multivariate statistics .
several applications of 
 multispectral brain mri classification have been presented in the 
 literature [ 2527 ] .
( 2 ) deformable models the second family of segmentation methods deals with 
 geometric deformable models , including active surfaces 
  and level - set - based deformable models .
the level set 
 implementation framework for surface propagation offers the 
 advantages of easy initialization , computational efficiency , and 
 the ability to capture deep sulcal folds .
two coupled level set 
 surfaces were proposed by zeng et al .   for cortex 
 segmentation from 3-d mr images , assuming a constant thickness 
 range of the cortical mantle . a combination of joint - prior shape 
 appearance models with a level set deformable model was proposed 
 by yang and duncan .
this method was motivated by the observation that the 
 shapes and gray levels variations in an image had some consistent 
 relations building a map shape - appearance prior model provided 
 some configurations and context information to assist the 
 segmentation process .
the model was formulated in a level set 
 framework rather than using landmark points for parametric shape 
 description .
goldenberg   et al .   proposed a 
 geometric variational formulation for the propagation of two 
 coupled bounding surfaces , similar to zeng et al . . 
 
the authors put forward an 
 efficient numerical scheme for the implementation of a geodesic 
 active surface model .
several external driving forces , derived 
 from brain mri , have been proposed based on image gradient 
 , image intensity , and probability density 
 function   of tissue classes . combining
combining bayesian analysis and active surface method , and shen 
  et al .   combining geometric deformable model and 
 statistical information about the shapes of interest .
automation of the segmentation process is critical for 
 applications in clinical research where the number of cases to 
 process is large and the time available for experts to analyze the 
 data is very limited .
several of the aforementioned methods were 
 developed in a supervised or semi - automated framework still 
 requiring operators ' intervention .
full automation can be achieved 
 using automatic parameter tuning , automated 
 initialization , and combination with atlas 
 information .
segmentation methods : histogram - based thresholding , tissue 
 classification based on fuzzy connectedness , and 
 maximum - likelihood classification with hidden mrf using the 
 fsl - fast software tool .
we also present in this paper the 
 implementation of a new three - dimensional automated segmentation 
 method of brain mri using a four - phase three - dimensional active 
 contour implemented with a level set framework .
this multiphase 
 level set framework was initially proposed by vese and 
 chan   to simultaneously deform coupled level set 
 functions without any prior models or shape constraints . 
 
this framework achieves a global 
 partitioning of the image data into 2 homogeneous areas 
 using n level set curves , solely based on average gray values 
 measures .
these four segmentation methods were applied to a set of ten 
 clinical t1-weighted mris for segmentation of cortical tissues : 
 white matter ( wm ) , gray matter ( gm ) , and cerebrospinal fluid 
 ( csf ) .
the segmentation methods were also compared in a 
 statistical framework to assess their relative performance and 
 overall  quality . 
we present in this section the four segmentation methods that were 
 applied to ten brains t1-weighted mri and compared together . 

this  ground truth  was used to optimize the parameter settings 
 of the three  classical  segmentation methods , as detailed below .
intensity thresholding ( it ) is the easiest and fastest 
 segmentation method , often adopted for preprocessing of medical 
 images and preregistration problems .
segmentation of the 
 three brain cortical tissues is performed via thresholding of 
 voxel values within adjacent intervals .
the position of the 
 interval bounds was initialized as follows : we used the manually 
 labeled data to mask the mri data and compute the means of the 
 three cortical tissues of interest .
these mean values were used to 
 initialize the threshold values at the two interfaces csf / gm and 
 gm / wm .
finally , the simplex method 
  was applied to maximize the tanimoto value and identify 
 optimal interval bounds corresponding to a minimization of the 
 segmentation error .
classification of homogeneous objects using fuzzy connectedness 
 ( fc ) was introduced by falco et al . . 
 
this method can be used to 
 extract homogeneous tissues in a volume given an initial set of 
 seed points and prior statistics defining these tissues .
the 
 source code for this method was obtained from the 
 freely distributed insight segmentation and registration toolkit 
 ( itk )   sponsored by the national library of medicine . 

segmentation of the different tissues was performed via 
 thresholding of the real - valued fuzzy connectedness maps .
the 
 threshold levels were also determined via minimization of the 
 segmentation error using the simplex method .
the fast software tool   is a publicly available automated 
 segmentation tool for brain mri volumes .
this software tool is 
 part of the fsl comprehensive library of functional and structural 
 brain image analysis tools , distributed by the oxford centre for 
 functional magnetic resonance imaging of the brain ( fmrib ) at the 
 university of oxford , uk .
brain mri volumes were segmented into 
 three tissue types ( wm , gm , and csf ) with simultaneous correction 
 of rf bias - field .
the segmentation process is based on hidden 
 markov random field ( hmrf ) models .
fitting of the hmrf model to 
 the data was performed via maximum likelihood with expectation 
 maximization ( em ) .
the multiphase three - dimensional level set 
 ( m3dls ) segmentation method performing a minimal 
 partitioning of the image data into piecewise constant objects , 
 based on the mumford - shah functional , was introduced by 
  chan and vese .
this method uses a deformable model controlled by a 
 homogeneity - based energy functional to segment piecewise constant 
 or piecewise smooth volumetric data u0 . assuming a 
 piecewise - constant data with an object , of value c1 , and a 
 background , of value c2 , separated by the contours c , the 
 proposed energy functional is defined as 
 ( 1)f(c , c1,c2)=(length(c))+(area(inside  c ) )  + 1inside  ( c)|u0c1|2d  + 2outside  ( c)|u0c2|2d , where   0 ,   0 , 1 , 2 > 0 are fixed parameters .
segmentation of the data is performed via minimization of the 
 functional f with respect to ( c , c1 , c2 ) .
this energy 
 functional can be extended to the segmentation of multiple 
 homogeneous objects in the image by using several curves 
 { c1 , c2 ,  , ci}. in the case of two curves the following 
 energy functional is used : 
 ( 2)f(c1,c2,c00,c01,c10,c11 )  = 1  length  ( c1)+2  length  ( c2 )  + 1  area  ( inside  c1)+2area  ( inside  c2 )  + 1inside  c1 ,  inside  c2|u0c11|d  + 2inside  c1 ,  outside  c2|u0c10|d  + 3outside  c1 ,  inside  c2|u0c01|d  + 4outside  c1 ,  outside  c2|u0c00|d. 
 the set of parameters ( 1 , 2 , 3 , 4 , 1 ,   1 , 2 , 2 ) 
 takes real positive values .
the two closed curves ( c1 , c2 ) split the domain  into four phases defined by their 
 relative positions as   illustrated in 
 figure 1 . inside these four phases , u0 has mean 
 intensity values ( c00 , c01 , c10 , c11 ) .
minimization of this energy functional deforms simultaneously the 
 two curves and identifies four homogeneous areas defined by the 
 intersection of the two curves .
minimization of the functional in ( 1 ) and ( 2 ) can be performed within a level set framework .
this framework , 
 introduced by osher and sethian , provides an effective 
 implicit representation for evolving curves and surfaces , which 
 has found many applications in image segmentation , denoising and 
 restoration as reviewed in . in this framework ,
a given 
 curve c ( being now the boundary w of an open set 
    ) is represented implicitly , as the zero level 
 set of a scalar lipschitz function 
   :    ( called level set function ) , such that 
 ( 3)(x)>0  on  ,(x)<0  on  \,(x)=0  on  . 
 the level set function  is typically defined as the signed 
 distance function of spatial points defined on  to the 
 curve c. once  is computed , we define its associated 
 heaviside function h and dirac function    as 
 ( 4)h()={0 ,  if  0,1 ,  if  <0,()=ddh( ) . 
 using these two functions , the different components of the 
 functional in ( 1 ) , parameterized with the contour curve 
 c , can be reformulated as integrals , parameterized with the 
 level set function  and defined over the entire domain 
 . 
 length of the curve c : 
 ( 5)length(c)=length(=0 )  = |h()|d=()||d.area inside the curve c : 
 ( 6)area(c)=area(<0)=h( )  d.homogeneity of u0 inside and outside the curve c : 
 ( 7)inside(c)|u0c1|2d=|u0c1|2h()d,outside(c)|u0c2|2d=|u0c2|2(1h())d. 
 minimization of the energy functional leads to the expression of 
 the associated euler - lagrange equations for derivatives with 
 respect to ( c1 , c2 ,  ) .
the first two derivatives , with 
 respect to ( c1 , c2 ) , lead to the following results , used 
 for computation of the mean statistics : 
 ( 8)c1()=u0(x)h((x))dh((x))d,c2()=u0(x)(1h((x)))d(1h((x)))d. 
 the third partial derivative of   f( ,
c1 , c2 ) , with respect 
 to  , leads to the following dynamic equation : 
 ( 9)t=()[  div(||)+1|u0c1|22|u0c2|2 ] . 

this segmentation framework is extended to the detection of 
 multiple objects via the introduction of multiple level set 
 functions { 1 , 2 ,  } and the computation of mean 
 data values in areas of constant values defined via the 
 combination of their heaviside functions ( h(1 )  h(2 )   ) . in this study
we implemented the 
 segmentation functional with two level set functions generating 
 four phases .
we can introduce the binary functions  ( also 
 called characteristic functions ) that define the four regions : 
 ( 10)1()=h( ) ,  0()=(1h()),11(1 ,  2)=1(1)1(2 ) ,  10(1 ,  2)=1(1)0(2),01(1 ,  2)=0(1)1(2 ) ,  00(1 ,  2)=0(1)0(2 ) . 

the euler - lagrange equations , for the four - phase configuration , 
 lead to the following equations defining the mean values of each 
 phase : 
 ( 11)c11=u011(1 ,  2)d11(1 ,  2)d ,  c10=u010(1 ,  2)d10(1 ,  2)d,c01=u001(1 ,  2)d01(1 ,  2)d ,  c00=u000(1 ,  2)d00(1 ,  2)d. 
 the partial derivatives with respect to the two level set functions 
 define the following system of equations : 
 ( 12)1t=(1){  div  ( 1|1|)  + ( 1(u0c11)23(u0c01)2)1(2 )  + ( 2(u0c10)24(u0c00)2)0(2)},2t=(2){  div  ( 2|2|)  + ( 1(u0c11)22(u0c10)2)1(1 )  + ( 3(u0c01)24(u0c00)2)0(1)}. 
 we present in the next section details regarding the numerical 
 implementation of this system .
length of the curve c : 
 ( 5)length(c)=length(=0 )  = |h()|d=()||d. area inside the curve c : 
 ( 6)area(c)=area(<0)=h( )  d. homogeneity of u0 inside and outside the curve c : 
 ( 7)inside(c)|u0c1|2d=|u0c1|2h()d,outside(c)|u0c2|2d=|u0c2|2(1h())d. segmentation was performed via iterations of the system of 
 equations defined in ( 12 ) , deforming the two level set 
 fronts ( 1 , 2 ) until convergence to a stable position 
 was reached .
initialize the system for time n = 0 , with ( 1 , 2 ) defined as the distance functions from an initial set of curves.for time n > 0 , while the system is unstable , compute the average values in the four phases 
 ( c11 , c10 , c01 , c00);compute the curvature and homogeneity terms of the speed , defined 
 for each points in the spatial domain ;compute ( 1 , 2 ) from ( 12 ) with the speed term ( from ( b ) ) 
 and ( 1 , 2);evaluate the stability of the system . 

iterate at time n + 1 if the system is not stable.when the system is stable , extract the four phases as binary 
 volumes corresponding to 
 ( 11 , 10 , 01 , 00 ) .
initialize the system for time n = 0 , with ( 1 , 2 ) defined as the distance functions from an initial set of curves . for time
n > 0 , while the system is unstable , compute the average values in the four phases 
 ( c11 , c10 , c01 , c00 ) ; compute the curvature and homogeneity terms of the speed , defined 
 for each points in the spatial domain  ; compute ( 1 , 2 ) from ( 12 ) with the speed term ( from ( b ) ) 
 and ( 1 , 2 ) ; evaluate the stability of the system . 

when the system is stable , extract the four phases as binary 
 volumes corresponding to 
 ( 11 , 10 , 01 , 00 ) . in our implementation ,
the segmentation was initialized with two 
 level set functions defined as the distance function from two sets 
 of initial curves .
the curves were defined as 64 cylinders 
 centered at regularly spaced seed locations across the entire data 
 volume .
the two sets of cylinders were slightly shifted from each 
 other as illustrated in figure 2 .
note that such initialization does not use any a priori information on the location of particular tissues or on the 
 anatomy of the brain and does not require manual input by the 
 user . as commented by vese and chan
, this type of 
 initialization also brings some robustness to the method , limiting 
 the risks of convergence of the minimization process to local 
 minima .
the level set algorithm was implemented with a regularized 
 heaviside function and a semi - implicit scheme as proposed by vese 
 and chan   and extended to three dimensions . in this 
 scheme ,
the speed term at time n was computed with the curvature 
 derived semi - implicitly ( using ( i , i)i = 1 , 2 ) 
 while the homogeneity force term was computed explicitly with 
 ( ci)i = 1 , 2 .
this semi - implicit scheme provides 
 unconditional stability for any temporal and spatial 
 discretization parameters .
this means that we can set the time 
 increment to an arbitrary high value to speed up the segmentation 
 process without altering the stability of the minimization 
 process .
all four segmentation methods tested in this work perform a 
 partitioning of the volumetric data into three tissue classes and 
 a background relying on a strong assumption of tissue homogeneity 
 for wm , gm , and csf .
the notion of homogeneity is then translated 
 into a statistical framework ( hmrf - em , homogeneous tissue is 
 characterized by its mean and variance ) , a fuzzy connectivity ( fc , 
 homogeneous tissue is characterized by high affinity measures to a 
 prior seed point given a statistical model ) , or a distance measure 
 ( m3dls , homogeneous tissue is characterized by its mean intensity 
 value ; it , interval bounding values ) .
even though homogeneity - based segmentation methods are widely used 
 for brain mri segmentation , it is well known that there are strong 
 tissue inhomogeneities in mri volumes of the following four 
 origins : ( 1 ) biological tissues are inherently heterogeneous with 
 internal structures and multimolecular components which 
 contribution are integrated in the recorded mri signal , ( 2 ) the mri 
 acquisition system is degraded by inherent statistical noise , ( 3 ) 
 design of the mri acquisition system suffers from inhomogeneity of 
 the radio - frequency ( rf ) field ( this phenomenon is also referred 
 to as rf bias - field ) , ( 4 ) mri imaging systems have 
 limited spatial resolutions which generate partial volume 
 effects ( i.e. , mixing of mri signals from adjacent tissues ) at 
 tissue interface locations .
the first two sources are inherent to the modality but have minor 
 impacts on the cortical structures that we are trying to segment . 

the other two sources generate artifacts that tend to become 
 undetectable by simple visual inspection as mri scanner technology 
 improves , but that still constitute the major source of error and 
 failure of homogeneity - based segmentation algorithms . 

rf bias - field : elaborated algorithm has been developed to 
 correct rf bias - field [ 4649 ] .
the first 
 method was proposed by styner et al .   and is available in 
 the itk package .
the second method was proposed by 
 ashburner and friston   and is available in the 
 fsl software package .partial volume effect introduces inhomogeneities at 
 tissue interfaces that can be modeled in a statistical framework 
 by manipulating tissue mixture models as in [ 24 , 51 , 52 ] . 

mixed classes are created and if necessary assigned to one of the 
  true 
tissue classes constituting the mix . neither the 
 statistical hmrf - em method nor the three other methods tested in 
 this work included such correction and initial experiments on 
 clinical data sets revealed misclassification of tissue labels for 
 voxels located at tissue interfaces .
we therefore derived a 
 postprocessing sequence , based on morphological operators to 
 correct for misclassifications at interfaces between gm and wm and 
 between the csf and wm or gm .
rf bias - field : elaborated algorithm has been developed to 
 correct rf bias - field [ 4649 ] .
the second method was proposed by 
 ashburner and friston   and is available in the 
 fsl software package .
partial volume effect introduces inhomogeneities at 
 tissue interfaces that can be modeled in a statistical framework 
 by manipulating tissue mixture models as in [ 24 , 51 , 52 ] . 

mixed classes are created and if necessary assigned to one of the 
  true 
tissue classes constituting the mix . neither the 
 statistical hmrf - em method nor the three other methods tested in 
 this work included such correction and initial experiments on 
 clinical data sets revealed misclassification of tissue labels for 
 voxels located at tissue interfaces .
we therefore derived a 
 postprocessing sequence , based on morphological operators to 
 correct for misclassifications at interfaces between gm and wm and 
 between the csf and wm or gm .
a simple postprocessing scheme was designed to correct for pixel 
 assignment at tissue interfaces . after the level set segmentation 

was completed , wm , gm , and csf structures , corresponding to 
 separate phases , were saved as binary volumes .
these volumes were 
 then used as masks applied to the original data and a gaussian fit 
 of the histograms of each phase was performed .
each phase was then 
 characterized by the mean and standard deviation 
 ( phase , phase ) of the fitted 
 gaussian distribution .
first , the gm mask was dilated , to 
 correct for the under segmentation of this thin structure .
all new voxels , from the dilated mask , with 
 intensity values inside this interval were included in the gm 
 phase while voxels with intensity values outside the interval were 
 removed from the phase and assigned to the adjacent wm phase .
this 
 process was iterated until no new points were added to the gm 
 phase .
a similar process was then applied to the csf phase with 
 dilation of the binary mask .
finally , a 3d 
 connectivity algorithm was performed to remove isolated voxels in 
 the three phases .
this simple postprocessing approach has provided 
 very robust performance on the ten clinical mri cases segmented 
 for this study .
a second postprocessing correction was needed to compare results 
 from the different segmentation methods to manual labeling . 

indeed , manual labeling did not include sulcal csf and labeled the 
 outer most layer of cortical brain tissue as gray matter .
this was 
 an arbitrary assignment made by the manual expert , that does not 
 take into account partial volume effects in these voxels .
a 
 similar arbitrary classification had to be applied to the 
 segmented data , before comparison to manual labels .
this problem 
 was especially important for the fast hmrf - em classification which 
 is based on an anatomical model with sulcal csf , leading to 
 overestimation of the csf on the mri data .
this oversegmentation 
 of sulcal csf lead to very high false positive errors when 
 compared to manual labeling .
sulcal csf was removed from the 
 segmented data , after postprocessing of interface voxels , for all 
 segmentation methods , as follows .
given the spatial resolution of 
 the data ( 0.86 mm  3 mm slice thickness ) , a 
 2-voxels dilation of the manually labeled csf ventricles was 
 performed in axial views , and one 1-voxel dilation was performed 
 in the longitudinal direction .
the dilated ventricle masks were 
 applied to the segmented csf mask and segmented voxels outside the 
 dilated manual mask were assigned to the gm phase .
we applied our segmentation to one phantom and ten t1-weighted mri 
 volumes acquired on healthy young volunteers .
axial slices were 
 1.5 mm thick with an in - plane ( hyphen ) resolution of 
 0.86 mm .
these images were resliced coronally ( 3 mm slice 
 thickness ) and labeled via a labor - intensive ( 40 hours per brain ) 
 manual method in which expert raters with extensive training in 
 neuroanatomy choose histogram thresholds on locally hand - drawn 
 regions of interest .
this labeled data was used as a  ground 
 truth  for evaluation of the segmentation accuracy .
mri volumes were preprocessed to remove all noncortical brain 
 tissue by using the manually labeled data sets as binary masks . 

to 
 determine the practicality of masking out subcortical gray matter 
 on nave images , we constructed from a library of 
  labeled atlases two probabilistic atlases in which each 
 voxel was assigned a likelihood of being made of cortical gray 
 matter or subcortical gray matter .
less than 0.1% of the voxels 
 in the whole brain simultaneously had over 20% chance of being 
 made of both gray matter and subcortical gray matter .
such 
 statistical finding confirms that one can apply a 
 population - based method for masking out subcortical gray matter 
 for the purpose of applying cortical segmentation methods without 
 introducing significant errors .
a mathematical brain phantom was built from one manually labeled 
 mri data set to validate the performance of the multiphase level 
 set segmentation algorithm in an ideal case .
a constant intensity 
 value was applied for each tissue corresponding to the average 
 intensity value of the mri data under the manual mask .
prior to segmentation , we evaluated the homogeneity of the three 
 cortical tissues : wm , gm , and csf , computing the mean and variance 
 statistics of the intensity distribution for each tissue within 
 each slice .
these statistical measurements were performed on axial 
 slices across the entire data volumes for the ten mri cases 
 available for the study , masking the data with the manually 
 labeled data .
results , illustrated in figure 4 for three typical 
 cases , showed very stable estimates of intensity mean and variance 
  values for each tissue across the whole 
 volume , corroborating the accuracy of the homogeneity assumption 
 for the three tissues , and the absence of strong bias - field 
  inhomogeneity .
lower mean intensity values on extremity slices 
 were computed on small tissue samples ( less than 10 voxels ) 
 corresponding to small anatomical structures with relatively high 
 partial volume effects .
a gaussian fit was performed on the histogram of the entire gray 
 level distribution of the three brain tissues for each case .
we 
 observed that the three tissue types have well - separated average 
 values suggesting that global homogeneity measurements could 
 separate tissue types for each patient .
nevertheless , the 
 agreement between the volume histograms and the fitted gaussian 
 distribution was evaluated with a chi - squared test .
results for 
 the different tissues did not show a systematic agreement between 
 the data and the gaussian fit , at level 0.05 , except for the gray 
 matter . therefore , despite reasonable homogeneity , we need further 
 investigation before being able to introduce additional 
 constraints based on a priori gaussian statistics to the 
 method as proposed , for example , by baillard et al . .
segmentation errors were measured using the recent methodology 
 from udupa et al .   for comparison of 
 segmentation methods .
accuracy of the object contours obtained 
 with the proposed segmentation method , referred to as cseg , was evaluated by comparing the results to our  ground truth  segmentation of each object , using 
 manually labeled contours , referred to as cml .
the 
 overlap and difference between the two contours were measured via 
 counting the true positive ( tp ) , false positive ( fp ) , 
 and false negative ( fn ) voxels as illustrated in 
 figure 5 .
these quantities are reported as volume 
 fractions ( vf ) of the true delineated object volume cml as follows . 

fnvf = fraction of tissue defined in cml that was missed by the segmentation method cseg.fpvf = fraction of tissue defined in cseg    
 falsely identified by the segmentation method.tpvf = fraction of the total 
 amount of tissue in the segmentation cseg which 
 overlaps with the true object cml .
fnvf = fraction of tissue defined in cml that was missed by the segmentation method cseg .
fpvf = fraction of tissue defined in cseg    
 falsely identified by the segmentation method .
tpvf = fraction of the total 
 amount of tissue in the segmentation cseg which 
 overlaps with the true object cml .
we point out here that tpvf + fnvf = 1 , so that we 
 only need to observe tpvf and   fpvf measures to 
 evaluate the segmentation method .
in order to evaluate the two bias - field correction methods and 
 their effects on the homogeneity of the tissues , we applied the 
 intensity thresholding ( it ) segmentation before and after 
 bias - field correction , providing three groups of results for the 
 ten mri cases : segmentation of original data ( raw ) , segmentation 
 of corrected data with the itk software tool ( itk ) , and 
 segmentation of corrected data with the fsl software tool ( fsl ) . 

the it segmentation method was selected as it only relies on good 
 separation of the tissue range of intensities , without any 
 additional constraint or model .
the fsl bias - field correction did 
 not need any parameter setting , and the corrected data was 
 generated along with the segmentation via hmrf - em .
the bias - field 
 correction with the itk toolkit used the following parameters : 
 mean and standard deviation of the three tissue classes were 
 calculated from the  ground truth ,  the input and output mask 
 with nonbrain tissues removed were also generated from the 
  ground truth ,  the degree of the method was set to 2 , growing 
 parameter was set to 1.05 , shrinking parameter was set to 0.9 
 ( suggested values ) , and a maximum of 2000 iterations were 
  performed with 
 an initial step size of 1.02 .
effects of bias - field correction on it segmentation are 
  illustrated with bar plots of tpvf , in figure 6 .
we observed that the bias - field correction methods did not improve 
 and even degraded the tpvf accuracy of the it segmentation method . 

the itk method provided similar accuracy for all cases and all 
 tissues while the fsl method tended to degrade tp performance , 
 except for the csf on one mri case .
based on homogeneity results 
 presented in figure 4 , these results and the 
 difficulty involved with parameter settings of the bias - field 
 correction methods , we decided to exclude bias - field correction 
 for this study . the level set segmentation was initialized for the phantom and mri 
 clinical cases with two sets of regularly spaced cylinders as 
 illustrated in figure 2 .
segmentation of the brain phantom was performed with three phases 
 corresponding to the cortical tissues and the 4th phase to the 
 background . in table 1 , we observed almost perfect 
 segmentation of the brain structures , demonstrating the inherent 
 ability of the multiphase level set framework , initialized with 
 small regular cylinders , to the following . 

identify the four homogeneous phases.extract highly convoluted surfaces.perform topology splitting and 
 merging of the evolving front , enabling the identification of 
 separate structures within a single phase 
 ( such as the three ventricles for the csf ) that exist as different spatial objects .
perform topology splitting and 
 merging of the evolving front , enabling the identification of 
 separate structures within a single phase 
 ( such as the three ventricles for the csf ) that exist as different spatial objects .
computation of homogeneity - based speed terms for phases with small 
 structures , such as the csf , suffered from higher estimation 
 inaccuracy . in the absence of a curvature term , the 
 homogeneity - based speed terms generate a k - means classification 
 of the image data with poor quality extraction of these small 
 structures . in this context , the curvature term of the deformable 
 model plays a critical role in the m3dls to preserve small 
 structures shapes and sizes .
partial volume effect was not 
 simulated in the phantom and therefore no postprocessing was 
 applied .
we used this phantom to tune the parameters of the segmentation 
 method set to 
 ( 13)1=2=3=4=0.01 ,  =0,=4.108  volume_size / diagonal_distance,t=104 ,  x=y=z=1 . 

the parameter associated with the curvature term is defined 
 proportional to the data volume size ( volume_size ) and 
 inversely proportional to the diagonal distance of the volume data 
 ( diagonal_distance ) . by doing so , we consider this 
  diagonal distance as the unitary distance of our domain of 
 definition . setting the constant speed term  to zero eliminates the use of a constant inflating force on the 
 model .
this type of constant force should be used with caution as 
 it can override the homogeneity - based speed term when driving the 
 deformable contour and force it to move in only one direction 
 ( i.e. , constant inflation or deflation ) .
we present in figure 7 a three - dimensional rendering 
 of each segmented structure for one mri clinical case segmented 
 for this study .
visual rendering of the three cortical structures confirmed the 
 overall high performance of the multiphase segmentation method to 
 extract homogenous objects that correspond to distinct anatomical 
 tissues .
the segmentation method was able to handle multiple challenges 
 without any a priori information or shape constraints 
 that include the extraction of highly convoluted white matter 
 surfaces , the extraction of separate ventricular 
 structures for the csf , and handling of different volume sizes of 
 the three structures in a simultaneous segmentation scheme .
accuracy of the object contours obtained with the proposed 
 multiphase level set deformable model and the three other 
 segmentation methods was evaluated by comparing segmentation 
 results to the manual ground truth .
bar plots of error 
 measurements for the segmentation of the ten clinical cases 
 and average errors are reported in figures 8 , 
 9 , and 10 . because tpvf and fnvf are both 
 relative measures with respect to the ground truth volume , we have 
  tpvf + fnvf = 1 .
we can make the following observations based on the results 
 presented in these graphs for segmentation of gm , 
 characterized as a thin structure with large surface area : hmrf 
 generates the highest tpvf and fc the lowest .
it and m3dsl perform 
 well ( tpvf > 85% ) and stably over the 10 cases .
fc 
 generates over - segmentation of the gm with high fpvf ( up to 
 48% ) , followed by hmrf - em ( up to 25% ) .
m3dsl ( fpvf < 
 18% ) generates low and stable over - segmentation errors , and it 
  ( fpvf < 10% ) does not generate this type 
 of error .
overall , it provides the best performance for 
 tpvf and fpvf , followed by m3dls , with good performance .
we can make the following observations based on the results 
 presented in these graphs for segmentation of wm ( characterized as a large structure ) : it and m3dsl generated high 
 stable and similar tpvf ( tpvf > 90% ) , fc generated low 
 tpvf on 2 cases , while hmrf - em systematically generated lower 
 tpvf .
we can make the following observations based on the results 
 presented in these graphs for segmentation of the csf 
 ( characterized as small disconnected structures ) : hmrf - em was the 
 only segmentation method to show very high tpvf ( > 90% ) .
all 
 other three methods generated variable errors over the 10 cases 
 with m3dls generating the smaller range of errors 
 ( 50% < tpvf < 80% ) .
in particular , it performed very poorly , 
 certainly due to the absence of shape constraint .
hmrf showed 
 large over - segmentation errors ( above 10% and up to 30% ) , 
 while m3dsl showed very low fpvf errors ( < 5% ) . 
 
overall , hmrf - em showed perfect sensitivity along with 
 poor selectivity ( i.e. , over segmented the csf ) and m3dls showed 
 the best tradeoff of sensitivity versus specificity .
the four segmentation methods were compared statistically using a 
 characteristic index of their performance .
this index is a 
 quantitative parameter used to evaluate the segmentation results 
 and is defined as 
 ( 14)ti = tpvf1+fpvf . because ti populations do not follow a normal distribution , a 
 nonparametric analysis was performed for the four methods 
 over the 10 segmented cases , measuring the differences of the ti 
 indexes .
small   p values ( below 0.05 ) indicate a significant 
 statistical difference between the methods . 

distributions of the ti index over the 10 cases for each method 
 are plotted in figure 11 .
whiskers are 
 lines extending from each end of the box to show the extent of the 
 rest of the data ( 1.5 interquartile ranges ) .
outliers , identified 
 with red crosses , are data with values beyond the ends of the 
 whiskers .
if there is no data outside the whisker , a dot is placed 
 at the bottom whisker .
these results illustrate graphically the average performance and 
 the variability of performance of individual segmentation methods . 

we observed that m3dls and hmrf - em generated the most stable 
 performance over the three tissue types , with mean 
 performance higher with m3dls for wm and gm .
evaluation of the statistical differences of results from the four 
 segmentation methods was performed with a wilcoxon signed rank 
 test for paired data .
significance values for the three tissue 
 types for each pair of segmentation methods are reported in 
 table 2 . from these results
we see the following . 
 for gm , the two best methods , m3dsl and it were not statistically different.for wm , the two best methods , m3dsl and it were statistically different .
the 
 two weakest methods , fc and hmrf - em were not statistically different.for csf , the two best methods , m3dsl and hmrf - em were statistically different . 

therefore , for 
 csf , hmrf - em is superior to the three other methods . compared with it and fc
, m3dls 
 has a higher mean , and smaller variance , which shows better performance and more 
 stability . for gm , the two best methods , m3dsl and it were not statistically different . for wm , the two best methods , m3dsl and it were statistically different .
the 
 two weakest methods , fc and hmrf - em were not statistically different .
for csf , the two best methods , m3dsl and hmrf - em were statistically different
. 
 on the other hand , fc was not statistically different from m3dsl and it .
therefore , for 
 csf , hmrf - em is superior to the three other methods .
compared with it and fc , m3dls 
 has a higher mean , and smaller variance , which shows better performance and more 
 stability .
specific problems for segmentation of the csf phase as observed in the error plots and illustrated in 
 figure 11 , all methods except hmrf - em performed 
 significantly poorly on csf than on wm and gm , corresponding to 
 under segmentation of the ventricles , whose pixels were assigned 
 to white matter . on the other hand ,
very low resolution 
 at the ventricle borders explains in part this result . in 
 addition ,
manual labeling of the mri data for the ventricle can 
 also bear some error as localizations of its borders are difficult 
 even for an expert performing manual tracing .
we illustrated in 
 figure 12 csf segmentation with the different methods 
 with a three - dimensional rendering of the two lateral ventricles 
 and a section of the third ventricle .
we also recall that manual labeling is only used as a method of 
 reference , and does not provide a real  ground truth  to the 
 segmentation problem . in that context , kikinis et al . reported in   a variation in volumetric measurements 
 between manual observers in the order of 15% for wm , gm , and 
 csf .
this study used a 3d implementation of a segmentation framework 
 initially proposed by chan and vese .
a single 
 illustration of the ability of the method to segment a brain mri 
 slice was illustrated in their paper .
no group has previously 
 implemented and tested this method on a whole - brain mri dataset , 
 with fixed parameter settings , which is a critical aspect for the 
 demonstration of the ability of such a segmentation approach to be 
 an efficient clinical tool .
we chose to use a whole - brain mri data 
 set in order to establish a benchmark of the performance of 
 existing brain mri segmentation tools , as well as to develop a new 
 automated and robust segmentation tool that would alleviate the 
 need for a priori knowledge such as tissue statistics and 
 the need for manual initialization .
we can compare our segmentation error to results reported by zeng 
  et al .   and niessen et al . .
in zeng 
  et al . , the authors tested their algorithm for the 
 segmentation of frontal lobes on seven high - resolution mri 
 datasets from a randomly chosen subset of young autistic and 
 control adult subjects .
they ran a coupled - surfaces level set 
 algorithm to isolate the brain tissue and segment the cortex .
the 
 average tpvf and fpvf for the cortical gm in the frontal lobe were 
 86.7% and 20.8% .
in niessen et al . , a 
  hyperstack  segmentation method , based on multiscale pixel 
 classification , was tested for 3d brain mri segmentation .
a 
 supervised segmentation framework with manual postediting was 
 applied to a probabilistic brain phantom for estimation of 
 segmentation error .
first , a binary segmentation of the brain 
 phantom was performed to evaluate the minimal segmentation error 
 due to partial volume effects .
the study reported a volume 
 fraction of misclassified pixels ( fpvf + fnvf ) 
 around 20% for wm , gm , and csf .
optimal 
 ( fpvf + fnvf ) errors were obtained with the 
 probabilistic version reporting 10% for wm , 21% for gm , and 
 25% for csf . in our case
we proposed a fully automated segmentation method with 
 no a priori information , leading to the following 
 conclusions . 
 for gm
, we showed that the m3dls method was statistically 
 equivalent to the idealized it with average ( fpvf + fnvf ) 
 equal to 11.39% for it and 14.23% for m3dls.for wm , m3dls produced an 
 average ( fpvf + fnvf ) equal to 13.56% , superior 
 to hmrf - em with 21.64% and fc methods with 
 19.57%.for csf , m3dls , which produced an average ( fpvf + fnvf ) 
 equal to 32.56% , was found statistically superior to fc with 
 ( fpvf + fnvf ) equal to 45.0% and 48.92% for 
 it . 

the random initialization ensured robustness of the method to 
 variation of user expertise and biased or erroneous input 
 information that could be influenced by variation in image quality 
 or user expertise .
the method is automated , given that the set of 
 parameter values selected on the brain phantom is suited for the 
 clinical data . for gm
, we showed that the m3dls method was statistically 
 equivalent to the idealized it with average ( fpvf + fnvf ) 
 equal to 11.39% for it and 14.23% for m3dls . for wm
, m3dls produced an 
 average ( fpvf + fnvf ) equal to 13.56% , superior 
 to hmrf - em with 21.64% and fc methods with 
 19.57% .
for csf , m3dls , which produced an average ( fpvf + fnvf ) 
 equal to 32.56% , was found statistically superior to fc with 
 ( fpvf + fnvf ) equal to 45.0% and 48.92% for 
 it . comparing the m3dls method to the other ones , we see automation as 
 a major advantage compared to the supervision required by fc 
 segmentation and stability of performance over clinical cases and 
 tissue types as a major advantage over the hmrf - em method .
the it 
 was set up with ideal parameters and can not be considered as a 
 potential segmentation method .
noncortical structures using two simultaneous level set functions limits the segmentation 
 process to the extraction of four homogeneous phases , 
 corresponding to four tissue types at most .
segmentation of 
 several other noncortical brain structures including the thalamus , 
 the caudate , the putamen , the palladium , the hippocampus , and the 
 amygdale is important for detection of neurological diseases such 
 as schizophrenia or alzheimer 's disease .
as shown in a recent 
 study for manual labeling of whole brain mri data , these 
 structures have greatly overlapping grayscale distributions which 
 explains the failure of simple thresholding methods to perform 
 efficient extraction of these individual more subtle structures . 
 on the other hand , as long as the histograms show reasonable 
 compactness ( with a small standard deviation ) , the 
 homogeneity - based proposed method should be capable of segmenting 
 them into separate phases . for the present study , we only focused 
 on cortical structures .
subcortical structures were therefore 
 removed prior to segmentation using the manually labeled data to 
 mask them , as illustrated in figure 3 . for future 
 applications of the method to data sets without the benefit of 
 manually labeled data , we can either use the labeling software 
 used by clinicians to remove these structures or use an in - house 
  intelligent manual segmentation tool  such as the one designed 
 by barrett et al . .
this paper presented a novel clinical application and quantitative 
 evaluation of a recently introduced multiphase level set 
 segmentation algorithm using t1-weighted brain mris .
the 
 segmentation algorithm performed an optimal partitioning of a 
 three - dimensional data set based on homogeneity measures that 
 naturally evolved toward the extraction of different tissue types 
 and cortical structures in the brain .
experimentation 
 on ten mri brain data sets showed that this optimal partitioning 
 successfully identified regions that accurately matched wm , gm , 
 and csf areas .
this suggests that by combining the segmentation 
 results with fiducial anatomical seed points , the method could 
 accurately extract individual structures from these tissues . 

random initialization was used to speed up the numerical 
 calculation and avoid the need for a priori information 
 input .
a regular initial partitioning of the data added some 
 robustness to the presence of local minima during the optimization 
 process .
comparison to three other segmentation methods was 
 performed with individual assessment of segmentation performance , 
 statistical comparison of the performance , and evaluation 
 of the statistical difference between the methods .
future work will include postprocessing of the segmented volumes 
 to extract individual structures such as the ventricles as well as 
 incorporation of available coregistered t2-weighted mri and pet 
 data to improve the segmentation performance by running the 
 algorithm on vectorial - type data . such extension of the method has 
 been proposed for color images but never for multimodality 
 clinical data .
 this paper presents the implementation and quantitative evaluation 
 of a multiphase three - dimensional deformable model in a level set 
 framework for automated segmentation of brain mris . 
 the 
 segmentation algorithm performs an optimal partitioning of 
 three - dimensional data based on homogeneity measures that 
 naturally evolves to the extraction of different tissue types in 
 the brain . 
 random seed initialization was used to minimize the 
 sensitivity of the method to initial conditions while avoiding the 
 need for a priori information . 
 this random initialization 
 ensures robustness of the method with respect to the 
 initialization and the minimization set up . 
 postprocessing 
 corrections with morphological operators were applied to refine 
 the details of the global segmentation method . 
 a clinical study 
 was performed on a database of 10 adult brain mri volumes to 
 compare the level set segmentation to three other methods : 
  idealized  intensity thresholding , fuzzy connectedness , and an 
 expectation maximization classification using hidden markov random 
 fields . 
 quantitative evaluation of segmentation accuracy was 
 performed with comparison to manual segmentation computing true 
 positive and false positive volume fractions . 
 a statistical 
 comparison of the segmentation methods was performed through a 
 wilcoxon analysis of these error rates and results showed very 
 high quality and stability of the multiphase three - dimensional 
 level set method . 
1. INTRODUCTION
2. METHOD
3. Experiments
4. RESULTS
5. DISCUSSION
6. CONCLUSION