AUTOMATIC BRAIN MR IMAGE SEGMENTATION BY RELATIVE

THRESHOLDING AND MORPHOLOGICAL IMAGE ANALYSIS

Kai Li, Allen D. Malony

Department Of Computer and Information Science University of Oregon

1202 University of Oregon, Eugene, OR 97403 USA

Don M. Tucker

Electrical Geodesics, Inc.

Riverfront Research Park, 1600 Millrace Drive, Suite 307, Eugene, OR 97403 USA

Keywords:

Segmentation, brain, MR, intensity inhomogeneity, relative thresholding, mathematical morphology, skeleton-

based opening, geodesic opening, a priori knowledge, ﬁrst-order logic.

Abstract:

We present an automatic method for segmentation of white matter, gray matter and cerebrospinal ﬂuid in T1-

weighted brain MR images. We model images in terms of spatial relationships between near voxels. Brain

tissue segmentation is ﬁrst performed with relative thresholding, a new segmentation mechanism which com-

pares two voxel intensities against a relative threshold. Relative thresholding considers structural, geometrical

and radiological a priori knowledge expressed in ﬁrst-order logic. It makes intensity inhomogeneity transpar-

ent, avoids using any form of regularization, and enables global searching for optimal solutions. We augment

relative thresholding mainly with a series of morphological operations that exploit a priori knowledge about

the shape and geometry of brain structures. Combination of relative thresholding and morphological opera-

tions dispenses with the prior skull stripping step. Parameters involved in the segmentation are selected based

on a priori knowledge and robust to inter-data variations.

1 INTRODUCTION

Image segmentation of brain magnetic resonance

(MR) images, particularly T1-weighted MR images,

has found wide applications ranging from morphome-

tric or functional studies to neurosurgical planning.

Some brain MR image segmentation methods aim

at labeling voxels with different brain tissue types.

Typically, the tissue types of interest are white mat-

ter (WM), gray matter (GM) and cerebrospinal ﬂuid

(CSF). Based upon the voxel classiﬁcation, we can re-

construct some particular brain structures, such as the

outer and inner cortical surfaces. Some segmentation

procedures directly aim at the extraction of cortical

surfaces. These two goals roughly correspond to two

different categories of dominant segmentation meth-

ods in brain image segmentation: voxel classiﬁcation

methods based on the intensity of image voxels and

deformable model-based methods considering struc-

tural information derived from the image.

Traditional voxel classiﬁcation methods such as

Expectation-Maximization (EM) and fuzzy c-means

clustering are problematic in the presence of inten-

sity inhomogeneity (IIH). IIH, also referred to as bias

ﬁeld, is an image artifact induced by a number of fac-

tors in image acquisition and perceived as a smooth

variation of intensity across the spatial domain of

the image. IIH correction is often performed as a

preprocessing step (Sled et al., 1998; Cohen et al.,

2000) or proceeds simultaneously with voxel classiﬁ-

cation (Pham and Prince, 1999; Ahmed et al., 2002;

Wells et al., 1998). Brain MR image segmentation

based on deformable models (MacDonald et al., 2000;

Zeng et al., 1999; Xu et al., 1999) tracks cortical sur-

faces by exploiting information derived from the im-

age, typically edges, together with a priori knowledge

about the location, geometry, and even shape of these

surfaces. These methods usually involve minimiza-

tion of an object function that attempts to ﬁt the model

to the image data with the regularization typically on

the surface smoothness.

Automation robustness of these segmentation

methods is not satisfactory yet mainly due to three

reasons. First, a common disadvantage of these meth-

ods is that their performance depends on the initial

segmentation state and is prone to be trapped in lo-

cal optima. Second, in those methods with regular-

ization involved, one or more parameters that control

the relative importance of regularization often have

to be empirically selected for each input image and

354

Li K., D. Malony A. and M. Tucker D. (2006).

AUTOMATIC BRAIN MR IMAGE SEGMENTATION BY RELATIVE THRESHOLDING AND MORPHOLOGICAL IMAGE ANALYSIS.

In Proceedings of the First International Conference on Computer Vision Theory and Applications, pages 354-361

DOI: 10.5220/0001366103540361

Copyright

c

SciTePress

may be sensitive to inter-data variations. This is the

case in the regularization of surface smoothness in de-

formable model-based methods and regularization of

the IIH variation smoothness in some voxel classiﬁca-

tion methods. Third, some brain tissue segmentation

methods require a prior step of skull stripping, which

by itself is a difﬁcult problem for complete automa-

tion (Rehm et al., 04).

W present a brain T1-weighted MR image seg-

mentation method using relative thresholding (RT)

and morphological operations that aims to improve

automation robustness in all three aspects described

above. Relative thresholding is based on image mod-

eling in terms of local spatial relationships between

near voxels and exploits structural, geometrical and

radiological a priori knowledge expressed in ﬁrst-

order logic. RT makes the IIH problem transpar-

ent, avoids using any form of regularization, and en-

ables global searching for optimal solutions. Results

from relative thresholding are improved mainly us-

ing a series of morphological operations. The major

two morphological operations are what we refer to as

skeleton-based opening and geodesic opening. They

are designed to robustly remove unwanted structures

from brain structures motivated by the a priori knowl-

edge about their special shape and geometry. Parame-

ters involved in the segmentation are selected based

on a priori knowledge and robust to inter-data varia-

tions. The combination of RT and morphological op-

erations dispenses with the prior skull stripping step.

The paper is organized as follows. We ﬁrst give

some basic deﬁnitions in section 2. RT and the two

morphological operations are presented in section 3

and 4 respectively. The whole segmentation pipeline

is described in section 5. The results are presented in

section 6 and the paper concludes in section 7.

2 DEFINITIONS

We ﬁrst deﬁne some terms used throughout the paper.

A 3D image can be viewed as a set of cubes struc-

tured regularly, where each cube represents a volu-

metric pixel (voxel). Each voxel v has three types of

neighbors: 6 face neighbors,12edge neighbors and 6

point neighbors that share a face, an edge, or a point

with v respectively.

The 6 face neighbors are regarded as connected

to the central voxel v in 6-connectivity and form the

6-neighborhood N

6

(v) of v. The 6 face neighbors

and the 12 edge neighbors form the 18-neighborhood

N

18

(v) (18-connectivity). Finally all 26 neighbors

form the 26-neighborhood N

26

(v) (26-connectivity).

Corresponding to the three types of connectivity,

three types of distance between two voxels, D

6

, D

18

,

and D

26

, are deﬁned as the number of steps in the

minimal path between the two voxels. Finally, we de-

ﬁne a grid graph from an image taking each voxel as

a vertex and adding edges in terms of one of the con-

nectivies between voxels.

3 RELATIVE THRESHOLDING

Relative thresholding is characterized as differentiat-

ing the labels of near voxels by comparing their inten-

sities with respect to a relative threshold. RT exploits

various a priori knowledge in terms of a critical data

structure which we refer to as gradient graph, and is

justiﬁed by image modeling based on the spatial con-

straints on the intensities of near voxels. Optimal rel-

ative thresholds are found with a trial-and-evaluation

scheme.

3.1 A Priori Knowledge

Let g = ∇g(σ

∇

) be the gradient vector image of

g(σ

∇

). Throughout this paper, we use g(σ) to denote

the result image of performing Gaussian ﬁltering with

standard deviation σ on the input image y. We con-

struct a directed graph G =(V,E) from g such that

each vectex v

i

∈ V corresponds to the voxel x

i

in a

region of interest R and each directed edge e

i

∈ E

emanates from v

i

to v

j

, where v

j

is the one of v

i

’s

26-neighbors that is in the direction of the gradient

vector g

i

. When v

j

is outside R, e

i

is forced to be a

loop from v

i

to itself.

The structural, geometrical and radiological a pri-

ori knowledge that we use in RT is:

• K

1

: CSF, GM, and WM are organized as a layered

structure from outside to inside;

• K

2

: The average intensities of CSF, GM, and WM

are in ascending order in T1-weighted MR images.

• K

3

: The cortex thickness is nearly uniform;

Based on the gradient graph, this a priori knowledge

is formulated as the following ﬁrst-order logic: There

exist a suitable σ

∇

and p such that:

• For each GM voxel v

i

, there is a path in G of length

p from v

i

to a WM voxel;

• For each CSF voxel v

i

adjacent to GM, there is a

path in G of length ≤ p from v

i

to a GM voxel;

• There is no path from a WM voxel to a non-WM

voxel in G;

• There is no path from any non-brain voxels to WM

in G without passing CSF.

3.2 Image Modeling

We model images in terms of the spatial relationships

between voxels instead of as statistical distributions

AUTOMATIC BRAIN MR IMAGE SEGMENTATION BY RELATIVE THRESHOLDING AND MORPHOLOGICAL

IMAGE ANALYSIS

355

on the absolute voxel intensities. Our intuition is

that if the segmentation task is not beyond the hu-

man recognition capability, near voxels of the same

type should possess less difference in intensity with

respect to that of near voxels of different types. With

this type of image modeling, we attempt to avoid the

limitation imposed by the form of statistical distribu-

tions and provide a framework for introducing various

a priori knowledge into the segmentation task.

Suppose that there are K voxel types among a total

of N voxels in the space domain Ω, which represents

either the whole image or a region of interest in the

image. In brain MR image segmentation, K typically

equals to 3 and the three tissues of interest are CSF,

GM and WM. Let the coordinates of voxels be x

i

, 1 ≤

i ≤ N, and the variable and true label of each voxel

respectively be ω

i

(or ω(x

i

))∈ [1,K] and ω

i

∈ [1,K],

1 ≤ i ≤ N. Incorporating a multiplicative bias ﬁeld

b

i

and an additive noise ρ, the image intensity y

i

(or

y(x

i

)), 1 ≤ i ≤ N, is modeled as:

y

i

= b

i

K

k=1

δ

k

i

y

k

i

+ ρ, where δ

k

i

=

0

ω

i

= k

1

ω

i

= k

(1)

In equation 1, δ

k

i

y

k

i

represents the component given

by tissue k in the ideal image without inﬂuence from

noise and IIH and we refer their sum

K

k=1

δ

k

i

y

k

i

as

the ideal image. Here we don’t assume any particular

statistical form on the noise term. Equation 1 is our

initial image model and will be gradually transformed

to facilitate image segmentation.

The term y

k

in equation 1 can be seen as an arbi-

trary function over the space domain governed by the

constraints on the spatial relationship between near

voxels. Generally, we think the constraints should

consider a priori knowledge about the structure and

geometry of the objects in the image as well as the

inherent image properties related to the image acqui-

sition process. In brain T1-weighted MR images, we

consider a priori knowledge K

1

and K

2

and use the

following ﬁrst-order logic to describe a spatial con-

straint: any pair of near voxels from two adjacent tis-

sues differ more than any pair of near voxels in one of

the two tissues:

∀x

i

,x

j

∈ Ω ∀k ∈ [1,K] ∃T

k

∈ [0, 1)

d(x

i

,x

j

) ≤ p ⇒

(

ω

i

= k ∧ ω

j

= k +1⇒ r(y

k

i

,y

k+1

j

) <T

k

) ∧

(ω

i

= ω

j

= k ⇒ r(y

k

i

,y

k

j

) > max(T

k−1

,T

k

)) (2)

r(a, b)=

a/b a < b

b/a a ≥ b

(3)

In equation 2, T

0

and T

K

are forced to be 0, d(x

i

,x

j

)

represents the distance between voxel x

i

and x

j

, and p

is a distance threshold (a voxel cube is of unit dimen-

sion). Theoretically, any form of distance, including

Euclidean distance, can be used. However, D

6

, D

18

or D

26

distance is preferable because of the computa-

tional efﬁciency.

A reasonable assumption about the bias ﬁeld is that

it varies slowly across the space with respect to the in-

tensity variation between different tissues in the ideal

image. We use a ﬁrst-order logic to describe this as-

sumption in equation 4 without any constraints on the

variation patterns.

∀x

i

,x

j

∈ Ω ∃ ∈ (0, 1)

(1 − 1 − max(T

1

, ..., T

K−1

)) ∧

(d(x

i

,x

j

) ≤ p ⇒ r(b

i

,b

j

) >) (4)

Based on the low frequency property of the bias ﬁeld,

we can safely let the function y

k

absorb the bias ﬁeld

function and the bias ﬁeld term can thus be dropped

from equation 1 while validity of the constraint in

equation 2 is maintained. Therefore the image arti-

fact of IIH is made transparent in our image model.

Next, we apply Gaussian ﬁltering on the original

gray level image to counteract the noise and drop the

noise term from equation 1. Let z = g(σ

z

) be a spe-

ciﬁc blurred image. The new image model on z is:

z

i

=

K

k=1

δ

k

i

z

k

i

(5)

Here, z

k

corresponds to the contribution of tissue k to

the blurred image. After Gaussian ﬁltering, we want

to maintain the spatial relationships between voxels,

as described equation 6.

∀x

i

,x

j

∈ Ω ∀k ∈ [1,K] ∃T

k

∈ [0, 1)

d(x

i

,x

j

) ≤ p ⇒

(

ω

i

= k ∧ ω

j

= k +1⇒ r(z

k

i

,z

k+1

j

) <T

k

) ∧

(

ω

i

= ω

j

= k ⇒ r(z

k

i

,z

k

j

) > max(T

k−1

,T

k

)) (6)

It is well-known that Gaussian ﬁltering blurs both

homogeneous regions and edges. This leads to mainly

two types of violation to the constraint. First, for a

voxel pair (x

i

,x

j

) of different types on the opposite

sides of an edge, if they are too close to each other,

r(z

i

,z

j

) may be signiﬁcantly increased such that they

may be identiﬁed as the same type. We think this ad-

verse effect can be minimized by increasing the dis-

tance between voxels pairs for comparison in the rel-

ative thresholding procedure. The second type of vi-

olation may occur when the dimension of some parts

of the structure of interest is too narrow compared to

the Gaussian ﬁlter aperture (σ

z

). We found that for

current MRI techiques, the usual resolution (around

1mm

3

) is high enough so that this violation brings

very little negative inﬂuence.

VISAPP 2006 - IMAGE ANALYSIS

356

3.3 Applying Relative Thresholds

Given two relative thresholds T

gw

between GM and

WM, and T

cg

between CSF and GM, we use two al-

gorithms for the segmentation. The inputs to both in-

clude a thresholding applying region R,acomparing

image h from which voxels are compared (we use z),

a gradient graph G constructed on R, and a relative

threshold. Based on the formulation of the a priori

knowledge, we found that it sufﬁces for the segmen-

tation to compare only voxels pairs along paths in the

gradient graph. Each comparison involves an objec-

tive voxel and a reference voxel. In algorithm 1, the

reference voxel ref

gw

(v

i

) for the object voxel v

i

is

of distance p from v

i

along the path emanates from

v

i

. In algorithm 2, the reference voxel ref

cg

(v

i

) for

the object voxel v

i

is the farthest GM vertex along the

path emanates from v

i

. We found that it gave slightly

better results to temporarily change h(ref

cg

(v

i

)) to

h(ref

gw

(ref

cg

(v

i

)))(1−(1−T

gw

)∗2) for each com-

parison in algorithm 2. This is based on our observa-

tion that 1 minus the optimal threshold T

∗

gw

is roughly

half of 1 minus the ratio between average GM and

WM intensities.

It is easy to see that as long as our image model

(formula 6) and the a priori knowledge in 3.1 are valid

and optimal thresholds are given, the ﬁrst algorithm

guarantees the exact separation between GM and WM

and the second algorithm guarantees the exact separa-

tion between CSF and GM. In practice, the segmenta-

tion accuracy depends on the goodness of the thresh-

olds and the extent to which the image model and the

a priori knowledge are valid.

Algorithm 1 GW-thresholding(R, h, G, T

gw

).

initialization: ∀x

i

∈ Rω

i

←WM;

for all x

i

∈ R do

if ω(ref

gw

(x

i

)) = GM then

ω

i

← GM;

else if r(h

i

,h(ref

gw

(x

i

))) <T

gw

then

ω

i

← GM;

end if

end for

Algorithm 2 CG-thresholding(R, h, G, T

cg

).

for all x

i

∈ R s.t. ω

i

= GM do

if r(h

i

,h(ref

cg

(x

i

)) <T

cg

then

ω

i

← CSF;

end if

end for

3.4 Finding Optimal Thresholds

The optimal thresholds are found in a trial-and-

evaluation scheme, as illustrated in algorithm 3. The-

oretically, the range for both thresholds are (0,1) and

there are inﬁnite candidate thresholds. But in prac-

tice, we found that it sufﬁces to try thresholds that are

multpilers of 0.01 in [0.75, 1) and [0.1, 1) for T

gw

and

T

cg

respectively. The range of T

gw

is selected based

on our observation that 1 minus the optimal thresh-

old T

∗

gw

is roughly half of 1 minus the ratio between

average GM and WM intensities, which is typically

greater than 0.5 in T1-weighted MR images.

The ﬁrst input to the algorithm is the region of

interest that contain signiﬁcant amount of all three

brain tissues and non-brain tissues as less as possi-

ble so that the disturbance to threshold searching is

minimized. We will refer to this region as threshold

searching region R

s

. Other inputs include a compar-

ing image h, used in the above two thresholding algo-

rithms, an evaluation image l, used for segmentation

evaluation, and a gradient graph G

s

, constructed on

R

s

. We choose to use the original intensity image y

as the evaluation image and the blurred image z as

the comparing image. After the optimal thresholds

are found, they are applied on a threshold applying

region R

a

with graph G

a

constructed on R

a

. R

s

and

R

a

are obtained with a method described in section 6.

Algorithm 3 Threshold-searching(R

s

, h, l, G

s

).

for all candidate threshold T

gw

do

apply GW-thresholding(R

s

, h, G

s

, T

gw

);

for all candidate threshold T

cg

do

apply CG-thresholding(R

s

, h, G

s

, T

cg

);

evaluate the segmentation, record the best

end for

end for

return the optimal threshold T

∗

gw

and T

∗

cg

;

Segmentation is evaluated with an objective func-

tion J as deﬁned in equation 7:

J =

x

i

,x

j

∈R

s

ϕ

ij

|l

i

− l

j

|/

x

i

,x

j

∈R

s

ϕ

ij

, (7)

where ϕ

ij

is 1 if there is a path of length 2 in G

s

from x

i

to x

j

,or0 otherwise. Minimization of J en-

courages intensity uniformity and penalizes intensity

difference between voxels of the same type.

It is helpful to understand the mechanism of rel-

ative thresholding by comparing the result from op-

timal thresholds and those from non-optimal thresh-

olds. Figure 1(b) shows that over-low threshold T

gw

globally makes the white matter too fat while ﬁgure

1(c) shows that over-high T

gw

globally makes it too

AUTOMATIC BRAIN MR IMAGE SEGMENTATION BY RELATIVE THRESHOLDING AND MORPHOLOGICAL

IMAGE ANALYSIS

357

thin. In either case, the GM/WM boundary drifts

away the correct situation in different directions.

(a) Input image (b) T

gw

=0.68, T

cg

=0.7

(c) T

gw

=0.93, T

cg

=0.7 (d) T

gw

=0.86, T

cg

=0.7

Figure 1: Analysis of relative thresholding applied on a

phantom image with optimal T

∗

gw

=0.86 and T

∗

cg

=0.7.

Due to various violation to our image model and a

priori knowledge, segmentation errors arise even with

optimal relative thresholds. A major error involves

labeling non-brain voxels as brain voxels which are

connected to the true brain voxels. This type of er-

ror will be corrected by morphological operations de-

scribed in the next section. Another typical error is

missing certain ﬁne WM structures which are recov-

ered with some ridge detection operations. Examples

of such errors are circled in ﬁgure 3(b). The whole

segmentation pipeline will be presented in section 5

4 MORPHOLOGICAL IMAGE

ANALYSIS

In this section, we present the skeleton-based open-

ing operation for eliminating false WM voxels and

the geodesic opening operation for eliminating false

GM and CSF voxels. Both operations are motivated

by the a priori knowledge about the shape and geom-

etry of brain structures. Before describing these two

operations, we need to ﬁrst have a review of geodesic

dilation and geodesic erosion.

Geodesic dilation(erosion) (Soille, 2003) involves

two input images: a marker image and a mask image.

The geodesic dilation(erosion) of size 1 of the marker

image f with respect to the mask image m is deﬁned

as the intersection(union) between m and the elemen-

tary dilation(erosion) of f. A usual constraint here is

that the mask image m is the superset(subset) of the

marker image f in geodesic dilation(erosion). The

geodesic dilation(erosion) of size n of a marker im-

age f with respect to a mask image m is obtained by

performing n successive geodesic dilations(erosions)

of size 1 of f with respect to m.

4.1 Skeleton-based Opening

Skeleton-based opening considers the characteristics

of the special surface-like shape of the WM and the

strong connectivity of the WM structure. It provides

a large parameter range that gives rise to acceptable

results. The main idea is that a special erosion and

geodesic dilation are performed on the surface skele-

ton of the object instead of on the object itself. The

algorithm consists of the following 5 steps.

1). Compute the reversible surface skeleton S of a

3D object by iterative thinning of the distance trans-

forms with the method described in (Svensson, 2001).

The skeleton is at most of two-unit thickness and is

further thinned to two unit-thickness skeletons S

1

and

S

2

, where S

1

∪S

2

= S, with the sub-iteration method

described in (Borgefors et al., 1999).

2). Iteratively erode S

1

and S

2

based on the

classiﬁcation of the skeleton points. We use the

method presented in (Saha and Chaudhuri, 1996)

classifying the skeleton points into the following

types: I(Isolated), SE(edge point of surfaces), S(inner

point of surface), SS(junction point of surfaces),

SC(junction point of surfaces and arcs), CE(arc end

point), C(inner point of arcs), CC(junction point of

arcs). In each iteration, a skeleton point is eroded only

if its type is I, SE, CC, CE, C, or CC. The number of

iterations n is taken as the parameter of the opening

operation.

3). Find the largest components, C

1

and C

2

, of the

eroded skeletons respectively.

4). Compute the geodesic dilation of size n of C

1

and C

2

with respect to S

1

and S

2

respectively. The

dilated skeletons are D

1

and D

2

respectively.

5). Recover the binary object from D

1

∪ D

2

by re-

verse distance transforms. Any voxels in the original

object but not in the new object are discarded.

Our experiments of applying SBO on a set of WM

groundtruths with resolution around 1mm

3

demon-

strated the following common behaviors of SBO:

there are two parameter threshold d

1

(around 13) and

d

2

(around 27) such that: 1) when n<d

2

, all fore-

brain WM voxels remain almost intact in the result

of SBO; 2) when n>d

1

the brainstem and cerebel-

lum (BC) get dropped. Therefore, we can see that

SBO provides two large parameter ranges for elimi-

nating false WM voxels ([1,d

1

)) and BC WM voxels

VISAPP 2006 - IMAGE ANALYSIS

358

([d

1

,d

2

)) as long as the connectivity of the WM struc-

ture maintains well in the prior segmentation results.

4.2 Geodesic Opening

Geodesic opening also involves a marker image f and

a mask image m ⊂ f. An example is that m denotes

the WM voxel set and f denote the union of GM voxel

set and m. Another example is that m denotes the

union of GM and WM while f denotes the union of

CSF, GM and WM. A geodesic opening of size n con-

sists of three steps: a geodesic erosion of size n of f

with respect to m resulting in image f

e

, identifying

the largest connected component f

c

of f

e

and elim-

inating any other components in f

e

, and a geodesic

dilation of size n of f

c

with respect to f.

Geodesic opening has the effect of smoothing the

geodesic distances of points in f − m from m by

eliminating protrusion voxels from f. Here, geodesic

distance of a point v in f − m from m is the min-

imum length of paths from v to m via only points

in f . This smoothing effect is particularly helpful

for elimination of unwanted structures from the ini-

tial GM and CSF voxel sets given by relative thresh-

olding due to the following characteristics about GM

and CSF structures: the thickness of human cortex is

nearly uniform and CSF ﬁlls concave spaces forming

a smooth boundary of the whole brain volume.

5 BRAIN MRI SEGMENTATION

In this section, we present the integration of RT

and morphological operations into a whole brain

MRI segmentation pipeline, as shown in ﬁgure 5.

The pipeline starts with binarization of the input

image and computing the threshold searching re-

gion (R

s

) and the threshold applying region (R

a

)

from the foreground. We simply apply traditional

thresholding with a threshold T

b

for the binariza-

tion. The background intensities are modeled as a

constant ﬁeld a plus a noise term in Rayleigh sta-

tistics. The background pdf is thus f (x|a, b)=

((x − a)/b

2

)e

−(x−a)

2

/2b

2

. The parameters of the

pdf are estimated based on what we refer to as au-

tomatic training, which consists of following three

steps. First, N

b

voxels of lowest intensity in g(σ

b

)

are selected as the training set. Second, the original

intensities of the voxels in the training set are used

as the sample data to estimate the parameter a

∗

and

b

∗

that give rise to the maximum likelihood. Third, a

threshold T

b

= a

∗

+3b

∗

is applied on y = g(0), g(δ),

g(2δ), ..., g(σ

b

= nδ) for the binarization resulting

in a new training set of background voxels. Here, a

voxel is taken as background if its intensity in any of

g(kδ)|

0≤k≤n

is less than T

b

. The second and third

step can be repeated multiple times to get improved

results.

Ridge

detection

Relative

thresholding

Graph cut

operation

Morphology

thresholding

training−base

d

Automatic

R

CSF segmentation

GM segmentation

Binarization

WM segmentation

Relative thresholding

andCompute R

a

S

Figure 2: Segmentation pipeline.

We obtain a head mask image and take it as the

threshold applying region R

a

by applying SBO and

morphology closing on the foreground. R

s

is ob-

tained by taking the intersection between the result of

applying morphology erosion on R

a

and those vox-

els within a distance d

1

s

from the superior end of the

head. The erosion kernel ball radius d

2

s

is selected

based on the a priori knowledge on skull thickness.

A threshold searching region is shown in ﬁgure 3(a).

A threshold applying region is shown in ﬁgure 3(b) as

all colorful areas.

After T

∗

gw

and T

∗

cg

are found on R

s

, apply them on

R

a

with the comparing image being z = g(σ

z

). Let

S

z

w

, S

z

g

, S

z

c

, denotes the result voxel set of WM, GM,

and CSF respectively. We also use them to denote the

corresponding mask images. An example of the result

in shown in ﬁgure 3(b). Note that some missing ﬁne

WM structure are marked with cyan circles. Next,

we use two ridge detection operations to recover these

missing ﬁne WM structures.

For each voxel x

i

∈ S

z

g

, the ﬁrst ridge detection

operation checks the 26-neighborhood of x

i

. Among

all neighbor pairs of opposite direction from x

i

, the

pair (x

j

,x

k

) with maximum value of |z

i

− z

j

| + |z

i

−

z

k

| is selected. If z

i

>z

j

and z

i

>z

k

, then x

i

is

taken as a WM ridge point, removed from S

z

g

, and

put into S

r

.Ifz

i

<z

j

and z

i

<z

j

, then x

i

is taken

as a valley point and put into S

v

. Figure 3(c) shows

the augmentation of WM ridge points.

The second operation depends on a masked

Gaussian ﬁltering. We use

g(σ, m) to denote the re-

sult image of applying masked Gaussian ﬁltering on

the input image y with standard deviation σ and mask

image m. Only voxels in the mask are ﬁltered and

any voxels outside the mask are not counted in the

convolution calculation. A masked Laplacian ﬁlter-

ing is then performed on

z = g(σ

z

,S

z

w

∪ S

z

g

− S

v

)

with mask S

z

g

−S

v

to recover more missing WM vox-

els. The Laplacian of each voxel x

i

∈ S

z

g

− S

v

is

x

j

∈N

6

(x

i

)∧x

j

∈S

z

w

∪S

z

g

−S

v

(z

j

− z

i

). If the laplacian

AUTOMATIC BRAIN MR IMAGE SEGMENTATION BY RELATIVE THRESHOLDING AND MORPHOLOGICAL

IMAGE ANALYSIS

359

value is less than 0, x

i

is taken as a WM ridge point,

removed from S

z

g

, and put into S

lr

. The Laplacian

ridge points are mainly used to enhance the connec-

tivity of the WM structure and may incorporate many

false WM voxels. They will be eliminated as the ﬁnal

step of WM segmentation.

We apply SBO on S

z

w

∪ S

r

∪ S

lr

with parameter

d

1

w

resulting in a WM set S

1

w

with non-brain vox-

els removed. The same operation is performed with

a greater parameter d

2

w

resulting in S

2

w

with further

removement of brainstem and cerebellum(BC). Then

we get the forebrain WM set as S

f

w

by subtracting the

largest component in S

1

w

− S

2

w

from S

1

w

and the BC

WM set as S

bc

w

= S

1

w

− S

f

w

. Result of applying SBOs

is illustrated in ﬁgure 3(d).

The separation between forebrain WM and BC

WM is calibrated by a graph-cut based method. First

we construct a grid graph in 26-connectivity from

S

1

w

. Then all vertices in S

bc

w

are contracted into a

source terminal and all vertices in S

f

w

that are of geo-

desic distance greater than a threshold d

f

from S

bc

w

are contracted into a sink terminal. Finally we com-

pute the minimum graph cut between the source and

sink terminal. It is demonstrated in (Boykov and Kol-

mogorov, 2003) that as long as the edge weights in

the grid graph are set to be inversely proportional to

their length, the graph cut can approximate the area of

the cut surface. In this way, the initial forebrain/BC

separation is substituted with a smaller and more reg-

ular cut surface and S

f

w

and S

bc

w

are updated accord-

ingly. Finally, we remove Laplacian ridge points by

subtracting S

lr

from S

f

w

. All discarded WM voxels

from S

z

w

that have a path in G

a

to S

f

w

or S

bc

w

are de-

graded to S

z

g

. Final result of WM segmentation is

shown in ﬁgure 3(e).

Given the result of WM, GM segmentation is made

easier considering the a priori knowledge on the

nearly uniform cortex thickness. First, we obtain an

initial forebrain GM set S

f

ig

⊂ S

z

g

by extracting all

voxels in S

z

g

that have a path to S

f

w

in G

a

of length

less than d

1

g

. In a similar way, we can get the BC GM

voxel set S

bc

g

and the total BC set S

bc

= S

bc

g

∪ S

bc

w

.

False GM voxels in S

f

g

are removed by a series of

geodesic openings with size being 1, 2, ..., d

1

g

. Let the

result be S

f

g

. These operations may over-smooth the

geodesic distances of GM voxels to the WM set. We

improve the result by simply computing the geodesic

dilation of size d

2

g

d

1

g

of S

f

g

with respect to S

f

ig

.

Any discarded WM voxels from S

f

ig

that have a path

in G

a

to S

f

w

are degraded to CSF and put into S

z

c

.

To improve the CSF voxel set given by relative

thresholding, we ﬁrst extract the subset S

c

⊂ S

z

c

in

which each voxels has a path in G

a

to S

f

g

∪ S

f

w

∪ S

bc

.

From S

c

we build a normal distribution for CSF voxel

intensities using a similar automatic training method

to that used for image binarization (see the ﬁrst para-

graph of this section). The distribution then guides

incorporation of additional CSF voxels. The up-

dated CSF voxel set S

c

is ﬁnally processed by a se-

ries of geodesic openings smoothing the brain volume

boundary. A example of ﬁnal result is shown in ﬁg-

ure 3(f) and two real image segmentation results are

shown in ﬁgure 4.

(a) Input image and R

s

(b) RT result on R

a

(c) First ridge detection (d) Applying SBO on WM

(e) Graph cut result (f) Final result

Figure 3: Segmentation of a phantom brain (IIH level =

40%, noise level = 9%). WM, GM, CSF, and BC are col-

ored with red, green, blue, and yellow respectively.

6 RESULTS

We have tested the segmentation algorithm against

a set of phantom brains as well as 6 real images

from 4 subjects scanned at 3 different MR machines.

Some examples of the results are shown in ﬁgure 4.

The resolution of the images are between 1.3mm

3

and 1mm

3

. The average running time is less than

8 minutes on a single 1.3GHz Power4 processor.

Searching for the optimal relative thresholds takes

about 1.5 minutes. We use exactly the same set of

parameters for all these images. Some important pa-

rameters are selected as follows: σ

z

=1, σ

∇

=2,

p =8, d

1

s

=95, d

2

s

=25, d

1

w

=4, d

2

w

=18, d

1

g

=15,

and d

2

g

=2. We think this set of parameters work

VISAPP 2006 - IMAGE ANALYSIS

360

for all input images with resolutions around 1mm

3

,

which is typical for current MR techniques. The only

place where user intervention may be needed is deter-

mining the transaxial dimension and the superior side

to obtain the threshold searching region. This infor-

mation can also be obtained by computation or from

the input ﬁle.

To evaluate segmentation, we use a metric over-

lap, same as the Tanimoto coefﬁcient (Duda and

Hart, 1973), comparing two segmentations for a given

voxel type labeling as the sum of voxels that both

have the same label in each segmentation divided by

the sum of voxels where either segmentation has the

label. Table 1 list the overlaps of WM and GM re-

spectively between the groundtruth and our results

on phantom brains corrupted with partial volume ef-

fect, various levels of IIH and noise (Cocosco et al.,

1997). The noise level represents the percent ratio

of the standard deviation of the white Gaussian noise

versus the signal for a reference tissue. For a IIH level

x, the multiplicative IIH ﬁeld has a range of values of

1 − x...1+x over the brain area.

Figure 4: Results of two real images.

7 CONCLUSION

We have demonstrated a new automatic brain MRI

segmentation method with relative thresholding and

morphological operations. Relative thresholding

makes IIH transparent and allows global searching

for optimal solutions. Prior skull stripping is not re-

quired. Parameters are selected based on a priori

knowledge and are robust to inter-data variations.

Table 1: WM and GM overlaps on a set of phantom brains.

Type IIH N=3% N=5% N=7% N=9%

WM 20% .8931 .8899 .8897 .8812

WM 40% .8926 .8899 .8832 .8824

GM 20% .8100 .8043 .8060 .7967

GM 40% .8074 .8060 .8048 .7996

REFERENCES

Ahmed, M. N. et al. (2002). A modiﬁed fuzzy c-means

algorithms for bias ﬁeld estimation and segmentation

of mri data. IEEE Trans. Medical Imaging, 21(3).

Borgefors, G. et al. (1999). Computing skeletons in threee

dimension. Pattern recognition, 32(7).

Boykov, Y. and Kolmogorov, V. (2003). Computing geo-

desics and minimal surfaces via graph cuts. In Inter-

national Conference on Computer Vision.

Cocosco, C. A. et al. (1997). Brainweb: Online interface to

a 3d mri simulated brain database. NeuroImage, 5(4).

Cohen, M. S. et al. (2000). Rapid and effective correction

of rf inhomogeneity for high ﬁeld magnetic resonance

imaging. Human Brain Mapping, 10.

Duda, R. and Hart, P. (1973). Pattern Classiﬁcation and

Scene Analysis. Wiley, New York.

MacDonald, D. et al. (2000). Automated 3-d extraction of

inner and outer surfaces of cerebral cortex from mri.

NeuroImage, 12(3).

Pham, D. L. and Prince, J. L. (1999). An adaptive fuzzy c-

means algorithm for image segmentation in the pres-

ence of intensity inhomogeneities. Patt. Rec. Let.

Rehm, K. et al. (04). Putting our heads together: a consen-

sus approach to brain/non-brain segmentation in t1-

weighted mr volumes. NeuroImage, 22.

Saha, P. K. and Chaudhuri, B. B. (1996). 3d digital topology

under binary transformation with applications. Com-

puter vision and image understanding, 63(3).

Sled, J. G. et al. (1998). A nonparametric methods for au-

tomatic correction of intensity nonuniformity in mri

data. IEEE Trans. Medical Imaging, 17.

Soille, P. (2003). Morphological image analysis. Springer-

Verlag.

Svensson, S. (2001). Reversible surface skeletons of 3d ob-

jects by iterative thinning of distance transforms. In

G.Bertrand et al. (Eds.): Digital and Image geometry,

LNCS 2243.

Wells, W. M. et al. (1998). Adaptive segmentation of mri

data. IEEE Trans. Medical Imaging, 15.

Xu, C. Y. et al. (1999). Reconstruction of the human cere-

bral cortex from magnetic resonance images. IEEE

Trans. Medical Imaging, 18.

Zeng, X. et al. (1999). Segmentation and measurement of

the cortex from 3-d mr images using coupled-surfaces

propagation. IEEE Trans. Medical Imaging, 18(10).

AUTOMATIC BRAIN MR IMAGE SEGMENTATION BY RELATIVE THRESHOLDING AND MORPHOLOGICAL

IMAGE ANALYSIS

361