Retrospective Evaluation and Correction of Intensity Inhomogeneity - Advance Image Processing in Magnetic Resonance Imaging - page 147

These equations can be interpreted as follows (see Figure 5.2 ). Given an

estimate of the inhomogeneity field, a statistical classification of the image

voxels is performed by calculating the posterior tissue-class probabilities p

ik

(Figure 5.2b). From this classification and the Gaussian distribution parame-
ters, a prediction of the MR intensities without the inhomogeneity field is
reconstructed (Figure 5.2c). A residue image

y

is obtained by subtracting this
predicted signal from the original image, giving a local estimate of the inho-
mogeneity field in each voxel (Figure 5.2d), along with a weight image

r

that
reflects confidence of these local estimates (Figure 5.2e). Finally, the inhomo-
geneity field

w

β is estimated by applying the linear operator H , which, as will

be shown below, results in smoothing the residue image r while taking the

weights into account (Figure 5.2f ).

Recall that the exact choice of the covariance matrix

w

, governing the prior
distribution of the inhomogeneity field, has not been specified so far. In general,
this matrix is impracticably large, making the evaluation of Equation 5.6 com-
putationally intractable. In the case of equal covariances

Σ

β

σ 2 for all tissue types,

the confidence weights w are constant across the image, and Wells et al. showed

i

that can be chosen so that the linear operator in Equation 5.6 simplifies to

a shift-invariant linear low-pass filter. In general, however, this is not the case,
and the authors heuristically use a computationally efficient approximation by
low-pass-filtering the voxelwise product of the weights and the residues, and
dividing the result by the low-pass-filtered version of the weights:

Σ H

β

[ ]

[]

FWr
Fw

ˆ

β

i

i

i

with a low-pass filter.

It should be noted that the inhomogeneity field estimation involves knowledge
of the tissue classification (Equation 5.6), but that this tissue classification in turn
requires knowledge of the inhomogeneity field (Equation 5.7). Intuitively, both
equations can be solved simultaneously by iteratively alternating between these two
steps. Indeed, it can be shown that such an iterative approach is an instance of the
so-called expectation-maximization (EM) algorithm [21], often used in estimation
problems where some of the data is “missing.” In this case, the missing data is the
tissue type of each voxel; if these were known, estimation of the inhomogeneity
field would be straightforward. In effect, the algorithm of Wells et al. iteratively
fills in the missing tissue types based on the current inhomogeneity field estima-
tion during the E step (Equation 5.7) and updates the inhomogeneity field accord-

F

ingly during the M step (Equation 5.6); cf. Figure 5.3 . Such an EM scheme

guarantees increasingly better estimates of the inhomogeneity field with respect
to Equation 5.4 at each iteration [22], although there is no guarantee of finding the
global optimum. The iterative EM process is typically started on the E step, using
a flat initial inhomogeneity field, although Wells et al. also reported results by
starting on the M step using equal tissue-class probabilities. The authors reported