100	D.P. Costello and P.A. Kenny

Fig. 5.7 Images and their corresponding intensity histograms obtained using T1W TSE and WS b-SSFP sequences. Both images are acquired from the same anatomic slice of the same subject with breath hold. The T1W TSE image (a) is shown to have lower contrast between fat and nonfat. Water and partial-volume fat signal are also in the same gray level range as shown in (b), making automated fat quantification difficult. The WS b-SSFP image (d), however, shows negligible water signal, leading to fat-only images. The corresponding histogram (e) shows that signal from suppressed water, partial-volume fat, and full-volume fat are delineated. This makes it possible to perform automated, yet accurate fat quantification. This material is reproduced with permission of John Wiley & Sons, Inc. [44] with the addition of images (c) and (f). Image (c) is the result of OTSU thresholding applied to (a) and (f) is the result of the segmentation technique described by Peng et al. [44]

5.4.3 Gaussian Mixture Model

When overlapping peaks cannot be avoided the Gaussian mixture model (GMM) can be used to model complex probability density functions (PDF), such as image histograms, f (x), as k overlapping Gaussians. This is illustrated in Fig. 5.8 and can be expressed as:

k
f (x) = ∑ piN x|μi, σi2 ,	(5.6)

i=1

where μi and σi are the mean and standard deviation of the ith Gaussian, respectively. pi is the mixing proportion of the constituent Gaussians, N, used to model the PDF of the image. It satisfies the conditions:

k
pi > 0, ∑ pi = 1.	(5.7)

i=1

5 Fat Segmentation in Magnetic Resonance Images

101

Fig. 5.8 (a) T1-weighted GE image containing fat and soft tissue, (b) shows the threshold value T illustrates the threshold location using a red line. Voxels with gray level intensities higher than T are classified as fat. (c) Image histogram of (a) and its constituent Gaussians estimated using the GMM

Each Gaussian cluster is then modeled using the product of (5.8), the general equation for a Gaussian distribution and its probability (pi).

N(	,	σi	) =	1		exp		(x − μi)2	.	(5.8)
	μi			σi√2π			−	2σi2

Figure 5.8c, illustrates an image histogram containing 3 tissue classes and its constituent Gaussians calculated using a GMM. One of the most common algorithms used to calculate pi, μi and σi is the expectation maximization (EM) algorithm [45].

The GMM assumes that the PDF of all constituent tissue types in an image histogram are Gaussian and that tissue of the same class is uniform in intensity throughout the image or region of image to be segmented. When segmenting fat in MR images, k is usually estimated using a priori knowledge as 3 (fat, soft tissue and background).

The probability of an event (pixel/voxel) belonging to the ith distribution is given by

k
P(x|Θ) = ∑ piN(x|μi , σi ),	(5.9)

i=1

102	D.P. Costello and P.A. Kenny
where
	Θ = {μi=1, . . . , μi=k , σi=1, . . . , σi=k , pi=1, . . . , pi=k }.	(5.10)

The EM algorithm is a two-step iterative algorithm consisting of an expectation step (E-step) and a maximization step (M-step). The algorithm can be initialized using k- means clustering or by equal partitioning of the data into k regions or by automated seed initialization [46]. The expected log likelihood function for the complete data

is calculated using the E-step and is defined by Q(Θ, Θ(t)) using the estimated
parameters Θ (t).
		(5.11)
Q(Θ, Θ(t)) ≡ E[logN (X , Y |Θ)|X ,Θ (t)].

The function Q(Θ, Θ(t)) contains two arguments Θ denotes the parameters that

will be optimized in order to maximize the likelihood and Θ(t) corresponds to the estimated values. X is the observed data and remains constant while Y is the missing data, which is controlled by the underlying distributions.

The second step in the algorithm is the M-step. This step uses the maximized values from (5.11) above to generate a new set of parameters and is given by:

		(5.12)
Θ(t + 1) = arg max Q(Θ, Θ(t)).
	Θ

If we have an estimate of the means (μi) and standard deviation (σi) of the constituent Gaussians we can compute the probability (pi) of a point (gray-level) in the histogram belonging to the kth Gaussian. The maximization step updates the Gaussian parameters using (5.13), (5.14) and (5.15).

pinew =	1k		k
			∑ N( i|x, Θ(t)).
			i=1
μ new =	∑ik=1 xN( i|x, Θ(t)) .
i
	∑ik=1 N( i|x, Θ(t))

σ new =
	∑ik=1 N( i|x, Θ(t))(xi − μinew)(xi − μinew)T .
i
	∑ik=1 N( i|x, Θ(t))

(5.13)

(5.14)

(5.15)

The algorithm iterates between (5.11) and (5.12) until convergence is reached. Based on the Gaussian estimated using the EM algorithm an optimized threshold is calculated that minimizes misclassification error. The threshold for fat is set to the closest gray-level corresponding to the intersection of the fat and soft tissue Gaussians, as illustrated in Fig. 5.8b. Figure 5.13 in Sect. 5.5 is an example of an image which has been segmented successfully using the GMM. Lynch et al. [46]

5 Fat Segmentation in Magnetic Resonance Images

103

Fig. 5.9 Seed point selection voxels in a region of interest. Each of the nearest neighbor voxels are illustrated using yellow arrows

proposed a novel approach to initialize cluster centers based on histogram analysis. It used the GMM to obtain a more robust segmentation of fat in MR images. No objective assessment of this method was carried out.

5.4.4 Region Growing

Region growing is an image segmentation technique that creates regions based some predefined homogeneity criteria such as texture, gray-level or color. Gray-level is the characteristic most commonly used when segmenting fat in MR images. Region growing aims to merge voxels or small regions in an image into larger ones based on a homogeneity criteria. The first step in any region-growing algorithm is the selection of a seed point, as illustrated in Fig. 5.9. This can be a single voxel or a small region of voxels and can be selected manually or automatically. Next, a homogeneity criterion is set. For example, if the gray-level of the neighboring voxel is between two threshold values, merge the voxel (region) with the seed point. The seed voxels nearest neighbors are illustrated using the yellow arrows in Fig. 5.9. Each new voxel in the region becomes a growth point and is compared to its nearest neighbors using the same homogeneity criteria. Growth of the region ceases when

104	D.P. Costello and P.A. Kenny

no new growth points can be created within the confines of the homogeneity criteria. Growth or merging of regions will continue iteratively until the stopping criterion is reached.

Region Growing Algorithm

•Select a seed point/points

•Define a homogeneity/merging criteria

•Join all voxels connected to the seed that follow the homogeneity criteria to form a region

•Stop the algorithm when no adjacent voxels agree with the homogeneity criteria

Figure 5.10 shows an image containing a number of manually selected seed points and the resultant segmentation using region growing. Unlike thresholding, region growing can differentiate objects based on their spatial location, which enables the classification of different fat categories. In Fig. 5.10b, note also that bone marrow adipose tissue and the liver are not classified as fat. Fat classification as distinct from segmentation is an important issue and is discussed in detail in Sect. 5.5.

Region growing can be implemented using a number of algorithms, including region merging, region splitting and split and merge algorithms [39, 47]. Split and merge algorithms can be used to automate region growing, removing the need for subjective seed point selection [48]. Automated seed selection does not always guarantee complete segmentation due to the unpredictability of seed point selection and may therefore require manual addition of extra seed points after initial segmentation. In noisy images, it is possible for holes to appear in the segmented regions, post processing may be used to reduce these.

Siegel et al. [6] used a region growing technique in its simplest form to segment fat in transverse MR images of the abdomen. However, region growing is sometimes combined with edge detection to reduce the occurrence of overand under-segmentation. In a study by Yan et al. [49], a three-step segmentation process to isolate skeletal structures in CT images was employed. The first step was the application of a three-dimensional region-growing algorithm with adaptive local thresholds, which was followed by a series of boundary correction steps. This approach is both accurate and reliable for high contrast, low noise CT images. Subsequently, a similar technique was applied to the segmentation of fat in MR images by Brennan et al. [15].

Brennan et al. [15], used a four-step automated algorithm for the quantification of fat in MR images. Their algorithm was based on initialization of fat regions using conservative thresholding followed by steps of boundary enhancement, region growing and region refining (post processing). A weak correlation was found between total body fat quantified using MRI and BMI. The authors attribute this to the flawed nature of BMI measurements. A more appropriate measure of accuracy would have been comparison with manual delineation of fat by an expert radiologist (the gold standard). Combination of both region growing and edge–based segmentation algorithms contradicts work by Rajapakse and Kruggel [23], who state that the use of region and edge detection schemes are unsuitable in MR images due