A thesis presented
by
Yi-Ting Wang
to
Institute of Computer Science and Engineering
College of Computer Science
in partial fulfillment of the requirements for the degree of
Master
in the subject of
Computer Science
National Chiao Tung University Hsinchu, Taiwan
We proposed a multiresolutional brain extraction framework which utilize graph cuts technique to classify head magnetic resonance (MR) images into brain and non-brain re-gions. Our goal is to achieve both high sensitivity and specificity results of brain extrac-tion. Started with an extracted brain with high sensitivity and low specificity, we refine the segmentation results by trimming non-brain regions in a coarse-to-fine manner. The extracted brain at the coarser level will be propagated to the finer level to estimate fore-ground/background seeds as constraints. The short-cut problem of graph cuts is reduced by the proposed pre-determined foreground from the coarser level. In order to consider the impact of the intensity inhomogeneities, we estimated the intensity distribution locally by partitioning volume images of each resolution into different numbers of smaller cubes. The graph cuts method is individually applied for each cube. The proposed method was compared to four existing methods, Brain Surface Extractor, Brain Extraction Tool, Hybrid Watershed algorithm, and ISTRIP, by using four data sets, the first and the second IBSR data set of the Internet Brain Segmentation Repository, BrainWeb phantom images from the Montreal Neurological Institute, and healthy subjects collected by Taipei Veterans General Hospital. The performance evaluation for brain extraction, our method outperforms others for the first/second IBSR data set and BrainWeb phantom data set, and performs compara-bly with the BET and ISTRIP methods when using the VGHTPE data set.
List of Figures v
List of Tables vii
1 Introduction 1
1.1 Brain extraction . . . 2
1.1.1 Related works of brain extraction . . . 4
1.2 Max-flow and min-cut . . . 9
1.2.1 Ford and Fulkerson max-flow algorithm . . . 9
1.2.2 Graph cuts . . . 11
1.3 Motivation of this work . . . 17
2 Multiresolutional graph cuts for brain extraction 19 2.1 Flowchart of multiresolutional brain extraction . . . 20
2.2 Removal of bright voxels in non-WM regions . . . 22
2.3 Graph cuts for brain extraction . . . 25
2.3.1 Edge weights in the graph cuts method . . . 25
2.3.2 Multiresolutional graph cuts for brain extraction . . . 30
2.4 Postprocessing to fill holes that dark brain region be classified as back-ground by graph cuts . . . 33
3 Results for brain extraction 41 3.1 Image data sets for performance evaluation . . . 42
3.2 Evaluation metrics for brain extraction . . . 43 iii
1.1 Anatomical image from left to right are original image, GM, WM and CSF 3
1.2 An example of brain extraction . . . 3
1.3 The demonstration of intensity inhomogeneity . . . 3
1.4 Illustration of watershed . . . 5
1.5 Illustration of pre-flooded watershed . . . 5
1.6 Illustration of pre-flooding . . . 6
1.7 Temporal and cerebellum in T1-weighted MRI . . . 7
1.8 The skull stripping stages of BSE . . . 8
1.9 Illustration of HWA . . . 8
1.10 The example of Ford and Fulkerson max-flow algorithm . . . 10
1.11 Graph cuts . . . 12
1.12 Short-cut problem of graph cuts . . . 13
1.13 Lazy snapping . . . 14
1.14 Video snapping cut . . . 15
1.15 Multilevel banded graph cuts . . . 16
1.16 Illustration the influence of different threshold values on the quality of ini-tial mask for skull stripping using graph cuts . . . 17
2.1 Flowchart of the proposed brain extraction method . . . 21
2.2 An example of multilevel image pyramid. . . 22
2.3 The four categories by Otsu’s method of the initial brain extraction . . . 24
2.4 Demonstration of the start points at the center region for 3-D region grow-ing to rough segment WM . . . 25
2.12 Demonstration dark volume that need to be post-processed . . . 38 2.13 Illustration of the post-processing . . . 39 3.1 The discrete anatomical model of BrainWeb phantom image. . . 44 3.2 The example of brain volumes with apparent intensity inhomogeneity and
artifacts in the first IBSR data set. . . 47 3.3 The extracted brain images of subject in the first IBSR data set and using
the proposed method. . . 52 3.4 The extracted brain images of subject in the second IBSR data set and using
the proposed method. . . 53 3.5 The extracted brain images of subject in the VGHTPE data set and using
the proposed method. . . 54 3.6 The extracted brain images of subject in the BrainWeb phantom data set
and using the proposed method. . . 55 3.7 The result of extracted brain without using distance term at the highest level. 56 3.8 The result of extracted brain without using pre-determined foreground factor. 57 3.9 The result of extracted brain which perform graph cuts algorithm in a
slice-by-slice manner. . . 57 4.1 An example of the extracted brain with excluded cerebellum voxels by HWA. 60 4.2 An example of the possible disadvantage by using localized graph cuts. . . 61 4.3 An example of the possible disadvantage by using the distance term in the
3.1 Performance evaluation using the first IBSR data set and excluding the fail-ure results of each brain extraction method . . . 48 3.2 Performance evaluation using the first IBSR data set without including all
the failure cases . . . 48 3.3 Performance evaluation using the second IBSR data set . . . 49 3.4 Performance evaluation using the VGHTPE data set . . . 49 3.5 Performance evaluation using the BrainWeb phantom data set and only
considering WM and GM . . . 50 3.6 Performance evaluation using the BrainWeb phantom data set and
consid-ering WM, GM and CSF . . . 50 3.7 Performance evaluation of the proposed method without distance term/predetermined
term and in a slice-by-slice manner using the second IBSR data set . . . 51
processing. WM consists of myelinated axons and is responsible for delivering message between different brain regions of GM. The function of CSF, which is a transparent bodily fluid inside the brain and ventricle, is to protect the brain.
Accurate brain extraction is important to voxel-based morphology (VBM) [1] and neu-roimaging studies. VBM of magnetic resonance (MR) image is a powerful and non-invasive method to investigate human brain in vivo, it compares the concentration of GM or WM regionally in cross-sectional (the same time between groups) or longitudinal (the same group at different times) manner to describe ageing [2] or diseases, such as schizophrenia, bipolar disorder, and Alzheimer’s [3] [4]. The VBM procedure can be separated into fol-lowing steps: brain extraction, intensity inhomogeneity correction, tissue segmentation, registration, and voxel based analysis. Brain extraction is the work of excluding non-brain tissues such as skull, muscle, and fat, as shown in Figure 1.2. Accurate brain extraction can improve the above procedures. Intensity inhomogeneity correction is to correct the intensity variation due to magnetic inhomogeneity. Theoretically, the intensity of the same tissue in brain should be similar. However, the same tissue may brighten or darken at dif-ferent parts because of magnetic inhomogeneity, as shown in Figure 1.3, and may affect intensity based method in the following sequence. Tissue segmentation is to separate WM, GM and CSF. Registration is to warp each subject into the same stereotactic space. In this way, VBM can analysis the difference of specific tissue between inter-subjects in the same space. The accuracy of above processes is important to neuroimaging studies, such as cortical thickness analysis [5]. Intensity inhomogeneity correction, tissue segmentation, registration, and voxel based analysis are all benefit from brain extraction.
Figure 1.1: Anatomical image from left to right are original image, GM, WM and CSF.
Figure 1.2: An example of brain extraction.
Figure 1.3: The demonstration of intensity inhomogeneity. Intensity of GM and WM in red window is dark than other parts.
Intensity-based approaches
Intensity-based approaches use the threshold of brain/non-brain to achieve brain extrac-tion. This method is sensitive to noise or intensity inhomogeneity.
Watershed (WAT) assumes that brain anatomy is the connectivity of WM and WM is rounded by GM and CSF. This method describes the high intensities as hill and the low intensities as valley, as shown in left of Figure 1.4. If a valley between two regions then they are separated. But watershed may results in over-segmentation caused by noise that separate brain into too many regions. To solve over-segmentation, it’s necessary to merge those regions that belong to the same structure by postprocess [7]. Pre-flooded WAT [7] proposed a pre-flooding approach to tackle the problem of over-segmentation. They inverted the original data, as shown in Figure 1.5, and merged basin with height less than certain intensity value, Figure 1.6 show the pre-flooding manner. The result is shown in the right of Figure1.4. The disadvantage of this method is that sometimes the results contain some eye skull or exclude parts of the cerebellum.
Figure 1.4: Illustration of watershed. The gray value is transferred to height information. Left: the original image. Right: the skull stripped with pre-flooding manner. (The figure is cited from [7].)
Figure 1.5: Illustration of pre-flooded watershed. The gray level of Figure 1.4 is inverted. Hills in the Figure 1.4 are now represented by basins. Left: axial view. Right: sagittal view. (The figure is cited from [7].)
Figure 1.6: Illustration of pre-flooding. A basin is merged into a deeper basin, if and only if its depth relative to the current voxel intensity is less than or equal to the preflooding height. (The figure is cited from [9].)
Edge-based approaches
Edge-based approaches use gradient to locate the boundary between brain/non-brain, but may cause some errors in temporal since there are strong gradient to separate temporal and cerebellum, Figure 1.7 show the location of temporal and cerebellum).
Brain Surface Extractor (BSE) [8] uses anisotropic diffusion filtering to smooth noise, Marr-Hildreth edge detector to produce closed contours, and refines the result by mathe-matical morphology. The steps of BSE shown in Figure 1.8.
Deformable approaches
Deformable approaches is a method of brain extraction by deforming a rough surface to the brain surface.
temporal cerebellum
Figure 1.7: Temporal and Cerebellum in T1-weighted MRI. The region of yellow narrow is temporal. The region of green narrow is cerebellum.
BET [6] finds a brain/non-brain threshold, initial a triangular tessellation of a sphere’s surface at center of gravity (COG), and finally iteratively deforms to the brain surface.
ISTRIP [25] initializes an ellipsoid by finding the effective intensity range, then uses Wendland’s radial basis function (RBFs) deforming to the brain surface.
Hybrid approaches
Hybrid approaches combine methods of different types to segment brain/non-brain. For the purpose of correcting the method of WAT, HWA combines pre-flooded WAT and deformable surface model to get better accuracy of brain extraction. HWA first perform the pre-flooded WAT to obtain an rough brain extraction and this brain as initial of deformable surface model, then warp a accurately segmented brain atlas to the result of the pre-flooded WAT by rigid registration and deform to the accurate surface in the end, as shown in Figure 1.9. HWA is benefits of robust by using the pre-flooded WAT and geometric constraints of deformable surface model but is too conservative and still containing lots of non-brain.
Figure 1.8: BSE: Skull Stripping Stages. (a) A slice from the initial volume. (b) The same slice after anisotropic diffusion filtering three iterations with diffusion parameter 15. (c) The edge map following application of the MarrVHildreth operator. (d) The same slice from the extracted brain, following morphological processing of the edge map. (The figure is cited from [8].)
Figure 1.9: Illustration of HWA. The yellow surface is the initial template deformed to the segmented volume. The red line is the result after using deformable surface model. (The figure is cited from [9].)
1.2
Max-flow and min-cut
Max-flow min-cut theorem means that the value of max-flow is equal to the value of min-cut in the flow network. A flow network G = (V, E) is a graph that each edge (u, v)E has capacity c(u, v) ≥ 0 and with two specific vertices: a source s and a sink t. We can image each edge as a tube and compute maximum flow from source s to sink t per unit time. The min-cut presents the cut cross the edges whose summation is the minimum capacity to separate source s and sink t. Removing edges crossed by min-cut lead to no flow can pass from source s to sink t.
1.2.1
Ford and Fulkerson max-flow algorithm
Ford and Fulkerson proposed a method to find the value of max-flow and meanwhile find a min-cut. Below is the algorithm and Figure1.10 demonstrates the Ford and Fulkerson method.
Algorithm Ford and Fulkerson
Inputs Graph G with flow capacity c, a source node s, and a sink node t Output A flow f from s to t which is a maximum
1.f (u, v) ← 0 for all edges (u, v)
2.While there is a path p from s to t in Gf, such that cf(u, v) > 0 for all edges
1.Find cf(p) = min{cf(u, v) | (u, v)p}
2.For each edge (u, v)p
1.f (u, v) ← (u, v) + cf (Send flow along the path)
2.f (v, u) ← (v, u) − cf(p) (The flow might be ”returned” later)
1 Augmenting amount: 10 Initial flow : 0 Augmenting amount: 10 Augmenting amount: 10 10 s v2 v3 v1 t 50 50 30 10 30 10 10 s v2 v3 v1 t 10/50 50 30 10/10 30 10 10 s v2 v3 v1 t 20/50 50 10/30 10/10 30 10/10 10/10 s v2 v3 v1 t 20/50 10/50 20/30 10/10 30 10/10 Min cut Max flow=30
1.2.2
Graph cuts
Graph cuts is a method for separating foreground and background and is similar to max-flow min-cut theorem. The difference is that in graph cuts, each vertices has an edge (t-link) connected to source S, an edge connected to sink T and an edge (n-link) between neighboring vertices.
When we use graph cuts, the first step is to select foreground seeds and background seeds. And then we assign edge weights as capacity. The weight of n-link is according the difference of gray value or color between neighboring points. The weight of t-link is according the similarity to foreground seeds and background seeds. If the intensity of a node is like to foreground seeds than the weight of t-link connected to source S is big that means this node has more chance be classified as foreground and vise versa. Finally, according to the weights of t-link and n-link, we can get a cut to partition foreground and background by graph cuts algorithm. Figure1.11 demonstrates the graph cuts.
Short-cut is a well-known problem in graph cuts and mostly occur in thin and long region. Short-cut problem means that the cost of desired boundary is bigger than the short-cut path that misclassification happened. Figure1.12 is an example of the short-short-cut problem. This problem can be solved by combining geodesic information and graph cuts algorithm [12].
Graph cut in general image and video
Based on [11], many applications of graph cuts were published. Whatever no methods are perfect, therefore, user interactive tools often come with graph cuts [13][14].
Lazy Snapping [13], the picture is pre-segmented by watershed and user only need to use mouse to stroke off about foreground and background regions to define the seeds. Because the pre-segmented picture takes a watershed region as a node, therefore, nodes in
Figure 1.11: An example of graph cuts. (a) An Image with foreground/background seeds. (b) The constructed graph for graph cuts. The blue lines are t-links and the green lines are n-links. The thickness presents the weight of each link. (c) According to the weight of t-link and n-link, graph cut algorithm will find an optimal path to segment foreground and background. (d) The classified image.
Foreground seeds Background seeds
S
2S
1 Goal of boundary S1 Short-cut S2Figure 1.12: Short-cut problem of graph cuts. S1 is the goal of boundary and S2 is the
short-cut path. Although the edge weights of n-links on the goal of boundary boundary S1
are smaller than those on the short-cut path S2, however, the boundary S1 is to long cause
the cost is bigger than short-cut S2. Therefore, graph cuts may result in short-cut S2instead
Figure 1.13: Lazy Sanpping. (a)The original image. (b)User roughly define the foreground seeds (yellow line) and background seeds (blue line). (c) The blue polygon is the instant visual feedback and have UI tool to edit. (d) Synthesis the graph cuts result and another picture. (The figure is cited from [13].)
graph cuts is reduced. So Lazy Snapping can offer instant visual feedback and polygon tool to refine the result, as shown in Figure1.13. Furthermore, [15], [16], and [17] are also the segmentation methods using graph cuts and watershed.
Video Snap Cut [18] mentioned that graph cuts with global statistics can’t handle com-plexity scene therefore they segment with localized classifiers. The user will roughly define the foreground and background on key frame, and get a initial result. The localized clas-sifiers then surround along the initial result. Each localized clasclas-sifiers will combine local information to get a segmentation, as shown in Figure1.14. These localized classifiers will propagate to next frame based on optical flow and update their local statistic. For refining the result of graph cuts, [16] also provides a tool to compute local color distribution.
Graph cuts in medical image
Multilevel banded graph cuts (BGC) [19] and [21] are two graph cuts applications in medical image. Multilevel banded graph cuts (BGC) [19] run full grid graph cuts with the coarsest image and get an initial object , then refine the initial object in a coarse-to-fine
Figure 1.14: Video Snapping Cut. Illustration of local classifiers. (a) Overlapping clas-sifiers (yellow squares) are initialized along the object boundary (red curve) on frame t. (b) These classifiers are then propagated onto the next frame by motion estimation. (c) Each classifier contains a local color model and shape model, they are initialized on frame t and updated on frame t+1. (d) Local classification results are then combined to generate a global foreground probability map. (e) The final segmented foreground object on frame t + 1. Original video courtesy of Artbeats. (The figure is cited from [18].)
Figure 1.15: Multilevel Banded Graph Cuts.
manner and graph cuts only operate in a band which between eroded and dilated the object, as shown in Figure1.15. In this way, we can decrease memory usage and time cost but might exclude some thin structures such as blood vessels showed in the experiment of 3-D heart and pulmonary artery CT volume. To tackle the problem of BGC, [20] presented a accurate multilevel BGC by identifying thin structures and adding to BGC system when fining the lower resolution object.
In brain extraction topic, [21] first finding white matter (WM) by region growing and regarding as foreground seeds. Second, decide a threshold mask as background seeds. key point is that the threshold need to cause narrow connection between brain and non-brain tissue. For the graph cuts result, if the threshold is too low, it will lead to residual non-brain. On the contrary, choosing a too high value will exclude brain volume, Figure 1.16
Figure 1.16: Illustration the influence of different threshold values on the quality of initial mask for skull stripping using graph cuts. Too low threshold (second column) leads to in-sufficient separation between brain and non-brain structures, high threshold (right column) results in brain loss. (The figure is cited from [21].)
show the selection of threshold. The disadvantage of this method is the threshold mask may contains some partial brain voxels, hence after apply graph cuts algorithm the method performing closing operation (10mm voxel dilation and 10mm voxel erosion), however we always don’t know how much brain volume is included in background seeds and how much volume need to be recovered.
1.3
Motivation of this work
In the existing brain extraction methods, some have high sensitivities whereas some have high specificities. The high sensitivity stands for extracting large percent of the brain region, and the high specificity value rejects large percent of the non-brain region. The pursuits of high sensitivity and high specificity are generally conflicting. How to comprise between high sensitivity and high specificity is essential to accurate brain extraction. HWA method has the advantage of high sensitivity, but the disadvantage is its low specificity. In this work, we planned to develop an accurate brain extraction algorithm by using the brain region estimated by HWA method as the initial. While maintaining high sensitivity, our algorithm achieves high specificity by trimming non-brain regions from the initials with
ground seeds. The proposed multiresolutional graph cuts method extracts brain volume in a coarse-to-fine manner, The extracted brain at coarser level will propagate to finer level and as a constraint to estimate foreground seeds. To measure local histogram of brain tissue intensities, we partition volume images of each resolutions to different number of smaller cubes. The graph cuts method is individually performed for each cube at the same level and the results of these cubes are integrated to obtain the whole brain volume. Furthermore, we use a pre-determined foreground form coarser level to tackle the short-cut problem.
Multiresolutional graph cuts for brain
extraction
voxels brighter than the segmented WM tissue. This volume is denoted by VL, where VL
represents that this volume is at level L, in this thesis. The deleted voxels are usually non-brain voxels, such as fat.
The second step is the multiresolutional graph cuts and the advantage of multiresolu-tional framework to this work is that we can use the result of coarser levels as constraints to finer levels. We downsample VLimage into lower-resolution ones, and from fine to coarse
levels are VL−1, VL−2, . . . , V0, as shown in Figure 2.2. The eight voxels in L − l level (a
2 ∗ 2 ∗ 2 cube) are combined into one voxel in L − (l + 1) level, l = 0, 1, . . . , L − 1. We start the multiresolutional graph cuts from the coarsest level. The shape of brain in this level is much easier to obtain, because of the brighter voxels not belonging to brain are blurred. To prevent short-cut problem in graph cut, the best foreground seeds should cover the object as much as possible. Hence, we take the WM as our foreground seeds at the level 0. And the estimation of foreground/background seeds at finer levels are constrained by upper level extraction results.
At the final step, we fill holes, which are background areas surrounded by foreground areas in coronal view. Since the darker voxels in the ventricles might be classified as back-ground by the graph cuts algorithm and the deleted voxels in the first step might include some voxels of WM, and resulting in scattered holes.
There will be more details of removing the bright voxels in non-WM regions, multires-olutional graph cut and the final step in the following sections for the above steps.
Figure 2.2: The multilevel image pyramid.
2.2
Removal of bright voxels in non-WM regions
The foreground seeds are determined by intensity distribution, which may be affected by the voxels with high intensities. Therefore, before application of graph cuts algorithm, we delete those bright voxels according to the estimated WM intensity distribution. We use 3-D region growing and Otsu’s [24] method to estimate the WM tissue in the way similar to [21]. The other reason of removing bright voxels in non-WM regions is that they may be classified as foreground by graph cuts algorithm, because the intensities in non-WM brain regions should not be larger than those of WM regions.
To estimate WM, in our observation, we can partition the results of rough brain extrac-tion into four categories from darker to brighter with three thresholds which were found using Otsu’s method. The first category is composed of voxels with gray values like CSF and non-head voxels as shown in Figure 2.3 (c). The second category is composed of vox-els with gray valuess between those of CSF and GM as shown in Figure 2.3 (d). The third
category is majorly composed of voxels with gray values of GM as shown in Figure 2.3 (e). The fourth category is majorly composed of voxels with gray values of WM as shown in Figure 2.3 (f). When gathering statistics of the distribution of gray values, we only consider voxels with gray values lower than T , where T is the gray value higher than 98 percent of total nonzero voxels (similar to [6]). The reason to exclude those voxels with particularly high values is that they may dominate the third threshold of Otsu’s [24] method.
For each slice, we apply the Otsu’s method described above to determine the three thresholds and the largest threshold is used to segment WM voxels. We start 3-D region growing to find WM with voxels higher than the threshold between WM and GM for each slices located between 2/5 and 3/5 middle transverse slices. These starting points of region growing locate at the center region of each slices as shown in Figure 2.4. The center region at superior and inferior transverse slices may be non-brain region, therefore we limit start points between 2/5 and 3/5 middle transverse slices. Since the connectivity of 3-D region growing is higher than that of 2-D region growing, we use a strict threshold t as the stopping criterion of region growing, to avoid from growing WM region to brighter non-brain region. Thus threshold t represents the maximum intensity difference between 3-D neighboring voxels. To improve the accuracy of WM segmentation, we perform 2-D region growing starting from the results of 3D region growing, and uses the third threshold of Otsu’s method and T as the stopping criterion. Figure 2.5 demonstrates the results of WM segmentation.
To avoid the influence of bright voxels in non-WM region on the intensity distribution, we remove those voxels with gray values larger than the average of estimated WM sub-tracted by the estimated WM standard deviation. We call this volume VL, as shown in the middle row of Figure 2.6. The reason why we remove those voxels with high gray values is that those voxels most likely belong to a non-brain region, such as fat . These kind of voxels are determined as foregrounds while doing the graph cut. Removal of the brighter voxels in the non-WM areas can improve the estimation of WM as well as the brain extraction
(c) (d) (e) (f)
Figure 2.3: The four categories by Otsu’s method of the initial brain extraction. (a) The result of initial brain extraction. (b) The histogram of single slice and can be classified into four category by threshold 1, threshold 2 and threshold 3. The horizontal axis is the gray value and the vertical axis is the number of voxels. (c) The first category is below threshold 1 and composed of voxels with gray values like CSF and non-head voxels. This image is shown the nonzero voxels with gray value below threshold 1. (d) The second category is between threshold 1 and threshold 2 and composed of voxels with gray values between those of CSF and GM. (e) The third category is between threshold 2 and threshold 3 and composed of voxels with gray values of GM. (f) The fourth category is greater than threshold 3 and composed of voxels with gray values of WM.
Figure 2.4: Voxels overlaid with warm colors are the starting points at the center region for 3-D region growing. The background image is the result of rough brain extraction.
result, and the result can be more accurate.
2.3
Graph cuts for brain extraction
2.3.1
Edge weights in the graph cuts method
The proposed method construct a graph at each level and apply the graph cuts algorithm to extract the brain region. In the constructed graph, neighboring voxels have a connected edge, which is denoted by n-link. Besides, each voxel has edges, denoted by t-link, con-nected to two specific vertices: foreground terminal S and background terminal T as shown in Figure 2.7. There are two kinds of n-links that need to be assigned with weights, one is the inter-n-link and the other is the intra-n-link. Inter-n-link connects voxels in consecutive slices, shown as black lines in Figure 2.7, and the intra-n-link connects voxels in the same slice, shown as orange lines in Figure 2.7. For segmentation of foreground and background, we need to assign edge weight for each edge in the constructed graph and apply graph cut
Figure 2.5: The results of segmented white matter tissue in coronal slice. The subject is in the second IBSR data set. Top row are the results of the initial brain extraction. Bottom row are the results of the segmented white matter tissue.
Figure 2.6: The result of deleting non-WM voxels with gray values higher than estimated WM. Top row: the results of initial brain extraction. Bottom row: the result of deleting the voxels with gray values higher than estimated WM.
of the graph cut, and A(xi, xj)/B(xi) is the function to assign the weight of n-link/t-link,
respectively.
The finest level has different equation of edge weights from other levels. To assign the weight of n-link at coarser levels, the edge weight is assigned as:
A(xi, xj) = wn× exp(−
(Ixi − Ixj)
2
2α2 ) , (2.2)
where wnis a parameter representing the importance of n-link, Ixi and Ixj are the intensity
of two adjacent voxels, xi and xj, connected by this link. If the difference between Ixi and
Ixj is large, then the weight of this link becomes small. Thus, the cut path is more likely to
go through this link.
The weight of n-link at the finest level is assigned as:
A(xi, xj) = wn× exp(− (Ixi− Ixj) 2 2α2 ) + wd× ( D(xi) + D(xj) 2 ) . (2.3)
In the finest level, we add a distance term D(x) to Eq. (2.2). The D(x) is defined as the least-squares distance from x to the contour of the extracted brain which is upsampled from the coarser level to the finest level. In this work we use a library which is based on kd-trees and box-decomposition trees for approximate nearest neighbor searching (ANN) [23] to estimate the distance term.
The parameter value of wn and wd for inter slice is half of the value in intra slice, in
For the weight of t-link connected to foreground terminal is defined as:
B(xi) = wt×
PF(xi)
PF(xi) + PB(xi)
, (2.4)
where wt stands for the importance of this term in contrast to n-link. For each voxel x,
PF(xi) and PB(xi) are its likelihood functions to foreground seeds and background seeds,
respectively: PF(xi) = 1 √ 2πσf exp(−(Ixi− µf) 2 2σ2 f ) , (2.5) PB(xi) = 1 √ 2πσb exp(−(Ixi − µb) 2 2σ2b ) . (2.6)
The σf is the standard deviation of foreground seeds, whereas and σb is set as twice the
standard deviation of background seeds for increasing the separability of foreground and background by increasing the background range. The µf and µb are the average of
fore-ground seeds and backfore-ground seeds. If voxel Ixi is more likely a background seed, then
the t-link connected to the foreground terminal is much easier to be cut, then x has more chance to be classified as a background. The weights of foreground seeds connecting to foreground terminal are set to be infinity and they will be classified as foreground. Simi-larly, t-link connected to background terminal is defined as:
B(xi) = wt×
PB(xi)
PF(xi) + PB(xi)
. (2.7)
The weights of background seeds connecting to background terminal are set to be infin-ity and they will be classified as background. The edge weight calculation functions from Eq. (2.1) to Eq. (2.7) were modified from those in [11]. We added a distance term in Eq. (2.3) at the original resolution level. We use the graph cuts algorithm proposed in [22] and the implementation can be found in http://www.cs.ucl.ac.uk/staff/V.Kolm
Figure 2.7: Demonstration of n-link and t-link. Blue lines are the t-link connecting to background terminal. Green lines are the t-link connecting to foreground terminal. Orange lines are the intra slice n-link (link in the same slice). Black lines are the inter slice n-link (link in different slices).
2.3.2
Multiresolutional graph cuts for brain extraction
For the purpose of multiresolutional graph cut, we downsample the VL(the highest res-olution) at level L into coaser levels, VL−1, VL−2, . . . , V0. The equation of downsampling
is defined as: Vx,y,zL−l = 1 8 1 X i=0 1 X j=0 1 X k=0 V2x+i,2y+j,2z+kL−(l−1) , for 1 ≤ l ≤ L . (2.8) The Vx,y,zL−l is the intensity at coordinate (x, y, z) and at level L − l. The number of voxels at L − (l − 1) level are 8 times of voxels at L − l level.
Graph cuts at the coarsest level
To get the shape of brain, we start running graph cut with V0, which is at the lowest resolution. At this level, the shape of brain is much easier to obtain, because of the brighter voxels that are not belonging to brain were blurred. On the other hand, the brighter voxels
(WM) belonging to the brain area are set as foreground seeds, which always be classified as foreground. We use the region growing method, describe in Section 2.1, to segment the WM tissue in V0 as foreground seeds, as voxels overlaid with warm colors in the
middle row of Figure 2.8. On the other hand, we use the voxels with non-zero gray values and below the lowest threshold as background seeds, as shown as voxels overlaid with cool colors in the middle row of Figure 2.8. After the graph cut of the entire brain was completed, we get a foreground object, O0, as show in bottom row of Figure 2.8. Then
erode O0 one voxel to get a smaller volume, E0, as a constraint to the middle resolution level.
Graph cuts at the middle levels
The graph cuts at middle levels is to refine the extracted brain from the lowest resolution level. At middle dimension level l, 0 < l < L, to get a pre-determined foreground area we first upsample the eroded brain region, El−1, at coarser level l − 1 to level l the pre-determined foreground area El. The upsampling equation of the eroded volume is defined
as: Ex,y,z0l = El−1x 2, y 2, z 2 , (2.9)
where Ex,y,z0l is the intensity at coordinate (x, y, z) and at level l. The weights of t-link connecting to foreground terminal in E0l, as shown in third row of Figure 2.9, are set to be infinity. This ensures E0l must be classified as foreground. The pre-determined foreground area can reduce the short-cut problem and prevent cut at WM/GM boundary.
Considering local gray value distribution, we partition Vl into 2l× 2l× 2lcubes (each
dimension is divided into smaller overlapping region, as shown in Figure 2.10. Then we perform graph cuts algorithm for each cube individually. Each cube has 1/4 in length overlapped with its adjacent cubes, this reduce the inconsistent border between two cubes. Each cube has its own foreground and background seeds. We estimate the three thresholds
is defined as El. For the overlapped region between neighboring cubes, we integrate the
results by logical OR operation.
Graph cuts at the finest level
Graph cuts applied at the finest level is to get the final extracted brain. At the finest level, the pre-determined foreground area is the upsampled EL−1, as shown in the third row of Figure 2.11. We partition VLinto 2L× 2L× 2Lcubes individually. Then we perform graph
cut with each cube. The way of finding foreground seeds and background seeds is the same as in the middle levels, as shown in the second row of Figure 2.11. We also upsampling OL−1to this level, called O0L, as shown in top row cool colors of Figure2.11. The equation of upsampling OL−1is defined as:
O0Lx,y,z = P1 k=−1O L−1 x 2, y 2, z 2+k 3 , (2.10)
where Ox,y,z0L is the intensity at coordinate (x, y, z) and at level L. Different from middle and the coarsest levels, we add a distance term: wd× (D(xi)+D(x2 j)), in Eq. (2.3). The distance
term D(x) is defined as the least-squares distance from x to contour of O0L. The reason is that, the brighter voxels belonging to non-brain region are similar to foreground seeds, so they might be classified as foreground. But in the middle and the coarsest levels, those brighter voxels are blurred, hence those voxels can be judged as background seeds. There-fore, we use the contour of O0L as a factor of the edge weight in this level because when the n-link is closer to the contour, it should has more chance to be cut. For the overlapped
region between neighboring cubes, we integrate the results by logical OR operation.
2.4
Postprocessing to fill holes that dark brain region be
classified as background by graph cuts
The final step is to fill holes which are the dark brain areas classified as background regions by graph cuts algorithm, such as ventricles, as shown in Figure 2.12. The voxels surrounded by foreground regions, as shown in the yellow circle of Figure 2.13 also need to be filled, because these voxels were wrongly deleted in Step 1 of the system flowchart, since they failed to be included in WM region growing process due to their irregular intensities. The postprocessing is illustrated in Figure 2.13. We find the contour of the extracted brain from the highest resolution level, including the removed dark brain region, as shown in Figure 2.13 (b). Then we take the outermost contour, as shown in Figure 2.13 (c), and incorporate the voxels inside the contour to fill dark brain regions. Figure 2.13 (d) shows the result of postprocessing. Notice that this postprocessing is applied on coronal slices.
Figure 2.8: Demonstration of foreground seeds, background seeds and results at the lowest resolution level. The volumes are in the second IBSR data set. Top row: the volumes are V0, individually. Middle row: the warm colors are foreground seeds and cool colors are background seeds. Bottom row: the extracted brain at the lowest resolution level.
Figure 2.9: Demonstration of foreground seeds, background seeds and results at middle levels. Top row: the volume at middle levels. Second row: the warm color is foreground seeds and cool color is background seeds. Third row: the pre-determined foreground which is eroded and upsampled the result at the lowest resolution level. Bottom row: the extracted brain at middle levels.
Figure 2.10: Demonstration the unit of graph cut at level 1. The region surrounded by yellow, purple, green, and blue dotted line are the unit of graph cut, individually. Each unit perform graph cut once. The left image of this figure is the view in 2-D and the right one is the view in 3-D.
Figure 2.11: Demonstration of foreground seeds, background seeds and results at the finest level. Top row: the volume at the highest level and the cool color is contour of O0L. Second row: the warm color is foreground seed and cool color is background seed. Third row: the pre-determined foreground which is eroded and upsampled the result at the middle level. Bottom row: the extracted brain which is at the finest level.
Figure 2.12: The red circle is the dark volume may be classified as background that need to be post-processed
Figure 2.13: (a) The red circle is the dark voxels in brain but be classified as background, and yellow circle is the voxels deleted in step1of flowchart but belong to brain region, that both need to be post-processed. (b) We find contour of the extracted brain from the highest resolution level. (c) Take the outermost contour. (d) Fill the voxels which are inside the outermost contour and the dark voxels in brain are filled.
set.
The first IBSR data set consist of 20 MR images, each with around 60 coronal slices, matrix size 256 mm × 256 mm, FOV 256 mm × 256 mm, and slice thickness 3.1 mm. The most images in this data set involve apparent intensity inhomogeneity and the large region of neck. The large region of neck may bias the estimation of parameters in the HWA and BET method. Therefore, the large region of neck were manually removed before estimation [25].
The second IBSR data set is composed of 18 MR images, each with around 128 coronal slices, matrix size 256 mm × 256 mm, FOV 240 mm × 240 mm, and slice thickness 1.5 mm. Each image has brain region determined by manual segmentation. The IBSR data sets as well as the brain region images are from the Montreal Neurological Institute and are available at http://www.cma.mgh.harvard.edu/ibsr.
The VGATPE data set contains 18 MR volumes of healthy subjects collected by Taipei Veterans General Hospital, each with around 124 coronal slices, matrix size 256 mm × 256 mm, FOV 260 mm × 260 mm, and slice thickness 1.5 mm. The ages of these subjects are between 22 and 57 years old. The brain region images of all subjects are segmented manually.
BrainWeb phantom images from the Montreal Neurological Institute are simulated by various levels of noise and intensity non-uniformity. The noise levels are 0%, 1%, 3%, 5%, 7%, and 9%. Each noise level contains three levels of intensity non-uniformity, which are
0%, 20%, and 40%. Each image with around 217 coronal slices, matrix size 181 mm × 181 mm, FOV 181 mm × 181 mm, and slice thickness 1mm.
3.2
Evaluation metrics for brain extraction
This work used five coefficients to compare the accuracy of the brain regions extracted from existing methods.
The Jaccard similarity coefficient (JSC) is to measure the similarity of ground truth A and the extracted brain B by the following equation:
J SC(A, B) = |A ∩ B|
|A ∪ B| , (3.1)
where |x| denotes the voxel number of x. The value of JSC is between 0 and 1 and the larger JSC value means the better overlap of A and B.
The higher sensitivity Se and higher specificity Sp means more accurate extraction of
brain region. Sensitivity Se is to measure the ratio of extracted brain region to total brain
region:
Se=
T P
T P + F N . (3.2)
If more brain voxels are obtained, the Se value will be higher. Specificity Sp is to measure
the ratio of rejected non-brain region to total non-brain region: Sp =
T N
T N + F P . (3.3)
If more non-brain voxels are rejected the Sp value will be higher. The true positive rate,
TP, is the number of voxels extracted as brain correctly. The number of voxels incorrectly extracted as brain is denoted by FP. The true negative rate, TN, is the number of voxels extracted as non-brain correctly. The number of voxels incorrectly extracted as non-brain region is denoted by FN.
Figure 3.1: The left column is the discrete anatomical model of WM and GM. The right column is the discrete anatomical model of WM , GM and CSF.
Two probabilities of miss rate, pmand pf, were defined in [8]. The probability of missed
detection of the brain voxels is defined as:
pm =
|A − B|
|A ∪ B| . (3.4)
The probability of miss rejection of the non-brain voxels is defined as:
pf =
|B − A|
|A ∪ B| . (3.5)
The lower these two miss rates means the better brain extraction.
BrainWeb provides discrete anatomical model to label tissues and the CSF region pro-vided by BrainWeb contains ventricle areas as well as peripheral CSF areas, which are generally considered as non-brain regions. Therefore, we evaluate the brain extraction al-gorithms with this data set according to two criteria. One only considers WM and GM, as shown in left image of Figure 3.1, the other considers WM, GM, and CSF as the total brain region as shown in the right image of Figure 3.1.
3.3
Experimental results
Tables 3.1, 3.2, 3.2, 3.3, 3.4, 3.5, and 3.6 list the performance evaluation results of BET, BSE, HWA, ISTRIP and the proposed method for different data sets. In addition to the Table 3.1, in other cases, although HWA has the highest Se, its poor performance of Sp
results in its worse performance of JSC.
Since some of the brain volumes in the first IBSR data set has apparent intensity inho-mogeneity and artifacts, as shown in the Figure 3.2. Therefore, not all brain volumes have the satisfied result. As [25], Table 3.1 excludes the extracted brains which the JSC value is below 0.6 or the result is blank. The number of excluded brains are four, three, and five for the BSE, HWA, and the proposed method, respectively. Performance evaluation using the first IBSR data set without including all the failure cases as shown in Table 3.2, and the number of excluded cases is seven. In Table 3.1 and Table 3.2 we has the similar JSC value with ISTRIP, In terms of Seand Sp, the proposed method has better sensitivity than
ISTRIP. However, ISTRIP can extract the all brain volumes in the first IBSR data set.
The proposed method outperforms others for the second IBSR set, in terms of JSC, Se,
and Sp. As listed in Table 3.3, compared to HWA, our method maintains the high sensitivity
of HWA as much as possible and improves the low specificity of HWA to the highest one. Our method achieves the highest JSC, the second highest Se, and the highest Sp. In this
data set, the total amount of voxels in an MR volume is 8388168, and the brain region is about one-eighth of the total voxels. Therefore, 0.01 voxels Se of the value means there
are about 10000 of the extracted brains overlapped with the ground truth. The 0.01 voxels Spof the value means there are about 70000 of the extracted non-brain regions overlapped
with the non-brain regions of the ground truth. The 0.01 voxels JSC of the value means there are more than 10000 of the extracted brains overlapped with the ground truth.
The performance evaluation using the BrainWeb phantom data set and only considering WM and GM is listed in Table 3.5. Although BSE has the highest Sp, but has the lowest
Se. The proposed method has the best performance compared with BET and ISTRIP, in
terms of JSC, Se, and Sp.
Since the ground truth of the CSF region provided by BrainWeb contains peripheral CSF area, the values of Sp when considering CSF, as listed in Table 3.6, higher than those
in Table 3.5. In Table 3.6, the higher Se means that the segmentation result may contain
more peripheral CSF voxels. Therefore, we only concern the Sp value with this data set,
and we have the similar Sp of BET which has the highest value.
The parameters of the second IBSR data set are as same as those used in [25]. For each brain extraction method, except BSE, the parameters of the VGHTPE and BrainWeb phantom data sets are the same. The fractional intensity threshold of BET was set as 0.7; the parameters of HWA were set to the default values with surface-shrink option turned on; the intensity contrast of ISTRIP was set as 0.08 with ”further remove non-brain tissues with high intensity” option turned on; the weight of intra n-link was set as 0.35, t-link was set as 0.65 and the distance term was set as 0.07 in the proposed method. The edge constant, diffusion iteration, and diffusion constant, erosion size of BSE when using the VGHTPE (BrainWeb phantom) data set were set to be 25 (increased with the degree of non-uniformity), 3 ,0.62 and 1 (2), respectively.
The extracted brain images using the proposed method with the subject in the first/second IBSR data set, VGHTPE, and BrainWeb phantom image are shown in Figure 3.3, Figure
Figure 3.2: The example of brain volumes with apparent intensity inhomogeneity and arti-facts in the first IBSR data set.
BSE4 0.900 (0.025) 0.954 (0.035) 0.993 (0.003) 0.044 (0.034) 0.055 (0.017)
ISTRIP 0.910 (0.018) 0.986 (0.013) 0.991 (0.005) 0.013 (0.013) 0.077 (0.027) HWA3 0.752 (0.037) 0.974 (0.068) 0.970 (0.008) 0.022 (0.056) 0.226 (0.026)
Our method5 0.910 (0.021) 0.989 (0.012) 0.991 (0.003) 0.010 (0.011) 0.079 (0.020) The superscripts in the first column indicate the excluded cases.
Table 3.2: Performance evaluation using the first IBSR data set without including all the failure cases.
Method JSC Se Sp Pm Pf BET7 0.881 (0.017) 0.981 (0.026) 0.988 (0.003) 0.017 (0.024) 0.102 (0.021) BSE7 0.905 (0.025) 0.954 (0.035) 0.994 (0.002) 0.044 (0.034) 0.051 (0.010) ISTRIP7 0.911 (0.014) 0.988 (0.015) 0.991 (0.004) 0.012 (0.014) 0.077 (0.023) HWA7 0.762 (0.012) 0.999 (0.001) 0.967 (0.008) 0.001 (0.001) 0.237 (0.014) Our method7 0.910 (0.021) 0.991 (0.008) 0.991 (0.003) 0.009 (0.007) 0.081 (0.020) The superscripts in the first column indicate the excluded cases.
3.4, Figure 3.5, and Figure 3.6, respectively. In our method we bring up three points, dis-tance term at the highest level, pre-determined foreground and application of graph cuts algorithm localized by using smaller cubes to local estimation of intensity histogram. The left column of Figure 3.7 illustrations the extracted brains without using the distance term,
Table 3.3: Performance evaluation using the second IBSR data set. Method JSC Se Sp Pm Pf BET 0.891 (0.052) 0.959 (0.042) 0.989 (0.005) 0.038 (0.038) 0.071 (0.031) BSE 0.838 (0.083) 0.957 (0.042) 0.973 (0.030) 0.041 (0.041) 0.119 (0.104) ISTRIP 0.915 (0.018) 0.978 (0.011) 0.990 (0.003) 0.021 (0.011) 0.064 (0.022) HWA 0.814 (0.040) 0.9997 (0.0003) 0.965 (0.016) 0.0002 (0.0002) 0.186 (0.040) Our method 0.930 (0.011) 0.981 (0.016) 0.991 (0.005) 0.018 (0.015) 0.052 (0.017)
Table 3.4: Performance evaluation using the VGHTPE data set.
Method JSC Se Sp Pm Pf BET 0.921 (0.009) 0.971 (0.013) 0.994 (0.001) 0.028 (0.012) 0.051 (0.011) BSE 0.933 (0.008) 0.958 (0.011) 0.997 (0.001) 0.041 (0.011) 0.026 (0.008) ISTRIP 0.915 (0.011) 0.984 (0.012) 0.991 (0.002) 0.015 (0.011) 0.07 (0.015) HWA 0.847 (0.019) 0.998 (0.002) 0.979 (0.003) 0.002 (0.002) 0.151 (0.019) Our method 0.920 (0.010) 0.982 (0.009) 0.992 (0.002) 0.017 (0.009) 0.063 (0.017)
and the right column is the extracted brain using distance term. The difference is obvi-ous, the border of the extracted brains in the right column are more smooth. Figure 3.8 shows the extracted brain of the proposed method overlaid by another without considering pre-determined foreground. The brain region of proposed method is larger than the one without using pre-determined foreground. Some brain extraction methods segment brain volume in a slice-by-slice manner, however we use smaller cubes to do this work. The left column of Figure 3.9 shows the extracted brain of the proposed method overlaid with a brain extracted by graph cuts algorithm in a slice-by-slice manner, and the right column
HWA 0.696 (0.005) 0.9995 (0.0001) 0.876 (0.003) 0.0003 (9.00947E-05) 0.304 (0.005) Our method 0.840 (0.025) 0.997 (0.001) 0.946 (0.011) 0.003 (0.001) 0.157 (0.026)
Table 3.6: Performance evaluation using the BrainWeb phantom data set and consid-ering WM, GM and CSF. Method JSC Se Sp Pm Pf BET 0.928 (0.009) 0.944 (0.011)) 0.993 (0.001) 0.055 (0.011) 0.017 (0.003) BSE 0.870 (0.031) 0.897 (0.009) 0.988 (0.018) 0.0999 (0.012) 0.030 (0.041) ISTRIP 0.934 (0.001) 0.963 (0.004) 0.988 (0.001) 0.036 (0.004) 0.031 (0.003) HWA 0.850 (0.003) 0.993 (0.002) 0.936 (0.002) 0.006 (0.002) 0.144 (0.004) Our method 0.912 (0.007) 0.934 (0.017) 0.991 (0.005) 0.064 (0.0178) 0.024 (0.013)
shows the extracted brain of proposed method. The brain region of proposed method is larger than the one using the slice-by-slice manner. Since the upper part of the image in the right column of Figure 3.9 is darker than the bottom part of image in transverse view, therefore, applying graph cuts algorithm localized by using the local estimated histogram is preferable than slice-by-slice manner.
The performance evaluations of the proposed method without distance term/predetermined term and in a slice-by-slice manner using the second IBSR data set are listed in Table 3.7. Although the proposed method without predetermined term has the highest JSC value and specificity, but it has lower sensitivity than the proposed method with the three factors.
According to the value of JSC, Se, and Sp, distance term, pre-determined foreground and
application of graph cuts algorithm localized by using smaller cubes to local estimated intensity histogram are important factors in our method.
Table 3.7: Performance evaluation of the proposed method without distance term/predetermined term and in a slice-by-slice manner using the second IBSR data set.
Method JSC Se Sp Pm Pf
Our method 0.930 (0.011) 0.981 (0.016) 0.991 (0.005) 0.018 (0.015) 0.052 (0.017) Our method without distance term 0.911 (0.023) 0.975 (0.017) 0.989 (0.007) 0.024 (0.016) 0.066 (0.026) Our method in slice-by-slice manner 0.904 (0.020) 0.953 (0.024) 0.991(0.006) 0.045 (0.023) 0.051 (0.022) Our method without predetermined term 0.933 (0.012) 0.979 (0.016) 0.992 (0.004) 0.020 (0.016) 0.047 (0.016)
Figure 3.3: The extracted brain images of subject in the first IBSR data set and using the proposed method are shown in (a) transverse, (b) coronal and (c) sagittal views.
Figure 3.4: The extracted brain images of subject in the second IBSR data set and using the proposed method are shown in (a) transverse, (b) coronal and (c) sagittal views.
Figure 3.5: The extracted brain images of subject in the VGHTPE data set and using the proposed method are shown in (a) transverse, (b) coronal and (c) sagittal views.
Figure 3.6: The extracted brain images of subject in the BrainWeb phantom data set and using the proposed method are shown in (a) transverse, (b) coronal and (c) sagittal views.
Figure 3.7: The result of extracted brain without using distance term at the highest level. The left column is the extracted brain without considering the contour factor. The right column is the extracted brain and consider the contour factor corresponding to the left column.
Figure 3.8: The extracted brain of proposed method overlaid by another without consider-ing pre-determined foreground factor.
Figure 3.9: The left column shows the extracted brain of proposed method overlaid with a brain extracted by graph cuts algorithm in a slice-by-slice manner, and the right column shows the extracted brain of proposed method.
Figure 4.1: The left column is the extracted brain with excluded cerebellum voxels by HWA and the right column is the corresponding ground truth volume segmented manually.
In this work we use an initial brain with high sensitivity and low specificity segmented by HWA and trim the brain into a more accurate one. If the initial brain extraction fails then our method fails too. In most cases HWA is a method with high sensitivity for brain extraction, but it sometimes excludes parts of the cerebellum, The example is shown in Figure. 4.1.
The program execution speed of the existing methods for brain extraction which com-pared in this work are fast (less than a minute). The program execution speed of the pro-posed method is about 5 minutes for an MR volume. Regarding to the performance eval-uation for brain extraction, our method outperforms others for the second IBSR data set. Therefore, the time cost of the proposed method is not expensive. The compared method
Figure 4.2: An example of the possible disadvantage by using localized graph cuts. The red window shows the inconsistent border.
in this work, BET, HWA, and ISTRIP are performed on Genuine Intel CPU 2.93GHz×16 processor running Linux, BSE and our method are performed on Intel Core2 Quad XP 2.83GHz processor on Windows system. Most of the running time of our method spent on the Otsu’s method, since it was performed frequently. Although the speed of multires-olutional graph cuts with narrow band should be faster, but we use multiresmultires-olutional graph cuts to extracted brain instead of using multiresolutional graph cuts with narrow band. The reason is that the multiresolutional graph cuts with narrow band will limit the refinement of the results in a band and may miss the thin structure.
In our method we perform graph cuts algorithm individually for each cube, the advan-tage is to consider the influence of the intensity inhomogeneity. In the finest level, when the cube contains only a small portion of brain region. The intensity of foreground seeds may be similar to the non-brain region. Therefore the results may contain non-brain and have irregular border, as shown in the red window of Figure. 4.2. The distance term can correct this problem, but if the distance term is too big then the result may not be smooth.
According to our experiences, the result of adding the WM voxels to the pre-determined foreground will be better than pre-determined foreground only. Because the pre-determined foreground is eroded and upsampled from the extracted brain at the coarser level, in the
Figure 4.3: An example of the possible disadvantage by using the distance term in the finest level. (a) The extracted brain at the highest level and the red windows show the non-brain voxles. (b) The yellow line shows the contour used in distance term at the highest level, and the positions corresponding to the red windows in (a) are inconsistent with the brain area. (c) The extracted brain at middle levels. (a) and (b) are shrunk as the size of (c).
process of erosion may lose some thin structure or other details. But considering the time cost of WM region growing, we did not apply the WM region as pre-determined foreground at middle levels.
The advantage of the distance term at the finest level is its capability to constrain the shape of brain boundary. However, it has disadvantages. In some cases the extracted brain at middle levels are better than that at the highest level. Since the contour used in distance term is upsampled from the coarser level, the contour may not be entirely consistent with the border of the brain at the highest level, as illustrated in Figure 4.3.
We have compared to the existing methods and compared with our method by using four data sets. Regarding the performance evaluation for brain extraction, our method outperforms others for the first/second IBSR data set and BrainWeb phantom data set, and comparably with the BET and ISTRIP methods for the VGHTPE data sets.
[1] John Ashburner and Karl J. Friston. Voxel-Based Morphometry-The Methods. NeuroImage, 11:805-821, 2000.
[2] Catriona D. Good, Ingrid S. Johnsrude, John Ashburner, John Ashburner, Karl J. Friston, and Richard S. J. Frackowiak. A voxel based morphometric study of ageing in 465 normal adult human brains. NeuroImage, 14:21-36, 2001.
[3] M. Kubicki, M. E. Shenton, D. F. Salisbury, Y. Hirayasu, K. Kasai, R. Kikinis, F. A. Jolesz, and R. W. McCarley. Voxel-Based Morphometric Analysis of Gray Matter in First Episode Schizophrenia. NeuroImage, 17:1711-1719, 2002.
[4] Andrew M. McIntosh, Dominic E. Job, T. William J. Moorhead, Lesley K. Harri-son, Karen Forrester, Stephen M. Lawrie, and Eve C. Johnstone. Voxel-based mor-phometry of patients with schizophrenia or bipolar disorder and their unaffected relatives. Biological Psychiatry, 56:544-552, 2004.
[5] Paul M. Thompson, Agatha D. Lee, Rebecca A. Dutton, Jennifer A. Geaga, Ki-ralee M. Hayashi, Mark A. Eckert, Ursula Bellugi, Albert M. Galaburda, Julie R. Korenberg, Debra L. Mills, Arthur W. Toga, and Allan L. Reiss. Abnormal Cortical Complexity and Thickness Profiles Mapped in Williams Syndrome. Neuroscience, 25(16):4146-4158, 2005.
tenberg, and Richard M. Leahy. Magnetic Resonance Image Tissue Classification Using a Partial Volume Model. NeuroImage, 13:856-876, 2001.
[9] F. S´egonne, E. Busa, M. Glessner, D. Salat,A.M. Dale, H.K. Hahn, and B. Fischl. A hybrid approach to the skull stripping problem in MRI. NeuroImage, 22:1060-1075, 2004.
[10] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. Introduction To Algorithms, Second Edition. 2001.
[11] Yuri Y. Boykov and Marie-Pierre Jolly. Interactive Graph Cuts for Optimal Bound-ary & Region Segmentation of Objects in N-D Images. Internation Conference on Computer Vision,vol.I,p.105, 2001.
[12] Brian L. Price, Bryan Morse, Scott Cohen. Geodesic Graph Cut for Interactive Image Segmentation. Conference on Computer Vision and Pattern Recognition, 2010.
[13] Yin Li, Jian Sun, Chi-Keung Tang, Heung-Yeung Shum. Lazzy Snapping. ACM, 2004.
[14] Carsten Rother, Vladimir Kolmogorov and Andrew Blake. ”GrabCut”-Interactive Foreground Extraction using Iterated Graph Cuts. ACM Transactions on Graphics, 23:309-314, 2004.
[15] Jean Stawiaski, Etienne Decenci`ere and Franc¸ois Bidault. Interactive Liver Tumor Segmentation Using. MICCAI, 2008.
[16] Yin Li, Jian Sun, Heung-Yeung Shum. Video Object Cut and Paste. ACM Transac-tions on Graphics, 24:595-600, 2005.
[17] J. Reese and W. Barrett. Image Editing with Intelligent Paint. EUROGRAPHICS, 2002.
[18] Xue Bai, Jue Wang, David Simons, Guillermo Sapiro. Video SnapCut: Robust Video Object Cutout Using Localized Classifiers. ACM SIGGRAPH, 2009.
[19] Herve Lombaert, Yiyong Sun, Leo Grady, Chenyang Xu. A Multilevel Banded Graph Cuts Method for Fast Image Segmentation. ICCV, 259-265, 2005.
[20] Ali Kemal Sinop and Leo Grady. Accurate Banded Graph Cut Segmentation of Thin Structures Using Laplacian Pyramids. 896-903MICCAI, 2006.
[21] Suresh A. Sadananthan, Weili Zheng, Michael W.L. Chee, Vitali Zagorodnov. Skull stripping using graph cuts. NeuroImage, 49:225-239, 2010.
[22] Yuri Boykov and Vladimir Kolmogorov. An Experimental Comparison of Min-Cut/Max-Flow Algorithms for Energy Minimization in Vision. IEEE Transactions on Pattern Analysis and Machine Intelligence, 2004.
[23] David M. Mount and Sunil Arya. ANN: A Library for Approximate Nearest Neigh-bor Searching, Internet: www.iceengg.edu/staff.html, Jan 27, 2010 [Sep. 17, 2010]. [24] Ostu N. Threshold selection method from gray-level histograms. IEEE Trans Syst
Man Cybern, 9:62-66, 1979.
[25] Jia-Xiu Liu, Yong-Sheng Chen, Li-Fen Chen. Accurate and robust extraction of brain regions using a deformable model based on radial basis functions. Joural of Neuroscience Methods, 183:155-166, 2009.
