Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC3318956

Formats

Article sections

Authors

Related links

IEEE Trans Med Imaging. Author manuscript; available in PMC 2012 April 4.

Published in final edited form as:

Published online 2011 April 19. doi: 10.1109/TMI.2011.2143724

PMCID: PMC3318956

NIHMSID: NIHMS365011

Kio Kim, Piotr A. Habas, Member, IEEE, Vidya Rajagopalan, Julia A. Scott, James M. Corbett-Detig, Francois Rousseau, A. James Barkovich, Orit A. Glenn, and Colin Studholme, Senior Member, IEEE

Kio Kim, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA;

The publisher's final edited version of this article is available at IEEE Trans Med Imaging

See other articles in PMC that cite the published article.

A common solution to clinical MR imaging in the presence of large anatomical motion is to use fast multi-slice 2D studies to reduce slice acquisition time and provide clinically usable slice data. Recently, techniques have been developed which retrospectively correct large scale 3D motion between individual slices allowing the formation of a geometrically correct 3D volume from the multiple slice stacks. One challenge, however, in the final reconstruction process is the possibility of varying intensity bias in the slice data, typically due to the motion of the anatomy relative to imaging coils. As a result, slices which cover the same region of anatomy at different times may exhibit different sensitivity. This bias field inconsistency can induce artifacts in the final 3D reconstruction that can impact both clinical interpretation of key tissue boundaries and the automated analysis of the data. Here we describe a framework to estimate and correct the bias field inconsistency in each slice collectively across all motion corrupted image slices. Experiments using synthetic and clinical data show that the proposed method reduces intensity variability in tissues and improves the distinction between key tissue types.

Ultrafast multi-slice imaging sequences, such as single shot fast spin echo (SSFSE) [1] or half-Fourier acquisition single shot turbo spin echo (HASTE) [2] are increasingly popular in clinical imaging of moving anatomy, allowing the clinician to view 2D slices of anatomy in challenging clinical applications such as in utero fetal brain studies [3], [4]. However, accurate 3D representation or measurement of the anatomy is not possible from these studies due to between-slice motion. Recently, techniques that combine multiple multi-slice acquisitions into a single geometrically correct volumetric image via full 3D slice motion estimation have emerged as a viable route to true 3D imaging of moving anatomy. These techniques estimate full 3D motion of a region of anatomy within each individual slice by a retrospective image registration, followed by reconstruction into a true volumetric 3D image from the spatially scattered slice data. These techniques make use of the common clinical practice of planning multiple approximately orthogonal stacks of slices as part of an imaging study. These views in different axes can be used to provide complementary geometric constraints to correct relative slice positioning in 3D. (See Fig. 1.)

(left) MR image stacks of a human fetal brain, planned and acquired *in utero* in axial, sagittal and orientations. (right) A full 3D volume of the fetal brain is reconstructed from the 2D slice image stacks, after the motion of the local rigid anatomy **...**

The first 3D motion corrected fetal brain imaging approach was described by Rousseau et al., who proposed an iterative slice to volume matching process [5], [6]. Here normalized mutual information (NMI) was used to iteratively refine the alignment of individual slices to a previously estimated 3D volume reconstruction. The repeated steps of volume reconstruction interleaved with slice alignment estimation form a refined 3D volume. An alternative approach, termed the Slice Intersection Motion Correction (SIMC) [7], [8], avoids the intermediate 3D reconstruction step by considering slice to slice matching directly, to improve the final subvoxel alignment accuracy and reduce the computational expense of repeated volume reconstruction. This is achieved by minimising the intersection mismatch of all possible combinations of intersecting slice pairs within an efficient least squares minimization framework.

Both of these approaches then make use of a final reconstruction step which combines and resamples the spatially scattered slice data onto a regular voxel lattice to form a single 3D image through various forms of interpolation [6], [9] or super-resolution [10], [11]. The key here is that to combine the images all methods assume that the contrast in each slice is consistent. However, it is common in many imaging situations, including abdominal imaging used for fetal studies, to have significant bias field inhomogeneity across the slices. The cause of this bias field inhomogeneity includes the spatially varying coil sensitivity [12] and the interaction between the imaged object and the scanner [13]. The influence of the inhomogeneous sensitivity field has been reported especially in case of using phased-array coils for the image acquisition [4], [14]. The signal received from the target region will have different levels depending on the positioning of the region with respect to the coil. If this is combined with significant motion of the anatomy of interest, then the intensity bias of the same region of tissue when seen in different slices will not be consistent. In the presence of such bias field inconsistency, the resulting reconstruction of the regions may contain artifacts. An example of this occurring in fetal imaging is illustrated in Fig. 2 which shows a T2 weighted axial slice and a coronal slice through the fetal brain which have visibly inconsistent image brightness as shown in Fig. 2(A). This difference is highlighted in the comparison of the intensity profile along the intersection of the spatially aligned slices in Fig. 2(B).

(A) An axial slice (left panel) and a coronal slice (right panel) with different bias field strength, displayed in the same contrast and scale. The intersection lines of the slice pair are drawn, with markers for relative locations. (B) The intensity **...**

Bias field inhomogeneity is usually surmountable to a human observer, but images with an inhomogeneous bias field often impose challenges in post-acquisition processing, including segmentation, registration and 3D surface extraction [15], [16]. However, in addition, in motion correction of multislice studies, 2D voxels acquired with inconsistent bias field can also produce image artifacts when combined into a 3D volume. For a successful 3D volume reconstruction, it is crucial to ensure that the bias fields of 2D slices are spatially consistent with one another [6].

There has been extensive research into the problem of retrospective bias inhomogeneity correction of MRI [15]–[17]. However, the studies of the bias field correction have focused on the elimination of the bias field inhomogeneity within a single 2D or 3D image, and the number of studies on the correction of bias field inconsistency between multiple images is limited. One related clinical example is the problem of resolving differences between a body coil and a surface coil [18], [19]. In another related problem, a histogram based method was proposed to jointly resolve the different bias fields present between MR images of different subjects [20]. In this work, MR images acquired with bias fields that are inconsistent across the subject populations are corrected by an entropy minimization approach. The most closely related work is the early work of Rousseau et al. on fetal motion correction [6], which incorporated a correction for the bias field inconsistency between slice stacks into a 3D reconstruction process. This approach explicitly assumes that one entire stack of slices has minimal motion and can be used as a 3D ‘bias reference’ to which all other slices are corrected. Although this has the advantage of simplicity, it is suboptimal for correcting bias field inconsistency of 2D slices acquired during motion, as it assumes at least one stack is acquired with no motion and no bias field variability.

In this paper we present a framework to collectively estimate and correct the bias field inconsistency between all individual slices of a multi-stack imaging study that does not assume a motion free 3D estimate of the bias field or underlying anatomy is available. Given a set of starting slice alignment estimates, accounting for relative motion of individual slices, the proposed method then brings the bias field of all slices into mutual agreement so that they can be used to reconstruct a 3D volume image that will be free of bias-related artifacts. The bias field inconsistency correction process is applied after the slice locations are aligned using the SIMC method [8], and before the final 3D volume reconstruction. After the bias field inconsistency resolved reconstruction, the absolute bias field inhomogeneity of the 3D volume can be removed using conventional retrospective correction techniques.

The goal of this work is to estimate the geometrically correct 3D MR image of a region of target anatomy, represented by a vector *z*_{whole} * ^{n}*, where

When the signal from the anatomy is acquired by the scanner, the signal strength in the rigidly moving tissues is influenced by the spatially varying bias field
and current location of the rigid anatomy described by an *n×n* permutation matrix T for rigid transformation;

$${\mathit{z}}_{\text{acq}}=\mathcal{B}{\text{T}}^{\top}{\mathit{z}}_{\text{rigid}},$$

(1)

where
is an *n × n* diagonal matrix.

When a 2D slice image consisting of *m _{i}* pixels,

$${\text{R}}_{i}{\mathit{y}}_{i}={\text{D}}_{i}{\text{T}}_{i}{\text{S}}_{i}{\mathcal{B}}_{i}{{\text{T}}_{i}}^{\top}{\mathit{z}}_{\text{rigid}}+\eta ,$$

(2)

where R* _{i}* is the masking operator
in the

Matrix representations of MR image slice acquisition of two slices, slice *i* and slice *j*, accounting for the rigid motion T_{*}, spatially non-uniform bias field
, slice selection S_{*}, and 2D downsampling D_{*}.

The true 3D image *z*_{rigid} can be then estimated by solving an inverse problem

$${\widehat{\mathit{z}}}_{\text{rigid}}={\left[\begin{array}{c}{\text{D}}_{1}{\text{T}}_{1}{\text{S}}_{1}{\mathcal{B}}_{1}{\text{T}}_{1}^{\top}\\ \vdots \\ {\text{D}}_{N}{\text{T}}_{N}{\text{S}}_{N}{\mathcal{B}}_{N}{\text{T}}_{N}^{\top}\end{array}\right]}^{+}\left[\begin{array}{c}{\text{R}}_{1}{\mathit{y}}_{1}\\ \vdots \\ {\text{R}}_{N}{\mathit{y}}_{N}\end{array}\right],$$

(3)

where superscript + denotes the Moore-Penrose pseudo inverse.

Matrices D* _{i}* and S

$${\widehat{\mathit{z}}}_{\text{rigid}}={\left[\begin{array}{c}{\text{D}}_{1}{\text{T}}_{1}{\text{S}}_{1}{\text{T}}_{1}^{\top}\\ \vdots \\ {\text{D}}_{N}{\text{T}}_{N}{\text{S}}_{N}{\text{T}}_{N}^{\top}\end{array}\right]}^{+}\left[\begin{array}{c}{\text{R}}_{1}{\mathit{y}}_{1}^{\prime}\\ \vdots \\ {\text{R}}_{N}{\mathit{y}}_{N}^{\prime}\end{array}\right].$$

(4)

Since the slice selection and the bias field distortion are point-to-point multiplications, the matrices S* _{i}* and
are diagonal and thus commutative. Therefore, (2) can be rewritten as

$$\begin{array}{l}{\text{R}}_{i}{\mathit{y}}_{i}={\text{D}}_{i}{\text{T}}_{i}{\mathcal{B}}_{i}{\text{S}}_{i}{\text{T}}_{i}^{\top}\phantom{\rule{0.16667em}{0ex}}{\mathit{z}}_{\text{rigid}}+\eta ,\\ =({\text{D}}_{i}{\text{T}}_{i}{\mathcal{B}}_{i}{\text{T}}_{i}^{\top}{\text{D}}_{i}^{\top}){\text{D}}_{i}{\text{T}}_{i}{\text{S}}_{i}{\text{T}}_{i}^{\top}\phantom{\rule{0.16667em}{0ex}}{\mathit{z}}_{\text{rigid}}+\eta ,\\ ={\text{B}}_{i}{\text{D}}_{i}{\text{T}}_{i}{\text{S}}_{i}{\text{T}}_{i}^{\top}\phantom{\rule{0.16667em}{0ex}}{\mathit{z}}_{\text{rigid}}+\eta ,\\ ={\text{R}}_{i}{\text{B}}_{i}{\mathit{y}}_{i}^{\prime}+\eta ,\end{array}$$

(5)

where the 2D bias field operator B* _{i}* is defined by
${\text{B}}_{i}\equiv {\text{D}}_{i}{\text{T}}_{i}{\mathcal{B}}_{i}{\text{T}}_{i}^{\top}{\text{D}}_{i}^{\top}$, which is an

For any given pair of slices *i* and slice *j* that have a line of intersection, it is assumed that the intensities along the intersection should match such that:

$${\left|\right|{\text{T}}_{j}{\text{S}}_{j}{\text{T}}_{j}^{\top}{\text{D}}_{i}^{\top}{\text{R}}_{i}{\mathit{y}}_{i}^{\prime}-{\text{T}}_{i}{\text{S}}_{i}{\text{T}}_{i}^{\top}{\text{D}}_{j}^{\top}{\text{R}}_{j}{\mathit{y}}_{j}^{\prime}\left|\right|}^{2}=0,$$

(6)

where the slice selection operators *S _{i}* and

For the *i*-th image slice, we define a 2D windowing vector *r** _{i}* = diag{R

$${r}_{i}({\mathit{x}}_{i}){y}_{i}({\mathit{x}}_{i})={r}_{i}({\mathit{x}}_{i}){b}_{i}({\mathit{x}}_{i}){y}_{i}^{\prime}({\mathit{x}}_{i})+\eta ,$$

(7)

for the voxel located at *x** _{i}*. We further define a bias field correction term
${b}_{i}^{c}\equiv {b}_{i}^{-1}$ as a functional of the parameter vector

$$\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}\cap {\mathrm{\Omega}}_{j}}{\{{r}_{i}({\mathit{x}}_{i}(\mathit{\omega})){b}_{i}^{c}({\mathit{x}}_{i}(\mathit{\omega});{\mathit{\theta}}_{i}){y}_{i}({\mathit{x}}_{i}(\mathit{\omega}))-{r}_{j}({\mathit{x}}_{j}(\mathit{\omega})){b}_{j}^{c}({\mathit{x}}_{j}(\mathit{\omega});{\mathit{\theta}}_{j}){y}_{j}({\mathit{x}}_{j}(\mathit{\omega}))\}}^{2}=0,$$

(8)

where Ω* _{i}* is the set of 3D points covered by the

In order to find the bias field parameters that minimize the discrepancy of the image intensity of all the slices, we define a global energy function, and minimize it with respect to the bias field parameters. For the energy function, we use the mean square intensity difference between two intersection profiles, evaluated for all orthogonally planned slice pairs. Given 2D slice images *y*_{1}, …, *y _{N}* and bias field correction terms
${b}_{1}^{c},\cdots ,{b}_{N}^{c}$, the energy function is

$$E(\mathrm{\Theta})=\sum _{i=1}^{N}\sum _{j\in {S}_{i}}\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}\cap {\mathrm{\Omega}}_{j}}{r}^{2}(\mathit{\omega})\times {\{{b}_{i}^{c}({\mathit{x}}_{i}(\mathit{\omega}),{\mathit{\theta}}_{i}){y}_{i}({\mathit{x}}_{i}(\mathit{\omega}))-{b}_{j}^{c}({\mathit{x}}_{j}(\mathit{\omega}),{\mathit{\theta}}_{j}){y}_{j}({\mathit{x}}_{j}(\mathit{\omega}))\}}^{2},$$

(9)

where *S _{i}* is the set of slices in the stacks that are orthogonally planned to the

Given that the bias field correction term
${b}_{i}^{c}$ is a linear functional of the bias field parameter vector *θ** _{i}* and the 2D spatial location

$${b}_{i}^{c}(\mathit{x})=\mathit{u}{(\mathit{x})}^{\top}{\mathit{\theta}}_{i},$$

(10)

where ** u**(

$$\begin{array}{l}\mathit{u}(\mathit{x})={[1,{x}_{1},{x}_{2}]}^{\top}\\ {\mathit{\theta}}_{i}={[{\theta}_{i,0},{\theta}_{i,1},{\theta}_{i,2}]}^{\top},\end{array}$$

(11)

for the first degree bias correction model, at the voxel location ** x** = (

$$\begin{array}{l}\mathit{u}(\mathit{x})={[1,{x}_{1},{x}_{2},{x}_{1}^{2},{x}_{2}^{2},{x}_{1}{x}_{2}]}^{\top}\\ {\mathit{\theta}}_{i}={[{\theta}_{i,0},{\theta}_{i,1},{\theta}_{i,2},{\theta}_{i,3},{\theta}_{i,4},{\theta}_{i,5}]}^{\top},\end{array}$$

(12)

for the second degree model. The energy function in (9) then can be rewritten in a quadratic form

$$E(\mathrm{\Theta})=2{\mathrm{\Theta}}^{\top}\mathbf{M}\mathrm{\Theta}.$$

(13)

The bias covariance matrix **M** is a block matrix;

$$\mathbf{M}=\left(\begin{array}{cccc}{\text{M}}_{11}& {\text{M}}_{12}& \dots & {\text{M}}_{1N}\\ {\text{M}}_{21}& {\text{M}}_{22}& \dots & {\text{M}}_{2N}\\ \vdots & \vdots & \vdots & \ddots \\ {\text{M}}_{N1}& {\text{M}}_{N2}& \dots & {\text{M}}_{NN}\end{array}\right),$$

(14)

where

$${\text{M}}_{ij}=\{\begin{array}{l}\sum _{k\in {S}_{i}}\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}\cap {\mathrm{\Omega}}_{k}}{r}^{2}(\mathit{\omega}){y}_{i}^{2}({\mathit{x}}_{i}(\mathit{\omega}))\times \mathit{u}({\mathit{x}}_{i}(\mathit{\omega}))\mathit{u}{({\mathit{x}}_{i}(\mathit{\omega}))}^{\top}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}i=j\\ -\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}\cap {\mathrm{\Omega}}_{j}}{r}^{2}(\mathit{\omega}){y}_{i}({\mathit{x}}_{i}(\mathit{\omega})){y}_{j}({\mathit{x}}_{i}(\mathit{\omega}))\times \mathit{u}({\mathit{x}}_{i}(\mathit{\omega}))\mathit{u}{({\mathit{x}}_{j}(\mathit{\omega}))}^{\top}.\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}i\ne j\end{array}$$

(15)

For a *d*-th degree model, **M** is a square matrix with
${\scriptstyle \frac{1}{2}}(d+1)(d+2)N$ rows and columns.

In order to eliminate undesired degrees of freedom in the spatial intensity distribution, we impose a set of constraints that require the following quantities are preserved;

$$\mathit{\mu}=\sum _{i=1}^{N}\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}}{y}_{i}({\mathit{x}}_{i}(\mathit{\omega}))\mathit{k}(\mathit{v}(\mathit{\omega}))$$

(16)

$$=\sum _{i=1}^{N}\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}}{b}_{i}^{c}({\mathit{x}}_{i}(\mathit{\omega})){y}_{i}({\mathit{x}}_{i}(\mathit{\omega}))\mathit{k}(\mathit{v}(\mathit{\omega})),$$

(17)

where *x** _{i}*(

$$\mathit{k}(\mathit{v})={[1,{v}_{1},{v}_{2},{v}_{3}]}^{\top},$$

(18)

for the first degree model, and

$$\mathit{k}(\mathit{v})={[1,{v}_{1},{v}_{2},{v}_{3},{v}_{1}^{2},{v}_{2}^{2},{v}_{3}^{2},{v}_{1}{v}_{2},{v}_{2}{v}_{3},{v}_{3}{v}_{1}]}^{\top},$$

(19)

for the second degree model, etc. We define the intensity distribution matrix $\mathit{A}={[{\text{A}}_{1}^{\top},\cdots {\text{A}}_{N}^{\top}]}^{\top}$, where

$${\text{A}}_{i}=\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}}{y}_{i}({\mathit{x}}_{i}(\mathit{\omega}))\mathit{u}({\mathit{x}}_{i}(\mathit{\omega})){\mathit{k}}^{\top}(\mathit{v}(\mathit{\omega})).$$

(20)

Equation (17), then, can be also rewritten in a matrix-vector product form,

$$\mathit{\mu}={\mathbf{A}}^{\top}{\mathrm{\Theta}}_{0}={\mathbf{A}}^{\top}\mathrm{\Theta},$$

(21)

where Θ_{0} is the vector of parameters for no bias correction.

The solution of the constrained quadratic equation

$$\widehat{\mathrm{\Theta}}=arg\underset{\mathrm{\Theta}}{min}{\mathrm{\Theta}}^{\top}\mathbf{M}\mathrm{\Theta}$$

(22)

$$\text{subject}\phantom{\rule{0.16667em}{0ex}}\text{to}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}{\mathbf{A}}^{\top}\mathrm{\Theta}=\mathit{\mu}$$

(23)

can be found using the method of Lagrange multipliers [22] where the Lagrangian function is given by

$$L(\mathrm{\Theta},\lambda )={\mathrm{\Theta}}^{\top}\mathbf{M}\mathrm{\Theta}-\lambda ({\mathbf{A}}^{\top}\mathrm{\Theta}-\mathit{\mu}),$$

(24)

and *λ* is the Lagrange multiplier. The solution of the stationary conditions of (24) is

$$\widehat{\mathrm{\Theta}}={\mathbf{M}}^{-1}\mathbf{A}{({\mathbf{A}}^{\top}{\mathbf{M}}^{-1}\mathbf{A})}^{-1}\mathit{\mu}.$$

(25)

Note that (25) is the exact and general solution to Eqs. (22)–(23), for any multiplicative bias field model.

The original image stacks are acquired approximately in axial, sagittal and coronal planes of the subject coordinate system. The multislice MR acquisitions are prepared for alignment using the SLIMMER tool [23]. The motion that occurred during the acquisition of the slices is corrected using the SIMC method, where the fetal motion was estimated by comparing the spatial gradient of the image intensity and the intensity difference between the slices images. Slices corrupted by in-plane motion were excluded automatically from the reconstruction by thresholding the mean signal intensity difference of a slice along the slice intersections, or manually by a visual inspection. After the motion correction, for comparison, a 3D volume was reconstructed without bias correction using the gradient weighted Gaussian averaging [8].

The least square error solution of (22–23) is determined such that it minimizes the square difference of intensity between all slices. However, there are cases that can be mistaken for a bias field mismatch arising from fine scale differences in structure. Examples include partial voluming at sharp intensity boundaries, slight misalignment of images, different noise levels, or subtle differences in image quality. Since such mismatches cannot be resolved by adjusting the bias parameters, the solution can overfit the curve and consequently considers it as a darkened region. In order to avoid such overfitting caused by fine scale differences in structure and leave only the influence of the bias field, the high frequency components of the intersection profiles can be suppressed through low-pass filtering. A 2D low-pass filtering, however, is not legitimate for the purpose, because the 2D filters used for low-pass filtering are defined within individual slice planes that have different orientations. The comparison of the intersection profiles of slice images filtered with different 2D filters is suboptimal. In this consideration, the slice intersections are filtered in 1D along the direction of the intersection line to avoid the overfitting. We defined a modified bias field matrix **M**′, where the entries are redefined from (15) such that;

$${{\text{M}}^{\prime}}_{ij}=\{\begin{array}{l}\sum _{k\in {S}_{i}}\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}\cap {\mathrm{\Omega}}_{k}}{r}^{2}(\mathit{\omega}){\{[{y}_{i}\ast {g}_{ij}^{\sigma}]({\mathit{x}}_{i}(\mathit{\omega}))\}}^{2}\times \mathit{u}({\mathit{x}}_{i}(\mathit{\omega})){\mathit{u}}^{\top}({\mathit{x}}_{i}(\mathit{\omega}))\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}i=j\\ -\sum _{\mathit{\omega}\in {\mathrm{\Omega}}_{i}\cap {\mathrm{\Omega}}_{j}}{r}^{2}(\mathit{\omega})\{[{y}_{i}\ast {g}_{ij}^{\sigma}]({\mathit{x}}_{i}(\mathit{\omega}))\}\times \{[{y}_{j}\ast {g}_{ji}^{\sigma}]({\mathit{x}}_{j}(\mathit{\omega}))\}\mathit{u}({\mathit{x}}_{i}(\mathit{\omega})){\mathit{u}}^{\top}({\mathit{x}}_{j}(\mathit{\omega})),\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}\phantom{\rule{0.38889em}{0ex}}i\ne j\end{array}$$

(26)

where * stands for convolution and
${g}_{ij}^{\sigma}$ a 1D Gaussian filter defined in the *i*-th slice along the intersection with the *j*-th slice, with the standard deviation *σ*.

For the implementation, first and second degree bias correction models are used, as given in (11) and (12). The dimensions of the matrix **M** vary depending on the number of slices in a study. Each entry of **M** is computed using a pair of slice intersection profiles, smoothed with a 1D Gaussian filter as in (26), where *σ* = 7.5 mm was empirically chosen. Since the true masking operator
is unavailable in the clinical setup, we used a parametric model consisting of an ellipsoid within the field of view of the subject as in [8]. The matrix inversion of (25) is computed using the LU decomposition.

The correction for the bias field inconsistency between slices brings the bias field of the individual slices into collective agreement. This allows the artifact-free contributions when combined to reconstruct a 3D volume image. However, the global inhomogeneity of bias field that would arise even with no motion present is still present across the scattered slice data. When the volume image has been reconstructed into a regular 3D voxel lattice it may thus contain significant bias field inhomogeneity. For accurate evaluation and further quantitative analysis, we remove this final absolute bias field inhomogeneity using the automated model-based bias correction method [24], using anatomically specific priors [25], for accurate evaluation and further quantitative analysis.

The proposed method was first evaluated using 50 simulated studies, generated using the Shepp-Logan phantom [26]. The simulation accounted for rigid target motion, inhomogeneous bias field, and partial volume effect. The original phantom was constructed with 256×256×256 voxels, with 0.333mm×0.333mm×0.333 mm voxel dimensions. Each simulated study consisted of 6 stacks of image slices, two stacks per each orientation. A simulated stack was formed consisting of 40 slices, where the voxel dimensions was 0.5mm×0.5 mm in-plane, and the slice thickness was 2.1 mm. Since the simulated slices have greater voxel dimensions than the original phantom, partial volume effect was simulated, assuming a Gaussian point spread function to model the 2D slice selection process. The target object motion was parameterized with three translational and three rotational parameters. The simulated acquisition employed a clinical setup, with an interleave of two, no gap and no overlap between slices. Two motion points were randomly added for each of six motion parameters in each stack, and the motion was temporally smoothed for physical reality. The simulated studies have a varying amount of motion, which was quantified by the slice intersection error before motion estimation [8]. A simulated coil sensitivity field was generated such that the signal is inversely proportional to the distance to a coil, which is located outside the field of view. Although the sensitivity field was static, the motion within the non-uniform sensitivity field induces bias field inconsistency between slices.

Three full 3D volumes were reconstructed for each study using the gradient weighted Gaussian averaging [8], before the bias field inconsistency correction, and after the correction with the first and second degree polynomial models. Figure 5 shows the reconstruction of one of the simulated studies before (A) and after (B) the bias field inconsistency correction using the second degree polynomial model. The reconstruction quality of the 3D volumes improved after the bias field inconsistency correction, which is drastically contrasted in the high contrast view.

Comparison of 3D reconstruction (A) without and (B) with the bias field inconsistency correction, shown in the original contrast (top) and in enhanced contrast (bottom).

The correlation between the amount of bias related reconstruction artifacts and the amount of corrected motion is plotted in Fig. 6. The coefficient of variation (CV; the standard deviation divided by the mean) of the intensity value in the main mid-gray region of the phantom was used to quantify the bias related artifacts [17]. Corrected motion was quantified by the root mean square (RMS) of intersection error (IE) [8]. CV of the signal intensity in the gray region of the volume reconstructed before the bias field inconsistency correction (marked by black diamonds) demonstrates the tendency of increasing bias field inconsistency with increasing motion. The CV appears linearly dependent on RMSIE in the plot, with the rate of 0.10 % per mm (*r*^{2}=0.80). When the bias field inconsistency was corrected, on the other hand, the CV of the reconstruction with severe motion recovery (higher RMSIE) was comparable to the CV of lower RMSIE studies, for both first and second degree polynomial models. The results clearly demonstrate that the proposed method effectively removes reconstruction artifacts caused by the bias field inconsistency between the image slices, and the first degree polynomial model is sufficient for the Shepp-Logan experiment.

The 18 subjects were healthy fetuses aged between 21 and 27 gestational weeks, as measured from the last menstrual period (LMP). All had MR imaging as part of a study of fetuses with mild ventriculomegaly, which was approved by our Institutional Review Board; written consent was obtained from the mother in all cases. No sedation was used during the examination. All scans were interpreted as normal by a pediatric neuroradiologist. Images were acquired using a single shot fast spin-echo (SSFSE) technique on a 1.5T MR scanner (Signa, GE Healthcare, Milwaukee, WI) using a torso phased array coil. For each study, slice stacks were planned in axial, sagittal, and coronal planes using voxel dimensions of 0.5mm×0.5mm×3 mm. Across the different subject studies, both the repetition time (TR) 6400±2600 ms and the echo time (TE) 90.6±0.6 ms varied. The motion between the acquisition of slices was corrected using a slice intersection motion correction method [7] and the full 3D volume was reconstructed using the motion corrected slices [8].

In order to quantify the improvement of the image quality, the 3D volume was reconstructed in three ways: First, the 3D volume reconstruction method in [8] was performed without the bias field inconsistency correction step, second, the reconstruction was done after the bias field inconsistency correction with the first degree polynomial model. and third, the reconstruction was done after the correction with the second degree polynomial model. All reconstructed 3D volumes were then corrected for bias field inhomogeneity using the model-based bias correction algorithm described in [24] and adapted for fetal imaging [27]. This used four pre-segmented labels—gray matter, fetal white matter (including subplate (SP) and intermediate zone (IZ)), germinal matrix, and cerebrospinal fluid (CSF).

The energy function in Eq. (9) was evaluated before and after the bias field inconsistency correction. It was then normalized and plotted against the gestational age in Fig. 7. The plot demonstrates that the amount of initial bias field inhomogeneity has no dependency on the gestational age of the fetus. In all cases, the proposed method improved the amount of bias field mismatch. The second degree model provided a small but consistent improvement over the first degree model. The energy value after the bias field inconsistency correction varied across the subjects, this may be attributable to different noise levels, within-slice motion artifacts and motion estimation error.

The normalized energy function plotted against the gestational age, before and after the bias field inconsistency correction. (Black diamond: before correction, white triangle: the first degree correction, white square: the second degree correction.)

For quantification, the CV of the intensity values in SP (*CV _{SP}*) and in IZ (

Figure 8(A) and (B) show *CV _{SP}* and

Before 3D reconstruction, the inspection of individual slices has shown that the proposed correction method removes visual differences in the bias field inhomogeneity. For example, the intensity profiles along the intersection lines of Fig. 2, which originally displayed significant discrepancies, are adapted to match each other more closely, as shown in Fig. 9.

After 3D reconstruction, the volumetric images were found to display less bias related artifacts and delineate structures with more subtle tissue contrast, after the bias field inconsistency correction. Figure 10 is the reconstructions of a fetal brain at GA=24.14 weeks, which demonstrated substantial decrease of CV both in SP (5.44%→4.91%) and IZ (5.33%→4.92%). Figure 10(A) and (D) are the views with a typical contrast, where the tissue boundaries at the interfaces of CP-SP and IZ-CSF are sharply resolved in axial (A) and coronal (D) planes. When displayed with greater contrast, the tissue boundary at the SP-IZ interface emerges (B,E). However, in this case, artifacts caused by the fusion of image slices with different bias fields becomes more apparent. The proposed bias field inconsistency correction method removed these artifacts, as shown in panels (C,F), where the second degree polynomial model was employed. The removal of the artifacts, consequently, improves the intensity based delineation of tissue boundaries by removing local variability of bias field strength. In panels (E,F), this improvement is apparent in the coronal view, where the boundary at the SP-IZ interface (white arrows) is better pronounced after the correction (F) than before the correction (E). Panels (G,I) are the division of uncorrected intensity values in (B,E) by the corrected image intensity values (C,F). This result demonstrates bright stripes (marked by black arrows) suggesting that the local inconsistency of the bias field in the uncorrected reconstruction is mainly caused by contributions from relatively hyperintense axial slices in this region, as indicated in the sampling density map (H,J).

Ultra-fast multi-slice imaging is increasingly important in imaging moving anatomy. When combined with retrospective slice motion correction it provides a viable route to full 3D imaging in the presence of regional tissue motion. The approach introduced in this paper contributes to this area by presenting a novel framework to estimate and remove the bias field inconsistency patterns between sets of multislice MR studies acquired from a moving region of anatomy. In such cases, the spatially varying bias field, coupled with the motion of the anatomy within the inhomogeneous bias field, induces inconsistent intensity levels for identical regions of tissue appearing in different slices. When these multislice studies are motion corrected and combined to form a single regularly sampled 3D volume, the underlying inconsistency in image intensity can introduce artifacts in the final volumetric image. The proposed method builds on earlier work to estimate relative slice motion using slice intersection matching, by deriving a least squares formulation that simultaneously resolves bias field inconsistencies between all intersecting slices.

We define an energy function to quantify the mismatch of the bias field between slice pairs. The overall energy minimization problem is then solved using a least squares approach with a set of constraint equations that conserve the spatial intensity distribution of the slices. The problem can be written into a quadratic equation, and the exact solution can be found using a Lagrange multiplier formulation.

The final solution of the equation is a set of bias field parameters for each slice that minimizes the mean square intensity difference between all slices. This solution resolves the bias field inconsistency between slices, and the residual bias field inhomogeneity can then be removed after reconstruction into a 3D volume using conventional techniques such as [24], [29]. We have shown that this postprocessing procedure successfully suppresses the resulting artifacts that appear when the scattered slice data are combined to form a regular 3D volumetric image.

The methodology presented makes use of the relative slice location information first computed using Slice Intersection Motion Correction (SIMC) method [8]. We have found that the initial motion estimation is relatively robust to the subtle changes in bias field inconsistency that can induce reconstruction artifacts, allowing the bias field inconsistency correction to be carried out as a post-hoc procedure. However, future work will examine optimal approaches to interleaving the motion correction and bias field inconsistency correction steps into a single formulation to deal with extreme cases of intensity bias.

This work was supported by NIH Grant R01 NS 055064 (PI: C. Studholme) and the data acquisition partially was funded by K23 NS52506 (PI: O.A. Glenn). The work of F. Rousseau was supported by European Research Council under FP7/2007-2013 Grant Agreement 207667.

Kio Kim, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA.

Piotr A. Habas, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA.

Vidya Rajagopalan, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA.

Julia A. Scott, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA.

James M. Corbett-Detig, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA.

Francois Rousseau, LSIIT, UMR 7005, CNRS-University of Strasbourg, Illkirch, 67412 France.

A. James Barkovich, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA.

Orit A. Glenn, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA.

Colin Studholme, Department of Radiology and Biomedical Imaging, University of California San Francisco, CA 94143 USA.

1. Levine D, Barnes P, Sher S, Semelka R, Li W, McArdle C, Worawattanakul S, Edelman R. Fetal fast MR imaging: Reproducibility, technical quality, and conspicuity of anatomy. Radiology. 1998;206:549–554. [PubMed]

2. Levine D, Hatabu H, Gaa J, Atkinson M, Edelman R. Fetal anatomy revealed with fast MR sequences. AJR. 1996;167:905–908. [PubMed]

3. Glenn O. Fetal central nervous system MR imaging. Neuroimaging Clinics of North America. 2006;16:1–17. [PubMed]

4. Prayer D, Brugger P, Prayer L. Fetal MRI: Techniques and protocols. Pediatric Radiology. 2004;34:685–693. [PubMed]

5. Rousseau F, Glenn O, Iordanova B, Rodriguez-Carranza C, Vigneron D, Barkovich A, Studholme C. A novel approach to high resolution fetal brain MR imaging. Proc Medical Image Computing and Computer Assisted Intervention (MICCAI’05) 2005:548–555. [PubMed]

6. Rousseau F, Glenn OA, Iordanova B, Rodriguez-Carranza CE, Vigneron D, Barkovich JA, Studholme C. Registration-based approach for reconstruction of high-resolution *in utero* MR brain images. Academic Radiology. 2006;13(9):1072–1081. [PubMed]

7. Kim K, Hansen M, Habas P, Rousseau F, Glenn O, Barkovich A, Studholme C. Intersection based registration of slice stacks to form 3-D images of the human fetal brain. Proc IEEE Int Symp Biomed Imag (ISBI 08) 2008:1167–1170.

8. Kim K, Habas P, Rousseau F, Glenn O, Barkovich A, Studholme C. Intersection based motion correction of multislice MRI for 3-D in utero fetal brain image formation. IEEE Trans Med Imag. 2010 Jan;29(1):146–158. [PMC free article] [PubMed]

9. Jiang S, Xue H, Glover A, Rutherford M, Rueckert D, Hajnal J. MRI of moving subjects using multislice snapshot images with volume reconstruction (SVR): Application to fetal, neonatal, and adult brain studies. IEEE Transactions on Medical Imaging. 2007 Jul;26(7):967–980. [PubMed]

10. Rousseau F, Kim K, Studholme C, Koob M, Dietmann J-L. On super-resolution for fetal brain MRI. Proc Medical Image Computing and Computer Assisted Intervention (MICCAI’10) 2010 Sep;6362:355–362. [PMC free article] [PubMed]

11. Gholipour-Baboli A, Estroff J, Warfield S. Robust super-resolution volume reconstruction from slice acquisitions: Application to fetal brain MRI. IEEE Transactions on Medical Imaging. 2010 Oct;29(10):1739–1758. [PMC free article] [PubMed]

12. Collins CM, Smith MB. Calculations of *B*_{1} distribution, SNR, and SAR for a surface coil adjacent to an anatomically-accurate human body model. Magnetic Resonance in Medicine. 2001;45:692–699. [PubMed]

13. Sled JG, Pike GB. Standing-wave and RF penetration artifacts caused by elliptic geometry: An electrodynamic analysis of MRI. IEEE Transactions on Medical Imaging. 1998 Aug;17(4):653–662. [PubMed]

14. Levine D, Smith AS, McKenzie C. Tips and tricks of fetal MR imaging. Radiologic Clinics of North America. 2003;41(4):729–745. [PubMed]

15. Arnold J, Liow J-S, Schaper K, Stern J, Sled J, Shattuck D, Worth A, Cohen M, Leahy R, Mazziotta J, Rottenberg D. Qualitative and quantitative evaluation of six algorithms for correcting intensity nonuniformity effects. NeuroImage. 2001;13:931–943. [PubMed]

16. Vovk U, Pernus F, Likar B. A review of methods for correction of intensity inhomogeneity in MRI. IEEE Transactions on Medical Imaging. 2007 Mar;26(3):405–421. [PubMed]

17. Belaroussi B, Milles J, Carme S, Zhu YM, Benoit-Cattin H. Intensity non-uniformity correction in MRI: Existing methods and their validation. Medical Image Analysis. 2006;10:234–246. [PubMed]

18. Lai S-H, Fang M. A dual image approach for bias field correction in magnetic resonance imaging. Magn Reson Imaging. 2003;21:121–5. [PubMed]

19. Fan A, Wells WM, Fisher JW, Cetin M, Haker S, Mulkern R, Tempany C, Willsky AS. A unified variational approach to denoising and bias correction in MR. Inf Process Med Imaging. 2003 Jul;18:148–159. [PubMed]

20. Learned-Miller EG, Ahammad P. Joint MRI bias removal using entropy minimization across images. In: Saul LK, Weiss Y, Bottou L, editors. Advances in Neural Information Processing Systems. Vol. 17. Cambridge, MA: MIT Press; 2005. pp. 761–768.

21. Ashburner J, Friston K. Unified segmentation. Neuroimage. 2005;26:839–851. [PubMed]

22. Fletcher R. Practical methods of optimization. John Wiley and Sons; 1987.

23. Kim K, Habas P, Rajagopalan V, Scott J, Rousseau F, Barkovic A, Glenn O, Studholme C. SLIMMER: Slice MRI Motion Estimation and Reconstruction tool for studies of fetal anatomy. Medical Imaging 2010: Image Processing, Proc. SPIE; 2011.

24. van Leemput K, Maes F, Vandermeulen D, Suetens P. Automated model-based bias field correction of MR images of the brain. IEEE Trans Med Imag. 1999 Oct;18(10):885– 896. [PubMed]

25. Habas P, Kim K, Corbett-Detig J, Rousseau F, Glenn O, Barkovich A, Studholme C. A spatiotemporal atlas of MR intensity, tissue probability and shape of the fetal brain with application to segmentation. Neuroimage. 2010;53(2):460–470. [PMC free article] [PubMed]

26. Shepp LA, Logan BF. The Fourier reconstruction of a head section. IEEE Transactions on Nuclear Science. 1974;21(3):21–43.

27. Habas P, Kim K, Rousseau F, Glenn O, Barkovich A, Studholme C. Atlas-based segmentation of developing tissues in the human brain with quantitative validation in young fetuses. Hum Brain Mapp. 2010 Sep;31(9):1348–1358. [PMC free article] [PubMed]

28. Corbett-Detig JM, Habas PA, Scott JA, Kim K, Rajagopalan V, McQuillen PS, Barkovich AJ, Glenn OA, Studholme C. 3D global and regional patterns of human fetal subplate growth determined in utero. Brain Structure and Function. 2011;215:255–263. [PMC free article] [PubMed]

29. Sled J, Zijdenbos A, Evans A. A nonparametric method for automatic correction of intensity nonuniformity in MRI data. IEEE Transactions on Medical Imaging. 1998;17(1):87–97. 2. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |