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

*n* is the number of voxels in the 3D image. A subregion of interest of

*z*_{rigid} undergoes rigid motion with respect to the imaging system. We consider this region

*z*_{rigid} ^{n} to be given by

*z*_{rigid} =

*z*_{whole}, where

is a masking operator which selects those points moving collectively as a rigid object from those deforming around it.

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;

where

is an

*n × n* diagonal matrix.

When a 2D slice image consisting of

*m*_{i} pixels,

*y*_{i} ^{mi} is acquired from a given slice selected by an operation S

_{i} with the spatial transformation of the rigid anatomy at that time taking the form T

_{i}, the 2D slice can be expressed as

where R

_{i} is the masking operator

in the

*i*-th slice,

an

*n × n* diagonal matrix for the 3D bias field distortion of the

*i*-th slice acquisition, D

_{i} an

*m*_{i} × n unitary operator which down-samples the 3D volume into the

*i*-th 2D image slice and,

*η* the additive noise. This matrix representation is explained in .

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

where superscript + denotes the Moore-Penrose pseudo inverse.

Matrices D

_{i} and S

_{i} can be directly obtained from the scanner setup and T

_{i} can be estimated using motion estimation algorithms [

6], [

8], [

9], The value of the individual slice intensity bias

, however, remains unknown, because motion of the coil array and the subject anatomy during image acquisition is unavailable. Since the undetermined

can corrupt the final reconstruction

_{rigid}, we eliminate

in (

3) by estimating the unbiased 2D slice image

;

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

where the 2D bias field operator B

_{i} is defined by

, which is an

*m*_{i} ×

*m*_{i} diagonal matrix. Although the true bias field

is not available, we can estimate B

_{i} that minimizes the intensity difference of the shared anatomy in the slices.

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:

where the slice selection operators

*S*_{i} and

*S*_{j} are cross-operated to extract the intersections of slices

*i* and

*j*.

For the

*i*-th image slice, we define a 2D windowing vector

*r*_{i} = diag{R

_{i}}, and a 2D bias field vector

*b*_{i} = diag{B

_{i}}, then (

5) is rewritten as

for the voxel located at

*x*_{i}. We further define a bias field correction term

as a functional of the parameter vector

*θ*_{i}, namely,

.

Eq. (6) is then rewritten in a functional form;

where Ω

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

*i*-th slice, and

*ω* is a point along the intersection line between slice

*i* and slice

*j* in 3D space. We estimate the bias field parameter vector

*θ*_{i} by simultaneously comparing all intersecting slice pairs along their intersection lines.

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

, the energy function is

where

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

*i*-th slice,

*r*(

*ω*) the 3D windowing function for the rigid object at

*ω*, and

*N* the total number of slices in the study. We define Θ, a vector representation of the bias field parameters of all slices, namely,

. This energy function assumes the noise is added before the signal scaling, to model the noise due to the tissue property variation. [

21]

Given that the bias field correction term

is a linear functional of the bias field parameter vector

*θ*_{i} and the 2D spatial location

*x*;

where

*u*(

*x*) is a vector of spatial basis functions at the location

*x*. The basis functions of

*u* can be any linear combination of linearly independent functions. For example, one can choose

for the first degree bias correction model, at the voxel location

*x* = (

*x*_{1}*, x*_{2}) in the

*i*-th slice, and

for the second degree model. The energy function in (

9) then can be rewritten in a quadratic form

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

For a

*d*-th degree model,

**M** is a square matrix with

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;

where

*x*_{i}(

*ω*) corresponds to the voxel location in the

*i*-th slice image and

*v*(

*ω*) the voxel location in the 3D reference space.

*k*(

*v*) is a polynomial vector, which is set to

for the first degree model, and

for the second degree model, etc. We define the intensity distribution matrix

, where

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

where Θ

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

The solution

of the constrained quadratic equation

can be found using the method of Lagrange multipliers [

22] where the Lagrangian function is given by

and

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

24) is

Note that (

25) is the exact and general solution to

Eqs. (22)–

(23), for any multiplicative bias field model.