Mixture of Gaussians Model

  

Figure: LEFT: Joint intensity histogram of the registered MR PD/T2 training images used to fit a mixture of Gaussians. RIGHT: Rough segmentation of a registered image pair. Each voxel is classified based on which Gaussian it most likely belongs to, based on the mixture of Gaussian model.

Given a pair of registered images from two different medical image acquisitions, we can assume that each voxel with coordinate x = [x₁, x₂, x₃]^T in one image, I₁, corresponds to the same position in the patient's anatomy as the voxel with coordinate x in the other image, I₂. Further, consider that the anatomical structure S_k at some position in the patient will appear with some intensity value i₁ in the first image and i₂ in the second image with joint probability $P(i_1, i_2 \,\vert\, S_k)$. We also define $P(S_k) = \pi_k$ to be the prior probability that a random point in the medical scan corresponds to structure S_k.

By making the assumption that voxels are independent samples from this distribution (and ignoring relative positions of voxels), we have

 

$$\prod_{x \in I} P(I_1(x), I_2(x))$$ P(I₁, I₂) = (2)

$$\prod_{x \in I} \sum_{k} \pi_k P(I_1(x), I_2(x) \,\vert\, S_k)$$ = (3)

We model the joint intensity of a particular internal structure S_kto be a two dimensional (dependent) Gaussian with mean $\mu_k$ and full covariance matrix $\Sigma_k$. Letting i be intensity pair [i₁, i₂]^T,

$$P(i \,\vert\, S_k) = \left( \frac{1}{2 \pi \vert\Sigma_k\vert... ...}} e^{-\frac{1}{2}(i-\mu_k)^T \Sigma_k^{-1} (i-\mu_k)} \right)$$ (4)

This model of the intensities corresponds to a mixture of Gaussians distribution, where each 2D Gaussian G_k corresponds to the joint intensity distribution of an internal anatomical structure S_k. Thus, the probability of a certain intensity pair (independent of the anatomical structure) given the model, M is

$$P(i \,\vert\, M) = \sum_{k} \left( \frac{\pi_k}{2 \pi \vert\S... ...} e^{-\frac{1}{2}(i-\mu_k)^T \Sigma_k^{-1} (i-\mu_k)} \right).$$ (5)

To learn the joint intensity distribution under this model, we estimate the parameters $\pi_k$, $\mu_k$, and $\Sigma_k$ using the Expectation-Maximization (EM) method [3].

The mixture of Gaussians model was chosen to represent the joint intensity distribution because we are imaging a volume with various anatomical structures that respond with different ranges of intensity values in the two acquisitions. We assume that those ranges of responses are approximately Gaussian in nature. Therefore, one might expect that each Gaussian in the mixture may correspond roughly to one type of anatomical structure. In other words, the model produces an approximate segmentation of the structures in the images. Figure 3b shows the segmentation of a registered pair of MR images using the Gaussian mixture model prior. Gerig, et al. [5] used similar methods of statistical classification to produce accurate unsupervised segmentation of 3D dual-echo MR data.

Segmentation of medical images based solely on intensity classification (without using position or shape information) is, in general, very difficult. Often different tissue types may produce a similar or overlapping range of intensity responses in a given medical scan, making classification by intensity alone quite challenging. MR images include nonlinear gain artifacts due to inhomogeneities in the receiver or transmit coils [4]. Furthermore, the signal can also be degraded by motion artifacts from movement of the patient during the scan.

The segmentation produced by this method shown in Figure 3b suffers from the difficulties described above. For example white matter and gray matter have overlapping ranges of intensities in both image acquisitions. Furthermore, note that the distinction between gray and white matter on the right hand side is not segmented clearly. This is most likely due to the bias field present in the image.

  

Figure: Two views of the joint intensity distribution function computed using Parzen estimation with a Gaussian windowing function.

Despite these difficulties, the segmentation does a reasonable job of picking out the major structures, although it is inaccurate at region boundaries. Therefore, we do not intend to use this method alone to compute an accurate segmentation of the underlying structures. Instead, we could use the mixture model in combination with a more sophisticated algorithm to solve for segmentation, or for registration purposes, as described in section 3.