Bias Field Correction through Modified EM Algorithm

One of the most successful methods for dealing with the bias field problem was developed by Wells et al. [31], in which the bias field $B=(b_1, \dots, b_N)$ is modelled as a multiplicative N-dimensional random vector with zero mean Gaussian prior probability density $p(B)=G_{\psi_B}(B)$, where $\psi_B$ is the $N \times N$ covariance matrix. Let $I=(I_1,\cdots, I_N)$ and $I^*=(I_1^*,\cdots, I_N^*)$ be the observed and the ideal intensities of a given image respectively. The degradation effect of the bias field at pixel $i, 1\leq i \leq N$ can be expressed as follows:

$$I_i = I_i^* \times b_i,$$ (29)

After logarithmic transformation of the intensities, the bias field effect can be treated as an additive artifact. Let y and y^* denote respectively the observed and the ideal log-transformed intensities: then y = y^* + B. Given the class labels x, it is further assumed that the ideal intensity value at pixel i follows a Gaussian distribution with parameter $\theta(x_i)=(\mu_{x_i} \sigma_{x_i})$:

$$p(y_i^*\vert x_i)=g(y_i^*;\theta(x_i)).$$ (30)

With the bias field b_i taken into account, the above distribution can be written in terms of the observed intensity y_i as

$$p(y_i\vert x_i,B)=g(y_i-b_i;\theta(x_i)) .$$ (31)

Thus, the intensity distribution is modelled as a Gaussian mixture, given the bias field. It follows that

$$p(y_i\vert B)=\sum_{j \in \mathcal L}\big\{g(y_i-b_i;\theta(j))P(j)\big\}.$$ (32)

The MAP principle is then employed to obtain the optimal estimate of the bias field, given the observed intensity values:

$$\hat B=\arg\max_B p({\mathbf y}\vert B)p(B),$$ (33)

A zero-gradient condition is then used to assess this maximum, which leads to (see [31] for detail):

  

$$\frac{p(y_i\vert x_i,\beta)p(x_i)}{p(y_i\vert\beta)}$$ W_ij = (34)

$$\frac{[FR]_i}{[F\psi^{-1}1]_i}, \text{ with } 1=(1,1,\cdots,1)^T,$$ b_i = (35)

where R is the mean residual for pixel i

$$R_i=\sum_{j \in \mathcal L}\frac{W_{ij}(y_i-\mu_j)}{\sigma_j^2},$$ (36)

$\psi$ is the mean inverse covariance

$$\psi_{ik}^{-1}=\left\{\begin{array}{ll} \sum_{j \in \math... ...& \text{if }i=k \\ 0 & \text{otherwise} \end{array} \right.$$ (37)

and F is a lowpass filter. W_ij is the posterior probability that pixel i belongs to class j given the bias field estimate. The EM algorithm is applied to Equations (34) and (35). The E step assumes that the bias field is known and calculates the posterior tissue class probability W_ij. In the M step, the bias field B is estimated given the estimated W_ij in the E step. Once the bias field is obtained, the original intensity I^* is restored by dividing I by the inverse log of B. Initially, the bias field is assumed to be zero everywhere. Wells et al.'s algorithm is found to be problematic when there are classes in an image that do not follow a Gaussian distribution. The variance of such a class tends to be very large and consequently the mean can not be considered representative[13]. Such situations are commonly seen in the regions of CSF, pathologies and other non-brain classes. Bias field estimation can be significantly affected by this type of problem. To overcome this problem, Guillemaud and Brady [13] unify all such classes into an outlier class, which is called ``other'', with uniform distribution. Let $\mathcal L_G$ denote the set of labels for Gaussian classes and l<sub>o</sub> the class label for the ``other'' class. The intensity distribution of the image is still a finite mixture except with an additional non-Gaussian class,

$$p(y_i\vert b_i)=\sum_{j \in \mathcal L_G} \big\{g(y_i-b_i;\theta(j))P(j)\big\}+\lambda P(l_o),$$ (38)

where $\lambda$ is the density of the uniform distribution. Due to the large variance of the uniform distribution, the bias field is only estimated with respect to the Gaussian classes. The same iterative EM method can be applied, except for a slight modification to the formulation of mean residual R_i (36)

$$R_i=\sum_{j \in \mathcal L_G}\frac{W_{ij}(y_i-\mu_j)}{\sigma_j^2}.$$ (39)

With such a modification, the performance of the EM algorithm can be significantly improved in certain situations. This approach is referred to as the modified EM (MEM) algorithm.