Stefano Cerri

Research Fellow @ Martinos

scerri at mgh dot harvard dot edu

Stefano Cerri

Research Fellow @ Martinos

scerri at mgh dot harvard dot edu

Modeling the Bias Field in Brain Segmentation

In a previous post, I described how to obtain brain segmentations using a Gaussian Mixture Model (GMM). One of the limitations of this brain segmentation method is that it doesn’t take into account image artifacts present in the scan. This is quite problematic since most of the MRI scans suffer from an image artifact called bias field.

A bias field artifact is characterized by a smooth variation of intensities in the MRI scan. Its effect is stronger in recent scanners due to the increase in the magnetic field strength (e.g. a 7 Tesla scan has a more marked bias field than a 3 Tesla scan).

If we look at the image and segmentation obtained in the previous post we can see that some areas of the image were misclassified by the GMM due to the effect of the bias field.

Example of areas where the GMM model produces an incorrect segmentation due to the bias field artifact.
Example of areas where the GMM model produces an incorrect segmentation due to the bias field artifact.

It is clear that in order to obtain better segmentations we need to take into account the bias field artifact in our model.

In the following, I’m going to describe how we can model the bias field together with our GMM model. This will allow us to simultaneously segment brain structures and account for the bias field present in the image. Lastly, I’m going to show an example of brain segmentation obtained with this model.

Modeling the bias field

We can model the bias field as a linear combination of spatially smooth basis functions \( \textbf{C} \boldsymbol{\phi}^{i} \) for each voxel \( i \), where

$$ \textbf{C} = \left( \begin{array}{c} \textbf{c}_1^T \\ \vdots \\ \textbf{c}_N^T \end{array} \right) , \quad \textbf{c}_n = \left( \begin{array}{cn} c_{n,1} \\ \vdots \\ c_{n,P} \end{array} \right) \quad \text{and} \quad \boldsymbol{\phi}^{i} = \left( \begin{array}{\phi} \phi_{1}^{i} \\ \vdots \\ \phi_{p}^{i} \end{array} \right) $$

Here \( P \) is the number of basis functions, \( N \) is the number of contrasts and \( \boldsymbol{\phi}_{p}^{i} \) is the basis function \( p \) evaluated at voxel \( i \). These basis functions can choose to be sine or cosine functions as well as uniform B-splines or polynomial functions.

Including the bias field model in the GMM

If we include the bias field model in the GMM, we obtain the following log-likelihood function:

$$ \text{ln} p ( \textbf{D} | \boldsymbol{\mu}, \boldsymbol{\Sigma}, \pi, \textbf{C} ) = \sum_{i}^{I} \text{ln} \{ \sum_{k}^{K} \pi_{k} \mathcal{N} ( \textbf{d}_{i} – \textbf{C} \boldsymbol{\phi}_{i} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k} ) \}. $$

Note that we modeled the bias field artifact as an additive artifact since we are working with log-intensities and the bias field is usually modeled as a multiplicative artifact due to the physics of the MR. (Recall the property of the logarithm: \( \text{ln} (ab) = \text{ln} (a) + \text{ln} (b) \)).

Inferring the bias field and GMM parameters using a GEM algorithm

We can now try to maximize the log-likelihood function as we did in the previous post using an Expectation-Maximization (EM) algorithm. However, the introduction of the bias field model makes the updates of the GMM parameters slightly more complicated: i.e. the M-step of the EM algorithm has no closed-form solution.

To solve this problem we can use a generalized EM algorithm (GEM), where we can update the GMM parameters, while keeping the current bias field coefficient estimates fixed. We can then update these bias field coefficients while keeping the GMM parameters fixed.

The E-step involves computing the responsibilities \( w_{k}^{i} \), defined as:

$$ w_{k}^{i} = \frac{\pi_{k} \mathcal{N}( \textbf{d}_{i} – \textbf{C} \boldsymbol{\phi}_{i} | \boldsymbol{\mu}_{k}, \boldsymbol{\Sigma}_{k})}{\sum_{j}^{I} \pi_{j} \mathcal{N} ( \textbf{d}_{i} – \textbf{C} \boldsymbol{\phi}_{i} | \boldsymbol{\mu}_{j}, \boldsymbol{\Sigma}_{j} )} $$

While the M-step gives these updates for the GMM parameters:

$$ \boldsymbol{\mu}_{k} \gets \frac{\sum_{i}^{I}w_{k}^{i} (\textbf{d}_{i} – \textbf{C} \boldsymbol{\phi}_{i})} {N_{k}}, $$

$$ \boldsymbol{\Sigma}_{k} \gets \frac{1}{N_{k}} \sum_{i}^{I} w_{k}^{i} ( \textbf{d}_{i} – \boldsymbol{\mu}_{k} – \textbf{C} \boldsymbol{\phi}_{i} ) ( \textbf{d}_{i} – \boldsymbol{\mu}_{k} – \textbf{C} \boldsymbol{\phi}_{i} )^{T} \text{ and } \pi_{k} \gets \frac{N_{k}}{I}, $$

where we have defined $$ N_{k} = \sum_{i}^{I}w_{k}^{i}. $$

For the coefficients of the bias field functions we obtain the following update:

$$ \left( \begin{array}{c} \textbf{c}_{1} \\ \vdots \\ \textbf{c}_{N} \end{array} \right) \gets \left( \begin{array} {@{}ccc@{}} \textbf{A}^T \textbf{S}_{1,1} \textbf{A} & \cdots & \textbf{A}^T \textbf{S}_{1,N} \textbf{A} \\ \vdots & \ddots & \vdots \\ \textbf{A}^T \textbf{S}_{N,1} \textbf{A} & \cdots & \textbf{A}^T \textbf{S}_{N,N} \textbf{A} \end{array} \right)^{-1} \left( \begin{array}{@{}c@{}} \textbf{A}^T \left( \sum_{n=1}^N \textbf{S}_{i,n} \textbf{r}_{1,n} \right) \\ \vdots \\ \textbf{A}^T \left( \sum_{n=1}^N \textbf{S}_{N,n} \textbf{r}_{N,n} \right) \end{array} \right).$$

Here

$$ \textbf{A} = \left( \begin{array}{@{}ccc@{}} \phi_{1}^1 & \cdots & \phi_{P}^{1} \\ \vdots & \ddots & \vdots \\ \phi_{1}^I & \cdots & \phi_{P}^I \end{array}\right) $$

is a matrix collecting the \( P \) basis functions evaluated at voxel \( i \), while \( \textbf{S}_{m,n} \) is the diagonal matrix defined as:

$$ \textbf{S}_{m,n} = \text{diag} \left( s_{i}^{m,n} \right), \quad s_{i}^{m,n} = \sum_{k=1}^{K} s_{i,k}^{m,n} \quad \text{and} \quad s_{i,k}^{m,n} = w_{i,k} \left( \boldsymbol{\Sigma}_k^{-1} \right)_{m,n}. $$

where each diagonal element contains the sum over the precision matrix for each mixture component weighted by the responsibility for each voxel.

Finally, \( \textbf{r}_{m,n} = \left( r^1_{m,n} , \cdots , r^I_{m,n} \right) \) is a vector containing the residue image computed as the difference between the original input data and the estimated bias-corrected data. Its elements have the form:

$$ r_{i}^{m,n} = d_{i}^{n} – \frac{\sum_{l=1}^K s_{i,k}^{m,n} \left( \boldsymbol{\mu}_{k} \right)_{n}}{\sum_{l=1}^K s_{i,k}^{m,n}}. $$

Example of brain segmentation and bias field correction

Let’s now have a look at how this model obtains better segmentations compared to the standard GMM model. For the following example, I’m using a generated scan from BrainWeb. Here is the image, where I show a slice in axial, coronal, and sagittal view:

Original image with the bias field artifact.
Original image with the bias field artifact.

We are using \( P=3 \) cosine basis function, initialized with \( C=1 \):

Initialization of the bias field with P=3 cosine functions and with C=1.
Initialization of the bias field with P=3 cosine functions and with C=1.

If we run now the GEM algorithm for 80 iterations we obtain the following segmentations:

Segmentation obtained after 80 iterations of the GEM algorithm using a GMM and a bias field correction model
Segmentation obtained after 80 iterations of the GEM algorithm using a GMM and a bias field correction model

We can see that the segmentation is better than the segmentation obtained without accounting for the bias field artifact.

The original image, the bias-corrected image and the estimated bias field are

Original image, bias-corrected image and estimated bias field.
Original image, bias-corrected image, and estimated bias field.

Final thoughts on modeling the bias field

In this post, we modeled the bias field artifact in a GMM brain segmentation model. This extension helps compute better brain segmentation at a small computational cost.

The bias field artifact is present in almost all MR scans and can alter the segmentation performance of most of the brain segmentation methods. To avoid segmentation errors, the bias field should be removed, either by modeling it during the parameter estimation of the segmentation model or by correcting the image as a pre-processing step.

Code

An implementation of this model is available on my Github repository (link), where you can replicate the same experiment and obtain the same figures. The code is written in a jupyter notebook that you can run on your laptop.

Leave a Reply

Your email address will not be published. Required fields are marked *

Top