Research Article  Open Access
Joint Brain Parametric Map Segmentation and RF Inhomogeneity Calibration
Abstract
We propose a constrained version of Mumford and Shah's (1989) segmentation model with an informationtheoretic point of view in order to devise a systematic procedure to segment brain magnetic resonance imaging (MRI) data for parametric Map and weighted images, in both 2D and 3D settings. Incorporation of a tuning weight in particular adds a probabilistic flavor to our segmentation method, and makes the 3tissue segmentation possible. Moreover, we proposed a novel method to jointly segment the Map and calibrate RF Inhomogeneity (JSRIC). This method assumes the average value of white matter is the same across transverse slices in the central brain region, and JSRIC is able to rectify the flip angles to generate calibrated Maps. In order to generate an accurate Map, the determination of optimal flipangles and the registration of flipangle images are examined. Our JSRIC method is validated on two human subjects in the 2D Map modality and our segmentation method is validated by two public databases, BrainWeb and IBSR, of weighted modality in the 3D setting.
1. Introduction
Brain structure segmentation is the apportionment of brain tissue into gray matter and white matter, based on the appearance of tissue in images produced by magnetic resonance imaging (MRI). Because manual tracing of the boundaries between tissues in the brain is labor intensive, difficult, errorprone, and unrealistic for large amounts of data, an automated or semiautomated segmentation technique is needed for either visualization or diagnosis. Different imaging modalities, such as weighted, weighted, or Proton Density (PD) images, have been used for different segmentation methods. weighted images, because of their good contrast [1], have been widely tested for various segmentation methods [2–6]. A Map is a parametric image of pure (spin lattice relaxation time), derived from the solution of an equation describing tissue relaxation, and a parametric Map which is different from a weighted image. The relationship between and several diseases, such as schizophrenia or sickle cell disease, has been studied [7–10], and may be used as a possible indicator of pathology. Change in of certain voxels in the brain over time may be an early indicator of possible pathology [8]. Therefore, the segmentation of a parametric brain Map may highlight pathology unseen by other imaging approaches.
Past research has studied the segmentation of cortex in the brain. The three tissues white matter (WM), gray matter (GM), and cerebrospinal fluid (CSF) constitute the main parts in the brain. The goal is to find their respective boundaries. Different methods have been proposed to achieve this goal, and they may be classified into various categories: fuzzy segmentation methods [4, 11], Markov random field (MRF) methods [12, 13], Bayesian methods [14], active contour methods [1, 5, 6], or the combinations of two or more techniques. Some of these combinations are as follows: Leemput et al. [3] used expectation maximization (EM) and MRF, Xu et al. [2] used fuzzy segmentation and deformable surfaces, and Zhang et al. [15] combined a hidden Markov random field and an EM algorithm. Our method falls in the active contour category in the regionbased formulation. It is an adaptive version of MumfordShah's model [16] to systematically segment different tissues in the brain.
Before the segmentation, several issues arise for obtaining a Map. A Map may be calculated by a rapid method known as the variable nutation method [17] which provides comparable precision but much faster speed over conventional methods [18]. This method requires the acquisition of a set of flip angle images and the information can be extracted therein. The problem of determining the set of optimal flip angles therefore was studied. Deoni et al. [18, 19] proposed methods to determine the optimal flip angles by basically maximizing the signaltonoise ratio (SNR) of a Map. Their first work [18] required the knowledge of average TR/ in advance though, and their second work [19] introduced a weighted least square method to estimate the angles. We take another approach to determine the flip angles to achieve the tradeoff between acquisition time and accuracy.
Since a Map is generated by a set of flip angle images, alignment of the images is important in order to obtain an accurate Map. We therefore propose a method to register the raw data. The registration of flip angle images is usually ignored though because the movement of the head in the coil is minute. However, slight registration errors affect the resulting Map dramatically, as will be shown.
Radiofrequency (RF)inhomogeneity is another unavoidable problem encountered in MR imaging. The nonuniform distribution of the RF field can cause the resulting images to have low contrast and inhomogeneous intensity, which makes quantitative description and segmentation of the image difficult [20]. RF inhomogeneity affects the generation of Map in the sense that the spins are not tilted by the predefined nominal angles [17, 20–22]. Hence, we focus on calibrating the RF nonuniformity in order to generate an accurate Map. Cheng and wright [17] calculated an analytical form of errors induced by RF nonuniformity and allowed simple correction of measurements. Both Wang et al. [20] and Venkatesan et al. [22] incorporate a scaling factor to rectify the nominal flip angles in their models. We propose a method, which assumes the average value among a transverse plane is the same across the central brain slices, to jointly segment a Map and calibrate the RF nonuniformity along the direction perpendicular to the transverse plane. This assumption may seem disputable because indeed exhibits regional heterogeneity in human cortex [23]. Taking into account that we average out the heterogeneity of , and doing so only within few central brain slices, this assumption help calibrate RFinhomogeneity and generate a quantitative Map for segmentation.
The paper is organized as follows. In Section 2 we state our adapted model, with a fast implementation method. In Section 3 we first illustrate some necessary preprocessing to generate an accurate Map, which includes the registration of flip angle images with the generation of a brain mask and the determination of optimal flip angles, to achieve a tradeoff between quality and efficiency. We then describe a systematic procedure to segment a brain Map. In Section 4 we propose a novel method to jointly segment a Map and calibrate RFinhomogeneity, and in Section 5 we show the results of registering the flip angle images, determining optimal flip angles, and segmenting the Maps. We also show validations, the resulting Map after RFcalibration, and the 3D segmentation for two sets of weighted data. At last in Section 6 we discuss our findings and offer some concluding remarks in Section 7.
2. Proposed Model for Segmentation
Active contour methods comprise a popular segmentation technique in which an initialized contour is driven by a partial differential equation (PDE) to minimize an energy functional designed to attract the contour toward image edges. Active contour methods can be classified into two categories: edgebased [24–27, 27–34] and regionbased [16, 35–44] methods. Edgebased methods examine the gradient information of the image, and stop the contour whenever the gradient is high. However, there are many situations when the edge is not clearly characterized by the gradient, and it has been shown that regionbased methods outperform edgebased methods [35–37, 41, 43, 44]. By examining some statistics of the region inside and outside of the active contour, and optimizing the separation of these two statistics, we may achieve a better segmentation performance, thus making regionbased methods more attractive.
Our proposed model is a modified Momford and Shah [16] functional, and falls in the regionbased category. Two adaptations of MumfordShah's model constitute the novelty of our proposed technique. First is the incorporation of an informationtheoretic view, which characterizes the statistical property of data. The second is a selective weighting, similar to the weight parameter in [35], which favors erring towards one tissue type or another, and thus make 3tissue segmentation possible.
2.1. Adapted MumfordShah Functional
We use a modified MumfordShah [16] energy functional: where approximate a function applied of the image for region , or . is smooth within each region , but not across the boundaries; denotes the region boundaries, which is the contour of interest in this paper, and are weights such that . By minimizing this functional we obtain a function which is faithful to the image (first term) and smooth within each region but not across the boundaries (second term), while penalizing excessive length of the boundaries (last term).
The first adaptation is the function applied on to image instead of the image itself. It is motivated by the work of Unal et al. [45]. Specifically, an informationtheoretic approach for maximizing a probabilistic disparity measure, JensenShannon (JS) divergence, was proposed. A constructed function characterizes the property of the probability density function (PDF) of the image intensity such as skewness (), or kurtosis (), relative to a Gaussian [46]. A proper choice of will capture the statistical characteristics of the data and will hence achieve a good segmentation.
The second adaptation, similar to that in [35], is the selective weight applied on the inside/outside energy terms. This provides a probabilistic assignment to the segmented regions, and also makes 3tissue segmentation possible, as will be shown in Section 3. Enhancing the weight is tantamount to penalizing both the error of the difference between the approximated function and the data fidelity term , as well as the degree of smoothness of . This would yield a smaller segmented region which is likely to be more faithful to the image, and of “purer” tissues. Figure 1 shows examples of white matter (WM) segmentation with different weights. With a larger insideweight , the segmented region is smaller and the segmented tissues are purer.
(a)
(b)
(c)
2.2. Fast MumfordShah Implementation
Minimizing the MumfordShah energy functional involves solving for the approximating functions and for the contour . The joint search for these infinite dimensional unknowns usually entails gradient descent flows. In particular, the approximated functions are typically modeled as a linear combination of a basis set whose dimension equals the number of pixels in the image, that is, each pixel is assumed to be independent [47]. The curve and the approximated functions are then evolved iteratively along the gradient flows. This implementation is computationally daunting and therefore a faster implementation method was in order. Alvino and Yezzi [48] proposed a fast implementation using a significantly smaller basis number to model the approximated functions, while still achieving sufficient resemblance to the obtained functions when using the pixelbypixel basis.
We adopt their socalled linear heat equation basis with the change from to to incorporate the statistical properties of the data from an informationtheoretic point of view in the previous section. We hence have, where is the average function, and the coefficients , , may be similarly derived in [48]. Our derivations are included in the appendix.
Substituting (2) into (1), and using a classical methods of variational calculus, the gradient descent evolution of the curve may be derived as where denotes the curvature of the contour , an artificial time evolution parameter, and the outward normal of the contour. are derived by substituting from (2) into (3) as shown in the appendix.
3. Segmentation of a Parametric Brain T_{1}Map
A brain parametric Map is the calculated result (variable nutation method [17]) from several flip angle images. Because the flip angle images are acquired at different times, and because the subject may move during image acquisition, registration should be carried out first to obtain an accurate Map. Moreover, in order to achieve a balance between the acquisition time and the resulting Map quality, we propose a method to determine optimal flip angles. In the following subsections we will first illustrate how we register the flip angle images and obtain a brain mask as a byproduct, then describe our method to determine a set of optimal flip angles, and at last describe our proposed procedure to segment a Map.
3.1. Registration of FlipAngle Images and Generation of Brain Mask
The value of is traditionally determined by acquisition methods such as Inversion Recovery (IR) or SaturationRecovery (SR). Other rapid methods, such as variable nutation (the DESPOT method) have been proposed [17], and require acquisition of several flip angle images, and calculation of . Since these flip angle images are acquired at different times, registration must first be accomplished. Even though the interval between consecutive scans may be as short as two minutes, and movement of the subject's head inside the receiver coil may be a few pixels (under the resolution of a image), the effects of such offregistration can be significant, as shown (Figure 5). Here we describe a method to register flip angle images and jointly obtain a mask, as a byproduct, to get rid of the skull and other structures around the brain.
We use a joint segmentation and registration (JSR) technique proposed by Yezzi et al. [49] with an additional tuning weight, to achieve registration and to obtain the mask. The theory of JSR technique consists of evolving two contours, with a enforced relationship (ex: rigid or affine transform) between them, in two images according to a partial differential equation (PDE) which is a result of optimizing, for example, the sum of two energy functionals.
For our particular task, we may choose a regionbased energy, such as Chan and Vese's model [35] incorporating weights, which arises as a special case of (1) with the data term and the approximating function inside and outside the contour. This model approximates the image by a simple piecewise constant function, which suits our goal here because it creates a mask that divides the image into two parts brain region and nonbrain region.
We observe that the boundary between the brain and nonbrain is visually more easily distinguished than the boundary between different tissues within the brain, for every flip angle image. We therefore put a very small weight on the inside energy in (1) to penalize very little the difference between the data inside the contour, , and the approximated function . Experimental results show that yields satisfactory results.
We carry out this registration technique pairwise for all flip angle images, and once the flipangel images are registered to each other, they are used to generate a Map.
3.2. Determination of Optimal FlipAngles
The determination of the flip angles that yields qualitative Map and saves time, is mainly by comparing the difference between the gold standard and the Maps generated from different combinatorial flip angles. We begin by acquiring images at a set of 19 flip angles which spans the range of standard angles [50] and will give an optimized Map. This is also confirmed by observing that the generated Map gives the best quality. We call this particular Map the gold standard, or , and denote this set of angles as . We then compare the Map generated from all combinations of the subset of with . Out of these combinations we select the optimal subset of angles. Since the data acquisition is slicebased, this study of optimal flip angle focuses on the central slice, which is least affected by RFinhomogeneity [51] and is routinely selected at the level of the basal ganglia, including both the genu and the splenium of the corpus callosum, and generally shows the putamen and lateral ventricle [52]. We denote the optimal angles , which is a subset of , as those that exhibit the smallest difference between and the Map generated by flip angles: where , , denotes the Map generated by , and the summation of is over the whole brain region at the central slice.
Equation (5) therefore gives 18 sets, with the number of elements ranging from 2 to 19, of optimal angles. For the determination of the optimal set of flip angles, in reaching a compromise between efficiency and Map quality, we examine the error between and and compute the error rate. The error rate is defined as where is the area of the brain, and is the error, defined as 1 if the difference between and is greater than some threshold , and 0 otherwise. The threshold is defined as the minimum of the standard deviation among the values of three brain tissues (WM, GM, and CSF) manually segmented by an expert. The summation of is over the whole brain at the central slice.
The plot of the error rate versus the number of flip angles is shown in Figure 6, where the inflection point of the fitted curve is at 10 flip angles (). We hereafter routinely use these 10 flip angles to generate the Map.
3.3. Brain T_{1}Map Segmentation Procedure
Once the brain mask is obtained (Section 3.1), it is used to segment away the skull, leaving only three major tissues in the image: WM, GM, and CSF. Notice that the curve evolution corresponding to the energy introduced in Section 2 always results in a “binary segmentation”, where we will have regions inside (foreground) and outside (background) the contour(s). We cannot simultaneously segment the three tissues, even though we will later talk about multiphase segmentation. We may, however, tune the weights (Section 2.1), to penalize the error between the data term and the approximated function (1) to segment one tissue at a time, with a progression analogous to that of “peeling an onion”.
We illustrate our Map segmentation procedure for a hard segmentation of three tissues, and the probabilistic segmentation is obtained by varying the weights around the value of the trained weight (Section 5.4). This procedure is applicable to both 2D and 3D data sets.
A Map segmentation procedure consists of three steps, where the first two steps are to evolve the contours by minimizing the energy in (1) with different weights, and the third step is just a simple subtraction. The procedure is as follows:
(1)Treat WM as the foreground, everything else as the background; let in (1) be the trained , and use it to segment WM in the interior region of the contour.(2) Treat WM and GM as the foreground, CSF and everything else as the background; let in (1) be the trained , and use it to obtain CSF in the exterior region of the curve filtered by the brain mask.(3) GM is obtained by subtracting WM and CSF from the whole brain.The procedure is based on the anatomical observation that GM is enclosed by CSF, and that CSF is separated from WM [1, 3], such that we may peel off one layer at a time.
The choice of function in (1), which is chosen to better capture the statistical property of Map (and for other modality), will be shown in Section 5. The values of and are determined through a training process. Suppose an expert's manual segmentation is regarded as the ground truth. If denotes the segmentation region by the expert, and denotes the segmentation region by the weight with some fixed , for some tissue, then the value of is determined by minimizing where denotes the number of pixels in region . The first term inside the bracket is an overlap metric (OM) [5] which will be discussed again in Section 4. It is commonly used in validating segmentation performance—the closer it is to 1, the better it has performed.
4. Joint T_{1}Map Segmentation and RFInhomogeneity Calibration
RFinhomogeneity is an unavoidable problem in MR imaging. The strength of the RF field varies within the MR scanner, such that the resulting image may be of low contrast or of nonuniform intensity. In what follows, we will first describe how a Map is calculated from a set of flip angle images, then how the Map is affected by RFinhomogeneity and then show our proposed correction method.
4.1. Variable Nutation Method with RFInhomogeneity
A Map is calculated by variable nutation [17]. Flipangle images are acquired using a FLASH sequence at different flip angles, and is calculated from the slope of a least square fitted line to the pair of data , where is the signal strength of the flip angle image expressed as a function of the flip angle .
RFinhomogeneity affects the computation of in that the spins are not tilted by the nominal angle because the strength of the RF field is not as predefined. This phenomenon appears more serious at the periphery of the receiver coil. Calibration of RFinhomogeneity, or, the correction of flip angles, produces a more accurate T_{1}Map. A spatially dependent scaling factor , therefore, is introduced to adjust the tilted angle [20, 22], that is, replacing by and is extracted from the slope of the line fitted to the above space. Our imaging setting is slicebased across different transverse planes (Figure 2), thus our proposed joint segmentation and RFinhomogeneity calibration (JSRIC) method is as an initial step to easily calibrate RFinhomogeneity vertically.
4.2. Flip Angle Rectification and Segmentation of T_{1}Maps
Our method is based on the assumption that the average of WM should be roughly the same across central brain transverse slices. exhibits regional heterogeneity in human cortex [23]. We however average out the heterogeneity of , and carry out the method only within few central brain slices. Our method therefore yields a better Map for segmentation and jointly rectifies flip angles. The advantage is that as the segmentation delineates more precisely the WM region, flip angle correction is therefore enhanced, and as flip angle correction is carried out more precisely, which gives better quality Map, segmentation is facilitated as well.
Our joint segmentation and RFinhomogeneity calibration (JSRIC) method requires first taking the average value for segmented WM at the central slice (which has relatively uniform B1 field [51, 53]) as the reference, and then iteratively segmenting and searching for the scaling factor in (8) for all other slices. This is a threestep iterative process (as shown inside the dashed box in flowchart 3).
(1)Segment WM for the central slice by the method proposed in Section 3.3, compute the average value of WM, and denote the average as (2)For slice , varies in (8), and calculate the corresponding Map, denoted as . If goes beyond , claim the current slice as uncalibratable and repeat this step for the next slice (3)Segment WM for , compute the average of WM, and denote it by . If , go back to step 2 and increase , otherwise claim it is done for the current slice. (4)Iterate steps 2 and 3 through all the slices.The variables , , and in flowchart 3q are all predefined parameters.
5. Experimental Results
In this section we show a series of results from segmenting a parametric brain Map, which includes the registration of flip angle images, generation of a brain mask as a byproduct, determination of the optimal flip angles, and joint Map segmentation and RFinhomogeneity calibration. In addition, we will apply our proposed segmentation method to another modality, weighted images, in a 3D setting, to show the generality of our segmentation method. Results for both modalities (Map and weighted) will be validated.
5.1. Subjects and Scan Protocol
The subjects being scanned were recruited from a clinic in the Department of Psychiatry at the University of North Carolina, under an IRBapproved protocol to image the brain. Informed consents were signed by subjects or their guardians. A transverse 3D FLASH sequence using different flip angles was acquired in a 3T Siemens MRI scanner with a quadrature head coil. The scan parameters were: TR = 25 msec, TE = 2.83 msec, 16 slices, and slice thickness = 5 mm. The center slice used for optimal flip angle study was routinely selected as described (Section 3.2).
5.2. Registration of FlipAngle Image and Brain Mask
The registration of a set of flip angle images is carried out pair wise, and Figure 4 shows two flip angle images and the resulting brain mask. This mask has otherwise to be generated manually or by other techniques [54]. Note that even though the flip angle images are offregistered by no more than four pixels in both the  and directions, the impact is obvious. Figure 5 shows two Maps generated by registered and unregistered flip angle images, and we can see that after registration the anatomical structure of the Map is more readily distinguished, and the high signal intensity artifact in the upper left of the unregistered Map map is reduced.
5.3. Determination of Optimal FlipAngles
The plot of the error rate (defined in (6) versus the number of flip angle is shown in Figure 6, and Section 3.2 already concluded that a set of 10 flip angles () is a good compromise between quality and efficiency. Figure 7 shows a comparison of generated Maps when using 2 optimized flip angles ([]), 6 optimized flip angles ([]) and 10 optimized flip angles, with an “error map” (error defined in Section 3.2) for each image. The Map generated by 2 angles exhibits substantially more error than the maps using 6 or 10 angles, as shown by the bright pixels in the error map. We conclude that 10 flip angles is an acceptable choice considering the tradeoff between accuracy and scan time.
5.4. FlipAngle Rectification and Segmentation of T_{1}Maps
We carried out JSRIC introduced in Section 4 on the Maps of two subjects. Out of a total of 16 slices, the bottom 4 slices were discarded because they did not cover sufficient WM for segmentation by our JSRIC method. The method follows the flowchart in the dashed box shown in Figure 3 from slice 5 to slice 15. Slice 16 could not be calculated, possibly because of its proximity to the boundary of the RF field and hence degraded.
(a)
(b)
(a)
(b)
(a)
(b)
(c)
The function introduced in (1) was empirically chosen as the cubic function , which characterizes the skewness of a PDF. Moreover, and are obtained by training according to (7) based on an expert's manual segmentation of one subject (training subject). These values are and . The same values are then applied to the other subject (testing subject). The segmentation of a Map is achieved by evolving contours according to (4). Since the curve evolution is based on a gradient flow, the final result may vary depending upon the initialization. A good initialization may avoid an undesirable local minimum. A Map does not necessarily have good contrast compared to a weighted image. A weighted image histogram is shown (Figure 8(a)), and a spectral analysis can be carried out to threshold the image as an initialization [5, 55], to achieve a good segmentation after fine tuning of the contour. The Map does not have the same level of contrast as a weighted image (Figure 8(b)). We therefore initialize the contours by either placing uniformly spaced squares or by manually seeding (by mouse clicking and dragging on the image). Both methods give similar performance with manual seeding slightly better. We will therefore show the results with manual seeding for Map.
(a)
(b)
Figure 9 shows the “map” ( versus slice number) for two subjects, where the parameter is as defined (Section 4.1). The strength of the RF field drops significantly at the top slices (slice 14 and 15) such that the corresponding flip angles have to be rectified by a scaling factor much smaller than 1. Figure 10 shows two Maps generated with and without flip angle rectification for the same subject at slice number 5, 6, 14 and 15. It is clearly seen that the top 2 slices of the Maps with calibration have much better quality than those without.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
Figure 11 shows some selective segmentation results of WM and GM for the test subject. Validation of these results requires comparison with manually segmented images. The commonly examined metrics which determine the performance of a segmentation are TP (true positive), FP (false positive), and OM (overlap metric) [1, 5]. The overlap metric is defined for a class assignment as the sum of the number of pixels that both have the class assignment in each segmentation divided by the sum of pixels where either segmentation has the class assignment. This metric approaches 1 for segmentations that are very similar, and is near 0 when they share no similarly classified pixels. The OM metric is usually compared for different segmentation methods. Figure 12 shows three overlap metric curves (OM versus slice number) of WM and GM segmentation for the training subject (the test subject exhibits a similar result). The first curve corresponds to the segmentation of calibrated Maps, with our tuned weights and cubic function , the second corresponds to the Maps generated by nonrectified flip angles (using the same segmentation parameters), and the third corresponds to the calibrated Map, with function and tuned weights ( and ). These results show that calibrated Maps enhance the segmentation performance compared to uncalibrated Maps, especially at the top two slices (slice 14 and 15) which are most seriously affected by RFinhomogeneity. The comparison of different functions shows that the performance of WM segmentation is comparable for the two functions. There is, however, a significant difference for CSF segmentation, thus affecting GM segmentation. The cubic function outperforms tremendously for GM segmentation. It also demonstrates that the cubic function better characterizes the statistical properties, the skewness, of the data. The OM for the other subject also has similar results.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
(a)
(b)
Figure 13 shows TP and FP for WM segmentation with different weights around the value of , to demonstrate the notion of our probabilistic segmentation. TP = , and FP = (, where and are the expert segmented regions and ours using respectively, and denotes the area of region . When the weight increases, so does the penalty for the error between the data term and the approximating function (1). Therefore TP and FP decrease correspondingly, and vice versa.
(a)
(b)
5.5. Segmentation of Weighted Images
In this section we show the generality of our proposed segmentation method in a 3D setting by applying it to a commonly exploited modality—weighted images. The same procedures are carried out as in Section 3.3, except that now the images are collated into volumes and the active contour is replaced by an active surface. We test it on two open databases accessible online— BrainWeb (http://www.bic.mni.mcgill.ca/brainweb/) [56] and IBSR (http://www.cma.mgh.harvard.edu/ibsr/). The former is a simulated brain MRI database, therefore ground truth is provided and the latter is genuine brain MRI data which also has been manually segmented by experts. We preprocessed the data to filter out everything except for the three main tissues, WM, GM, and CSF.
We first tested our segmentation method on two BrainWeb subjects (a training and a testing subject) using simulated brain T_{1}weighted data. The images are 1 mm slice thick, with 3% noise level, and 20% RF intensity nonuniformity (INU) and each data set was pixels. Because the contrast between different tissues was high, we did a histogram analysis and applied the threshold method similar to [5] for initialization. The function is still chosen as , and and are obtained by training as 0.3 and 0.2 respectively. The validation metrics for the testing subject are shown in Table 1. The results show that it achieves a high performance of OM being around 0.8 for three tissues. The computational time for one subject is less then 1 minute on a laptop with a 1.73 GHz CPU and 1 GB memory. Figure 14 shows the segmentation results for the testing brain.

(a)
(b)
We then tested our segmentation method on 20 data sets of weighted images provided by the Center for Morphometric Analysis at Massachusetts General Hospital on the IBSR website. The data set for each subject was , where ranges from 58 to 63 pixels. We arbitrarily chose one subject (subject 1_24) as the training data, empirically chose the function for WM and for CSF, , , and . We did not carry out any preprocessing to denoise the data or to decrease the RFinhomogeneity effect which seriously degraded the data. Therefore the spectral analysis could not be carried out and we used uniformly spaced cubes as the initialization. Out of 20 subjects we however still have 3 particularly unsuccessful cases (specifically subject 2_4, 16_3, and 111_2) due to RFinhomogeneity which we did not rectify. Including these three subjects, we have an averaged overlap metric around 0.651, and if excluding these 3 outliers, we achieved an averaged OM of around 0.707. The computational time is less then 5 minutes on the same PC as above. Figure 15 shows the OM for WM and GM segmentations compared to other techniques. The statistics show that our method outperforms most other methods even if the outliers are included; if excluding those, we have WM segmentation slightly poorer but GM segmentation slightly better than the current best segmentation method. Figure 16 shows the 3D segmentation for one subject.
(a)
(b)
6. Discussion
6.1. Registration of FlipAngle Images and Generation of Brain Mask
To register flip angle images we used the JSR technique, which has the advantage of jointly segmenting, registering, and generating a brain mask. Other techniques such as the informationtheoretic method [57] only register the images, and the brain mask has to be otherwise generated. Our JSR has been carried out pair wise on 10 flip angle images, registering 9 images to a reference image. In theory, however, it should be possible to integrate multipleimage in JSR's formulation. By summing the energy functionals of multiple images with a relationship enforced between each contour, that is, , where , , the evolution of the contours and registration parameters may be derived similarly.
6.2. Segmentation of Brain Data and Probabilistic Segmentation
Our segmentation of three brain tissues is based on the tuning of weights, to penalize differently the error of the approximated function, to obtain different regions of tissue. Segmentation uses the anatomical nature of brain tissue which has a layered structure such that we may peel off one layer at a time. Our method uses an active contour that is able to separate an image into two parts: the inside and the outside of the contour. However, multiphase active contour techniques exist [36, 58] and are able to evolve multiple contours simultaneously and segment multiple regions at once. These methods should theoretically be able to segment three brain tissues simultaneously. Our method is different in the sense that it has a probabilistic flavor that the tuning weights determine a puritylevel of segmented tissues.
6.3. Joint Segmentation and RFInhomogeneity Calibration
Our JSRIC method works by enforcing the average of white matter value to be homogeneous across different transverse planes in the central brain region, to find the scaling factors which affects the flip angles. Even though has regional heterogeneity, by taking the average of WM for the transverse plane we average out this heterogeneity, and doing so only in the central brain region. WM segmentation and RFcalibration enhance the precision of one another, as shown in the performance plot in the previous section.
7. Conclusion
In conclusion, we propose an adapted MumfordShah type energy functional for segmentation. The two adaptations are: (1) a function is able to characterize the statistical properties of the data to achieve better segmentation results, and (2) the tuning weights are able to segment brain tissues in a probabilistic fashion and achieve threetissue segmentation. We also propose a novel method (JSRIC) to jointly segment a Map and calibrate RFinhomogeneity. The whole procedure moreover includes the determination of optimal flip angles in achieving the balance between accuracy and efficiency, and joint registration of flip angle images and generation of brain mask. After RFcalibration, the top and bottom slices of Maps show better contrast and enhance the segmentation performance. The segmentation method has also been applied to weighted data, to show the generality of our segmentation method, and the results are validated by ground truth and by expert manual segmentations. The results show high performance with our method.
Appendix
Derivation of Gradient Flows for Our Proposed Energy
To derive the gradient flow of our energy in (1) with Alvino et al.'s [48] fast MumfordShah implementation, we first restate their result for the approximated function in scale space. Suppose the approximated function , or , is expressed as a linear combination of basis : then the coefficients can be solved by a linear system where and each matrices must be computed for each region. Also, are each vectors.
Therefore we may replace the basis by and as in (2), and solve for the coefficients : where
Then substituting the coefficients in (2) by the above results and then pluging into (3), we have in (4).
References
 X. Zeng, L. H. Staib, R. T. Schultz, and J. S. Duncan, “Segmentation and measurement of the cortex from 3D MR images using coupledsurfaces propagation,” IEEE Transactions on Medical Imaging, vol. 18, no. 10, pp. 927–937, 1999. View at: Google Scholar
 C. Xu, D. L. Pham, M. E. Rettmann, D. N. Yu, and J. L. Prince, “Reconstruction of the human cerebral cortex from magnetic resonance images,” IEEE Transactions on Medical Imaging, vol. 18, no. 6, pp. 467–480, 1999. View at: Publisher Site  Google Scholar
 K. Van Leemput, F. Maes, D. Vandermeulen, and P. Suetens, “Automated modelbased tissue classification of MR images of the brain,” IEEE Transactions on Medical Imaging, vol. 18, no. 10, pp. 897–908, 1999. View at: Google Scholar
 M. N. Ahmed, S. M. Yamany, N. Mohamed, A. A. Farag, and T. Moriarty, “A modified fuzzy Cmeans algorithm for bias field estimation and segmentation of MRI data,” IEEE Transactions on Medical Imaging, vol. 21, no. 3, pp. 193–199, 2002. View at: Publisher Site  Google Scholar
 H. Li, A. Yezzi, and L. D. Cohen, “Fast 3D brain segmentation using dualfront active contours with optional userinteraction,” in Proceedings of the International Workshop on Computer Vision for Biomedical Image Applications, vol. 3765 of Lecture Notes in Computer Science, pp. 335–345, 2005. View at: Google Scholar
 R. Goldenberg, R. Kimmel, E. Rivlin, and M. Rudzky, “Cortex segmentation: a fast variational geometric approach,” IEEE Transactions on Medical Imaging, vol. 21, no. 12, pp. 1544–1551, 2002. View at: Publisher Site  Google Scholar
 R. G. Steen, B. M. Koury, C. I. Granja et al., “Effect of ionizing radiation on the human brain: white matter and gray matter ${\text{T}}_{1}$ in pediatric brain tumor patients treated with conformal radiation therapy,” International Journal of Radiation Oncology Biology Physics, vol. 49, no. 1, pp. 79–91, 2001. View at: Publisher Site  Google Scholar
 R. G. Steen, C. Mull, R. McClure, R. M. Hamer, and J. A. Lieberman, “Brain volume in firstepisode schizophrenia: systematic review and metaanalysis of magnetic resonance imaging studies,” The British Journal of Psychiatry, vol. 188, pp. 510–518, 2006. View at: Publisher Site  Google Scholar
 R. G. Steen, R. J. Ogg, W. E. Reddick, and P. B. Kingsley, “Agerelated changes in the pediatric brain: quantitative MR evidence of maturational changes during adolescence,” American Journal of Neuroradiology, vol. 18, no. 5, pp. 819–828, 1997. View at: Google Scholar
 R. G. Steen and J. Schroeder, “Agerelated changes in the pediatric brain: proton ${\text{T}}_{1}$ in healthy children and in children with sickle cell disease,” Magnetic Resonance Imaging, vol. 21, no. 1, pp. 9–15, 2003. View at: Publisher Site  Google Scholar
 D. L. Pham and J. L. Prince, “Adaptive fuzzy segmentation of magnetic resonance images,” IEEE Transactions on Medical Imaging, vol. 18, no. 9, pp. 737–752, 1999. View at: Publisher Site  Google Scholar
 S. Ruan, B. Moretti, J. Fadili, and D. Bloyet, “Fuzzy Markovian segmentation in application of magnetic resonance images,” Computer Vision and Image Understanding, vol. 85, no. 1, pp. 54–69, 2002. View at: Publisher Site  Google Scholar
 K. Held, E. R. Kops, B. J. Krause, W. M. Wells, R. Kikinis, and H.W. MullerGartner, “Markov random field segmentation of brain MR images,” IEEE Transactions on Medical Imaging, vol. 16, no. 6, pp. 878–886, 1997. View at: Google Scholar
 J. L. Marroquin, B. C. Vemuri, S. Botello, F. Calderon, and A. FernandezBouzas, “An accurate and efficient Bayesian method for automatic segmentation of brain MRI,” IEEE Transactions on Medical Imaging, vol. 21, no. 8, pp. 934–945, 2002. View at: Publisher Site  Google Scholar
 Y. Zhang, M. Brady, and S. Smith, “Segmentation of brain MR images through a hidden Markov random field model and the expectationmaximization algorithm,” IEEE Transactions on Medical Imaging, vol. 20, no. 1, pp. 45–57, 2001. View at: Publisher Site  Google Scholar
 D. Mumford and J. Shah, “Optimal approximations by piecewise smooth functions and associated variational problems,” Communications on Pure and Applied Mathematics, vol. 42, 1989. View at: Google Scholar
 H.L. M. Cheng and G. A. Wright, “Rapid highresolution ${\text{T}}_{1}$ mapping by variable flip angles: accurate and precise measurements in the presence of radiofrequency field inhomogeneity,” Magnetic Resonance in Medicine, vol. 55, no. 3, pp. 566–574, 2006. View at: Google Scholar
 S. C. L. Deoni, B. K. Rutt, and T. M. Peters, “Rapid combined ${\text{T}}_{1}$ and ${\text{T}}_{2}$ mapping using gradient recalled acquisition in the steady state,” Magnetic Resonance in Medicine, vol. 49, no. 3, pp. 515–526, 2003. View at: Publisher Site  Google Scholar
 S. C. L. Deoni, T. M. Peters, and B. K. Rutt, “Determination of optimal angles for variable nutation proton magnetic spinlattice, ${\text{T}}_{1}$, and spinspin, ${\text{T}}_{2}$, relaxation times measurement,” Magnetic Resonance in Medicine, vol. 51, no. 1, pp. 194–199, 2004. View at: Publisher Site  Google Scholar
 D. Wang, K. Heberlein, S. LaConte, and X. Hu, “Inherent insensitivity to RF inhomogeneity in FLASH imaging,” Magnetic Resonance in Medicine, vol. 52, no. 4, pp. 927–931, 2004. View at: Publisher Site  Google Scholar
 E. M. Haacke, R. W. Brown, M. R. Thompson, and R. Venkatesan, Magnetic Resonance Imaging: Physical Principles and Sequence Design, John Wiley & Sons, New York, NY, USA, 1999.
 R. Venkatesan, W. Lin, and E. M. Haacke, “Accurate determination of spindensity and ${\text{T}}_{1}$ in the presence of RF field inhomogeneities and flipangle miscalibration,” Magnetic Resonance in Medicine, vol. 40, no. 4, pp. 592–602, 1998. View at: Publisher Site  Google Scholar
 R. G. Steen, W. E. Reddick, and R. J. Ogg, “More than meets the eye: significant regional heterogeneity in human cortical ${\text{T}}_{1}$,” Magnetic Resonance Imaging, vol. 18, no. 4, pp. 361–368, 2000. View at: Publisher Site  Google Scholar
 M. Kass, A. Witkin, and D. Terzopoulos, “Snakes: active contour models,” International Journal of Computer Vision, vol. 1, no. 4, pp. 321–331, 1988. View at: Publisher Site  Google Scholar
 L. D. Cohen and I. Cohen, “Finiteelement methods for active contour models and balloons for 2D and 3D images,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 15, no. 11, pp. 1131–1147, 1993. View at: Publisher Site  Google Scholar
 V. Caselles, “Geometric models for active contours,” in Proceedings of the IEEE International Conference on Image Processing, vol. 3, p. 3009, 1995. View at: Google Scholar
 V. Caselles, F. Catté, T. Coll, and F. Dibos, “A geometric model for active contours in image processing,” Numerische Mathematik, vol. 66, no. 1, pp. 1–31, 1993. View at: Publisher Site  Google Scholar  MathSciNet
 L. D. Cohen, “On active contour models and balloons,” CVGIP: Image Understanding, vol. 53, no. 2, pp. 211–218, 1991. View at: Google Scholar
 S. Kichenassamy, A. Kumar, P. Olver, A. Tannenbaum, and A. Yezzi Jr., “Conformal curvature flows: from phase transitions to active vision,” Archive for Rational Mechanics and Analysis, vol. 134, no. 3, pp. 275–301, 1996. View at: Google Scholar
 R. Malladi, J. A. Sethian, and B. C. Vemuri, “Shape modeling with front propagation: a level set approach,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 17, no. 2, pp. 158–175, 1995. View at: Publisher Site  Google Scholar
 H. Tek and B. B. Kimia, “Image segmentation by reactiondiffusion bubbles,” in Proceedings of the IEEE International Conference on Computer Vision, pp. 156–162, 1995. View at: Google Scholar
 D. Terzopoulos, A. Witkin, and M. Kass, “Constraints on deformable models: recovering 3D shape and nonrigid motion,” Artificial Intelligence, vol. 36, no. 1, pp. 91–123, 1988. View at: Google Scholar
 A. Yezzi Jr., S. Kichenassamy, A. Kumar, P. Olver, and A. Tannenbaum, “A geometric snake model for segmentation of medical imagery,” IEEE Transactions on Medical Imaging, vol. 16, no. 2, pp. 199–209, 1997. View at: Google Scholar
 C. Xu and J. L. Prince, “Generalized gradient vector flow external forces for active contours,” Signal Processing, vol. 71, no. 2, pp. 131–139, 1998. View at: Google Scholar
 T. F. Chan and L. A. Vese, “Active contours without edges,” IEEE Transactions on Image Processing, vol. 10, no. 2, pp. 266–277, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 A. Yezzi, A. Tsai, and A. Willsky, “A fully global approach to image segmentation via coupled curve evolution equations,” Journal of Visual Communication and Image Representation, vol. 13, no. 12, pp. 195–216, 2002. View at: Publisher Site  Google Scholar
 S. C. Zhu and A. Yuille, “Region competition: unifying snakes, region growing, and bayes/mdl for multiband image segmentation,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 18, no. 9, pp. 884–900, 1996. View at: Google Scholar
 R. Ronfard, “Regionbased strategies for active contour models,” International Journal of Computer Vision, vol. 13, no. 2, pp. 229–251, 1994. View at: Publisher Site  Google Scholar
 C. Samson, L. BlancFeraud, G. Aubert, and J. Zerubia, “A level set model for image classification,” International Journal of Computer Vision, vol. 40, no. 11, pp. 187–197, 2000. View at: Google Scholar
 M. Tang and S. Ma, “General scheme of region competition based on scale space,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 23, no. 12, pp. 1366–1378, 2001. View at: Publisher Site  Google Scholar
 G. Unal, A. Yezzi, and H. Krim, “Informationtheoretic active polygons for unsupervised texture segmentation,” International Journal of Computer Vision, vol. 62, no. 3, pp. 199–220, 2005. View at: Publisher Site  Google Scholar
 D. Cremers, T. Kohlberger, and C. Schnörr, “Nonlinear shape statistics in mumfordshah based segmentation,” in Proceedings of the 7th European Conference on Computer VisionPart II (ECCV '02), pp. 93–108, Springer, London, UK, 2002. View at: Google Scholar
 N. Paragios and R. Deriche, “Geodesic active regions and level set methods for supervised texture segmentation,” International Journal of Computer Vision, vol. 46, no. 3, pp. 223–247, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 J. Kim, J. W. Fisher III, A. Yezzi, M. Cetin, and A. S. Willsky, “A nonparametric statistical method for image segmentation using information theory and curve evolution,” IEEE Transactions on Image Processing, vol. 14, no. 10, pp. 1486–1502, 2005. View at: Publisher Site  Google Scholar
 G. Unal, H. Krim, and A. Yezzi, “Fast incorporation of optical flow into active polygons,” IEEE Transactions on Image Processing, vol. 14, no. 6, pp. 745–759, 2005. View at: Publisher Site  Google Scholar
 A. Stuart and J. K. Ord, Kendall's Advanced Theory of Statistics. Vol.1: Distribution Theory, Hodder Arnold, London, UK, 6th edition, 1994.
 A. Tsai, A. Yezzi Jr., and A. S. Willsky, “Curve evolution implementation of the MumfordShah functional for image segmentation, denoising, interpolation, and magnification,” IEEE Transactions on Image Processing, vol. 10, no. 8, pp. 1169–1186, 2001. View at: Publisher Site  Google Scholar  Zentralblatt MATH
 C. V. Alvino and A. J. Yezzi, “Fast MumfordShah segmentation using image scale space bases,” in Computational Imaging V, vol. 6498 of Proceedings of SPIE, pp. 1–10, San Jose, Calif, USA, January 2007. View at: Publisher Site  Google Scholar
 A. Yezzi, L. Zollei, and T. Kapur, “A variational framework for joint segmentation and registration,” in Proceedings of the Workshop on Mathematical Methods in Biomedical Image Analysis (MMBIA '01), pp. 44–51, 2001. View at: Google Scholar
 H. Z. Wang, S. J. Riederer, and J. N. Lee, “Optimizing the precision in ${\text{T}}_{1}$ relaxation estimation using limited flip angles,” Magnetic Resonance in Medicine, vol. 5, no. 5, pp. 399–416, 1987. View at: Google Scholar
 W. Lin, R. P. Paczynski, R. Venkatesan et al., “Quantitative regional brain water measurement with magnetic resonance imaging in a focal ischemia model,” Magnetic Resonance in Medicine, vol. 38, no. 2, pp. 303–310, 1997. View at: Publisher Site  Google Scholar
 R. G. Steen, S. A. Gronemeyer, P. B. Kingsley, W. E. Reddick, J. S. Langston, and J. S. Taylor, “Precise and accurate measurement of proton ${\text{T}}_{1}$ in human brain in vivo: validation and preliminary clinical application,” Journal of Magnetic Resonance Imaging, vol. 4, no. 5, pp. 681–691, 1994. View at: Google Scholar
 P. B. Kingsley, R. J. Ogg, W. E. Reddick, and R. G. Steen, “Correction of errors caused by imperfect inversion pulses in MR imaging measurement of ${\text{T}}_{1}$ relaxation times,” Magnetic Resonance Imaging, vol. 16, no. 9, pp. 1049–1055, 1998. View at: Publisher Site  Google Scholar
 M. Atkins and B. T. Mackiewich, “Fully automatic segmentation of the brain in MRI,” IEEE Transactions on Medical Imaging, vol. 17, no. 1, pp. 98–107, 1998. View at: Google Scholar
 Z. Y. Shan, G. H. Yue, and J. Z. Liu, “Automated histogrambased brain segmentation in ${\text{T}}_{1}$weighted threedimensional magnetic resonance head images,” NeuroImage, vol. 17, no. 3, pp. 1587–1598, 2002. View at: Publisher Site  Google Scholar
 BrwainWeb, “Mcconnell brain imaging center, montreal neurological institute,” http://www.bic.mni.mcgill.ca/brainweb/. View at: Google Scholar
 Y. He, A. B. Hamza, and H. Krim, “A generalized divergence measure for robust image registration,” IEEE Transactions on Signal Processing, vol. 51, no. 5, pp. 1211–1220, 2003. View at: Publisher Site  Google Scholar  MathSciNet
 L. A. Vese and T. F. Chan, “A multiphase level set framework for image segmentation using the Mumford and Shah model,” International Journal of Computer Vision, vol. 50, no. 3, pp. 271–293, 2002. View at: Publisher Site  Google Scholar  Zentralblatt MATH
Copyright
Copyright © 2009 PingFeng Chen et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.