DTIC ADA465665: A Statistically Based Surface Evolution Method for Medical Image Segmentation: Presentation and Validation

9 Pages·0.19 MB·English

by  Defense Technical Information Center

A Statistically Based Surface Evolution Method
for Medical Image Segmentation: Presentation
and Validation
Eric Pichon1 Allen Tannenbaum1 and Ron Kikinis2
1 Georgia Instituteof Technology, AtlantaGA 30332, USA
feric, [email protected],
http://users.ece.gatech.edu/~eric
2 HarvardMedical School, Boston, MA02115, USA
[email protected]
Abstract. Inthispaperwepresentanewalgorithmfor3Dmedicalim-
agesegmentation. Thealgorithmisfast, relativelysimpletoimplement,
and semi-automatic. It is based on minimizing a global energy de(cid:12)ned
from a learned non-parametric estimation of the statistics of the region
to be segmented. Implementation details are discussed and source code
is freely available as part of the 3D Slicer project. In addition, a new
uni(cid:12)ed set of validation metrics is proposed. Results on arti(cid:12)cial and
real MRI images show that the algorithm performs well on large brain
structures bothin terms of accuracy and robustnessto noise.
1 Introduction
The problem of segmentation, that is (cid:12)nding regions in an image that are ho-
mogeneous in a certain sense, is central to the (cid:12)eld of computer vision. Medical
applications, visualization and quanti(cid:12)cation methods for computer-aided diag-
nosis or therapy planning from various modalities typically involve the segmen-
tation of anatomical structures as a preliminary step.
In this paper wewill considerthe problem of (cid:12)nding the boundariesof onlyone
anatomical region with limited user interaction. Interactivity is very desirable
since the user will be given the opportunity to make use of often implicit but
absolutelynecessaryexternalknowledgetoguidethe algorithmtowardsaresult
that would make sense for her task. The segmentation process can be repeated
in order to identify as many di(cid:11)erent regions as necessary.
Manydi(cid:11)erentapproacheshavebeenproposedtoaddressthesegmentationprob-
lem which can be dually considered as (cid:12)nding regions or (cid:12)nding boundaries.
Focusing only on the boundaries is less complex computationally but also less
robust since information inside the region is discarded. Typically this is the
approach of the snake and active contours variational methods [7,16,17].
While the original region-growing algorithm [11] formalism is extremely crude,
interesting extensionshavebeen proposedin [9] where somestatistical informa-
tionisderivedfromtheregionasitexpands.Thesetechniqueshavebeenapplied
to medical image analysis [12,14]. The relation between region-growingand ac-
tive contours has been studied in [15] and more recently active contours have
been extended to an elegant active regions formalism [8] where the boundaries
of regions are deformed accordingto an evolution equation derived to minimize
an energy based on some statistics of the regions.
Report Documentation Page Form Approved
OMB No. 0704-0188
Public reporting burden for the collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information,
including suggestions for reducing this burden, to Washington Headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington
VA 22202-4302. Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to a penalty for failing to comply with a collection of information if it

does not display a currently valid OMB control number.
1. REPORT DATE 3. DATES COVERED
2003 2. REPORT TYPE 00-00-2003 to 00-00-2003
4. TITLE AND SUBTITLE 5a. CONTRACT NUMBER
A Statistically Based Surface Evolution Method for Medical Image
5b. GRANT NUMBER
Segmentation: Presentation and Validation
5c. PROGRAM ELEMENT NUMBER
6. AUTHOR(S) 5d. PROJECT NUMBER
5e. TASK NUMBER
5f. WORK UNIT NUMBER
7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 8. PERFORMING ORGANIZATION
Georgia Institute of Technology,School of Electrical and Computer REPORT NUMBER
Engineering,Atlanta,GA,30332
9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR’S ACRONYM(S)
11. SPONSOR/MONITOR’S REPORT
NUMBER(S)
12. DISTRIBUTION/AVAILABILITY STATEMENT
Approved for public release; distribution unlimited
13. SUPPLEMENTARY NOTES
14. ABSTRACT
15. SUBJECT TERMS
16. SECURITY CLASSIFICATION OF: 17. LIMITATION OF 18. NUMBER 19a. NAME OF
ABSTRACT OF PAGES RESPONSIBLE PERSON

a. REPORT b. ABSTRACT c. THIS PAGE 8
unclassified unclassified unclassified
Standard Form 298 (Rev. 8-98)
Prescribed by ANSI Std Z39-18
2 Basic Flow
In this section, we state the fundamental (cid:13)ow underpinning the segmentation
method.Let(cid:10) beanopenconnectedboundedsubsetofRn withsmoothbound-
ary @(cid:10). Let t : (cid:10) !Rn be a family of embeddings, such that 0 is the iden-
tity. Let (cid:30) : Rn ! R be a positive C1 function. We set (cid:10)(t) := t((cid:10)) and
S(t):= t(@(cid:10)). We consider the family of (cid:30)-weighted volumes
H(t):= (cid:30)( t(x))d t(x)= (cid:30)(y)dy:
Z Z
(cid:10) (cid:10)(t)
Set X = @ tj : then using the area formula [6] and then the divergence theo-
@t t=0
rem, the (cid:12)rst variation is dHj = div((cid:30)X) dx:=(cid:0) ((cid:30)X)(cid:1)N dy; where
dt t=0 (cid:10) @(cid:10)
N istheinwardunitnormalto@(cid:10).RConsequentlythecorRresponding(cid:30)-weighted
volume minimizing (cid:13)ow is
@S
=(cid:30)N:
@t
A di(cid:11)erent derivation of the same result has previously been proposed in [4].
3 Method
In what follows we will only consider the 3D case. A region R will be a subset
of R3 with smooth boundary S = @R: As above, N denotes the corresponding
inward unit normal vector to S.
Given an image I, a non-negative weighting function w((cid:1);(cid:1)) and a region R we
de(cid:12)ne the energy:
E(I;w;R):= w( I(x); krI(x)k ) dx: (1)
Z
R
E isthe weightedvolumeof the regionR.Theweightofavoxelxisdetermined
by the function w((cid:1);(cid:1)) of the local properties I(x) and k rI(x) k of the image.
Ideally, w should re(cid:13)ect the local properties of the region we want to segment.
As this is not known a priori we will heuristically estimate w as we evolve R to
maximize E.
Proposition 1. Notation as above. Then for a given weighting function w, the
evolution in which the energy E(I;w;R) is decreasing as fast as possible (using
@S
only local information) is =wN.
@t
Proof. Follows immediately from the discussion in Section 2. tu
Since w is a non-negative function, the (cid:13)ow is reversible. In particular, the (cid:13)ow
in the reverse direction,
@S
=(cid:0)wN; (2)
@t
gives the direction in which the energy is increasing as fast as possible (using
local information). In the context of segmentation, one may think of (2) as a
bubble and of the original (cid:13)ow as a snake.
GivenanapproximationR oftheregiontobesegmentedwecanuseamaximum
0
likelihood-like approach to determine the weighting function w which would a
0
posteriori justify the segmentation of R .
0
Proposition 2. For a given (cid:12)xed region R , the energy E(I;w;R ) is maxi-
0 0
mized by setting w to p the conditional probability on that region:
0
w =argmax E(I;w;R )=Pr( I(x); krI(x)k j x2R ): (3)
0 0 0
p
Proof. We can rewrite the energy as:
E(I;p;R )= N (u;v):w(u;v) dudv;
0 Z Z R0
I krIk
whereN (u;v)isthevolumeofthesetofpointsx2R suchthatI(x)=uand
R0 0
krI(x)k=v.Butthisisjustaconstantmultiple ofPr( I(x); krI(x)k jx2
R ) which by the Schwartz’s inequality is the maximizer of E. tu
0
As the region evolves, p will be periodically updated according to (3). This will
change the de(cid:12)nition of the energy (1) and therefore(2) can only be considered
a gradient (cid:13)ow for every time interval when w is (cid:12)xed.
4 Implementation
Weimplemented ourmethodasamoduleof the open-sourcesoftware3DSlicer.
It is freely available for download at http://www.slicer.org.
4.1 Surface evolution
As the (cid:13)ow (2) is unidirectional (the surface can only expand since w (cid:21)0) any
voxel x will eventually be reached at a time T(x). Knowing T is equivalent to
knowing R or S since by construction:
R(t)=f x; T(x)(cid:20)t g and S(t)=@R(t): (4)
Solvingthe (cid:13)ow(2)forS(t)isequivalenttosolvingthe Eikonalequation(5)for
T(x):
krT(x)k(cid:1)w(x)=1: (5)
This can be done very e(cid:14)ciently using the Fast Marching method [3]. Starting
from known seed points which de(cid:12)ne the initial surface, the algorithm marches
outwards by considering neighboring voxels and iteratively computing arrival
times T in increasing order. The seed points are set by the user inside the
structure to be segmented. By construction, when computing T(x), the surface
containsthevoxelxaswellasallvoxelsforwhichT hasalreadybeencomputed.
The algorithm will terminate when T is known for all points and using (4) we
willknowS(t)foranyt.Wewillthen letthe userdeterminewhattime t of the
0
evolution corresponds best to the region she wants.
Notethatwith averydi(cid:11)erentformalism,ourmethodis,initsimplementation,
very reminiscent of region growing. For example, the min-heap data structure
which makes Fast Marching e(cid:14)cient is the direct equivalent of the sequentially
sortedlistintheseededregiongrowingalgorithm[9].Infactouralgorithmcould
be made a direct non-parametric extension of seeded region growing simply
by arti(cid:12)cially forcing arrival times to zero for all points inside the surface S.
Relationsbetween regiongrowingand variationalschemeshavebeen previously
exposed in [15].
4.2 Estimation of probability density function
The probability has been modi(cid:12)ed to p = p (m)(cid:1)p (h) where M and H are
M H
the median and interquartile range (the di(cid:11)erence of between the (cid:12)rst and last
quartile) operators on a 3(cid:2)3(cid:2)3 neighborhood. M and H convey more or less
the sameinformationasI (graylevel)andkrI k(localhomogeneity)but their
non-linear nature makes them more robust to noise and allow them to respect
better the edges of the image.
WeuseParzenwindows[10]toestimatetheprobabilitydensityfunctions.Itisa
non-parametrictechnique and therefore no assumption is requiredon the shape
ofthedistributions.Givenawindowfunction(cid:30)andN samplesm ;:::;m and
1 N
h ;:::;h the densities are estimated by:
1 N
N N
1 1
p (m)= (cid:30)(m(cid:0)m ) and p (h)= (cid:30)(h(cid:0)h )
M i H i
N N
X X
i=1 i=1
5 Validation
Objective and quantitative analysis of performance is absolutely crucial (but
often overlooked) when proposing a segmentation algorithm. Since designing a
segmentation method is challenging (lack of unifying formalism, high diversity
in the applications, subjectivity, implicitness, etc.) it does not come as a sur-
prise that the validation of such an algorithm will also be challenging. Di(cid:11)erent
methods have been studied (see [20] and references therein). We will propose
a unifying framework for discrepancy measures based on the number and the
position of mis-segmented voxels and show how it relates to classical measures.
We will then apply it to the validation of segmentation of realistic synthetic
images(forwhichthe \groundtruth"i.e.perfectsegmentationis known)atdif-
ferentlevelsofnoiseforaccuracyandrobustnessassessmentaswellastomanual
expert segmentation of real datasets.
5.1 Classical discrepancy measures
Di(cid:11)erentmeasureshavebeenproposedtoassesstheresemblancebetweenapro-
posedsegmentationS andthe correspondinggroundtruthG.TheDiceSimilar-
ity Coe(cid:14)cient has been widely used and it can be derived as an approximation
of the kappa statistic (see [1]). It is de(cid:12)ned as:
V(S\G)
DSC(S;G):=
1(V(S)+V(G))
2
Where V((cid:1)) is the volume (number of voxels) of a set.
Onedisadvantageofthiscoe(cid:14)cientisthatitonlytakesintoaccountthenumber
ofmis-segmentedpixelsanddisregardstheirpositionandthereforetheseverities
of errors.This wascorrected in Yasno(cid:11)’s normalized discrepancymeasure(ND,
see [18]) and the Factor of Merit (FOM, see [5]):
N N
1 1 1
ND:= d(i)2 and FOM :=
N e N 1+d(i)2
X X
i=1 i=1
Where N is the number of mis-segmented voxelsand d(i) is the erroron the ith
voxel. Another popular measure is the Hausdor(cid:11) distance:
H(S;G):=maxf maxminks(cid:0)g k; maxminks(cid:0)g k g
s2S g2G g2G s2S
H(S;G) isthe maximumdistancewewouldhavetomovethe boundariesof one
set so that it would encompass completely the other set. As this is extremely
sensitivetoextremeerrors,thepartialHausdor(cid:11)distanceH (S;G)canbeintro-
f
duced(see[2])asthe maximumdistancewewouldhavetomovethe boundaries
of one set so that it would cover f% of the other set.
5.2 Proposed framework
Consider now the error-distance:
0 for x correctly segmented (x2S\G)
8
minkx(cid:0)sk for x under-segmented (x2GnS)
d(x):=>><s2S
minkx(cid:0)gk for x over-segmented(x2SnG)
g2G
Assumingthat>>:allpointsx2S[G areequally likelyd canbe seenasa random
variable D which describes completely the discrepancy between S and G. We
can study D using the standard statistical tools:
probability of error: PE:=Pr(D >0)
mean error: (cid:22) :=mean(D j D >0)
D>0
standard deviation of error: (cid:27) :=stdev(D j D>0)
D>0
partial distance-error: D :=f (cid:0)quantile(D)
f
These measures receive a natural intuitive interpretation.
{ PE is the probability for a voxel x 2S\G to be misclassi(cid:12)ed (either over-
or under-segmented).
{ An erroneous voxel is on average (cid:22) pixels o(cid:11). This value is or is not
D>0
typical depending on the standard deviation (cid:27) .
D>0
{ D is the error distance of the worst f% voxels or equivalently the max-
1(cid:0)f
imum distance we would need to move erroneous voxels for the error to be
improved to PE=f.
As an example, PE = 10%; (cid:22) = 3:1; (cid:27) = 0:3 and D = 14 would
D>0 D>0 0:99
meanthattheoverlapbetweenthegroundtruthandtheproposedsegmentation
is90%.The10%remaningpixelsareeither under-segmentedorover-segmented
pixels (\false positive" i.e. pixels that are in S and not in G). On averagethese
pixelsare3.1pixelso(cid:11).Thisvalueisverytypicalsincethe standardvariationis
low(0.3).Howeverthereisnoreasonforthe errortobe Gaussianand,here, the
tail probability is not negligible since the worst 1% pixels are at least 14 pixels
o(cid:11). This could be due to a thin, long (cid:12)nger of mis-segmented pixels.
The following proposition justi(cid:12)es the de(cid:12)nition of these new uni(cid:12)ed measures.
Proposition 3. ThesemeasuresarerelatedtothemeasurespresentedinSection
5.1 according to:
DSC
1(cid:0)DSC(cid:20)PE=(1(cid:0)DSC)=(1(cid:0) ) (6)
2
1
(cid:0)1(cid:20)((cid:22)2 +(cid:27)2 )=ND (7)
FOM D>0 D>0
e
H1(cid:0)f=(1(cid:0)PE) (cid:20)D1(cid:0)f (cid:20)H1(cid:0)f2 (8)
Proof. in future publication (in particular, D1 =H) tu
5.3 Results on simulated datasets
The publicly available Brain Web [19] datasets have been generated from a
knowngroundtruthusingasophisticatedphysicalmodeling[13]oftheMRIpro-
cess. We can assessin a perfectly objective way the performance of our method
by comparing the result of our segmentation with the underlying ground truth.
Note that even though these datasets are computer-generated they are very re-
alistic (see (cid:12)gure 1(b)) Another interesting aspect of this project is that from
the same ground truth, datasets with di(cid:11)erent levels of noise can be simulated
whichallowsustostudytherobustnessofourmethodwithrespecttonoise.We
segmented the lateralventricle, white matter (WM) and white matter and gray
matter (WM+GM) on 2 datasets:
{ Normal brain, T1, 1(cid:2)1(cid:2)1 mm (181(cid:2)181(cid:2)217 voxels), 3% noise, 20%
intensity non-uniformity ("RF") (standard parameters of the Brain Web
model).
{ Normalbrain, T1,1(cid:2)1(cid:2)1mm (181(cid:2)181(cid:2)217 voxels), 9%, 40%(highest
levels of noise available).
Our results (Table 1) show that the proposed algorithm gives very good results
on these structures (DSC > 0:7 has been described as a good agreement in
the literature, see for example [1]). The complex structure of the white matter
makes it more challenging and explains the somewhat mediocre performance
(in the case of the maximum noise dataset, the cerebellum was not perfectly
segmented). In the highest level of noise, connectivity between the lateral and
the thirdventricleswaslost(theintraventricularforamenof Monrodisappeared
in the noise). This increased the strength of the ventricle edges in the noisy
dataset and, paradoxically, simpli(cid:12)ed the segmentation. Overall the algorithm
appears extremely robust to noise.
DSC PE (cid:22)D>0 (cid:27)D>0 D0:95 D0:99
Ventricle 92.0% 95.1% 14.9% 9.4% 1.07 1.13 0.48 0.61 1.00 1.00 1.00 1.41
WM 91.9% 80.3% 15.0% 32.0% 1.59 2.03 1.58 1.94 1.00 2.83 3.61 8.25
WM+GM 96.2% 95.2% 7.4% 9.2% 1.42 1.40 1.25 1.15 1.00 1.00 1.41 2.00
Table 1. Performance measure on arti(cid:12)cial dataset. Left bold, with standard noise,
right,withmaximumnoise.Underlinedresultsareillustratedby(cid:12)gures1(b),1(d),1(f).
5.4 Results on real datasets
In this real case, the pathological diagnoses are meningiomas (brain tumor).
Patients’ heads were imaged in the sagittal and axial plane with a 1.5 T MRI
system3 withapostcontrast3Dsagittalspoiledgradientrecalled(SPGR)acqui-
sition with contiguous slices. The resolution is 0:975(cid:2)0:975(cid:2)1:5 mm (256(cid:2)
256(cid:2)124 voxels). These datasets were manually segmented by one expert.
Because of inter- and intra-expert variabilitywe should expect these results not
tobeasgoodasinthe syntheticcase.Itshouldalsobenotedthatthearbitrary
conventionsofthemanualsegmentationsareresponsibleforalotoftheobserved
error since for example the ventricle was labeled as gray matter, the medulla
oblongata and the spinal cord have been left out etc. (compare Fig. 1(a) and
1(c)). Overall, nonetheless, results are consistent with the arti(cid:12)cial case.
DSC PE (cid:22)D>0 (cid:27)D>0 D0:95 D0:99
Tumor 78.0% 88.0% 36.0% 21.4% 1.97 1.34 1.63 0.94 3.32 1.41 7.00 2.83
WM+GM 96.1% 92.4% 7.5% 14.2% 1.69 1.28 1.99 0.75 1.00 1.00 2.00 2.24
Table 2. Performance measure on 2 real datasets. Underlined results are illustrated
by(cid:12)gures 1(a),1(c),1(e).
3 Signa, GE Medical Systems,Milwaukee, WI.
(a) Sagittal slice of real (b)Axialsliceofnoisyar-
dataset and proposed ti(cid:12)cial dataset and pro-
white and gray matter posed ventricle segmenta-
segmentation (white) tion (white)
(c) Expert segmentation (d) Underlying ground
(gray)andproposedwhite truth (gray) and pro-
and gray matter segmen- posed ventricle segmenta-
tation (white) tion (white)
(e) Rendered surface of (f) Rendered surface of
proposed white and gray proposed ventricle seg-
matter segmentation mentation
6 Conclusion
Wepresentedanewcurveevolution(cid:13)owbasedonlearnednon-parametricstatis-
tics of the image. Implementation is straightforwardand e(cid:14)cient usingthe Fast
Marching algorithm and is freely available as part of the 3D Slicer project. An
extensivevalidationstudyaswellasanewuni(cid:12)edsetof validationmetricshave
also been proposed.
Future work will focus on extending our formalism into a purely variational
framework, adding some regularizing constraints and extending the validation
study.
Acknowledgements
EricPichonandAllenTannenbaumaresupportedbyNSF,NIH,AFOSR,MRI-
HEL and ARO.
RonKikinisissupportedbygrantsPO1CA67165,R01EB000304andP41RR13218.
References
1. ZijdenbosA.,DawantB.,andMargolinR. Morphometricanalysisofwhitematter
lesions inMR images: Method and validation. IEEE TMI,13(4):716{724, 1994.
2. HuttenlocherD.,KlandermanG., andRucklidgeW. Comparing images usingthe
Hausdor(cid:11)distance. PAMI,15(9):850{863, 1993.
3. Sethian J. Level Set Methods and Fast Marching Methods. Cambridge University
Press, 1999.
4. SiddiqiK.,LauziereY.,TannenbaumA.,andZuckerS.Areaandlengthminimizing
(cid:13)ows for shape segmentation. IEEE TMI,7:433{443, 1998.
5. StrastersK.andGerbrandsJ.Three-dimensionalsegmentationusingasplit,merge
andgroup approach. Pattern Recognition Letters, 12:307{325, 1991.
6. Simon L. Lectures on geometric measure theory. In Proceedings of the Centre for
Mathematical Analysis, Australian National University, Canberra, 1983.
7. Kass M., Witkin A., and Terzopoulos D. Snakes: Active contour models. Int. J.
Computer Vision, 1:321{332, 1988.
8. Paragios N.andDericheR. Geodesicactiveregions:Anewparadigmtodealwith
frame partition problems in computer vision. Journal of Visual Communication
and Image Representation, 13:249{268, 2002.
9. AdamsR. andBischof L. Seeded region growing. PAMI,16(6):641{647, 1994.
10. DudaR.,Hart P., andStorkD. Pattern Classi(cid:12)cation. Wiley-Interscience, 2001.
11. Gonzalez R.and Woods R. Digital Image Processing. Prentice Hall, 2001.
12. Justice R., Stokely E., Strobel J., Ideker R., and Smith W. Medical image seg-
mentation using 3-D seeded region growing. Proc. SPIE Symposium on Medical
Imaging Volume, 3034:900{910, 1997.
13. Kwan R., Evans A., and Pike G. MRI simulation-based evaluation of image-
processing andclassi(cid:12)cation methods. IEEE TMI,18(11):1085{1097, 1999.
14. Pohle R. and Toennies K. Segmentation of medical images using adaptive region
growing. In Proc. SPIE Medical Imaging.
15. Zhu S. and Yuille A. Region competition: Unifying snakes, region growing, and
bayes/MDLfor multibandimage segmentation. PAMI, 18(9):884{900, 1996.
16. McInerney T. and Terzopoulos D. Deformable models in medical image analysis:
Asurvey. Medical Image Analysis, 1(2):91{108, 1996.
17. Caselles V., Kimmel R., and Sapiro G. Geodesic active contours. In Proc. ICCV,
pages 694{699, 1995.
18. Yasno(cid:11)W.,MiuJ.,andBacusJ. Errormeasuresforscene segmentation. Pattern
Recognition, 9:217{231, 1977.
19. WorldWide Web. http://www.bic.mni.mcgill.ca/brainweb/.
20. ZhangY. Asurveyonevaluationmethodsforimagesegmentation. Pattern Recog-
nition, 29(8):1335{1346, 1996.