Motion Correction in MRI

07/03/2024

Recently I read some articles on motion correction. I listed them as follows:

Daniel Polak et al. Scout accelerated motion estimation and reduction (SAMER). MRM (2021)
Niek RF Huttinga at al. MR-MOTUS: model-based non-rigid motion estimation for MR-guided radiotherapy using a reference image and minimal k-space data. PMB (2020).
Veronike Spieker et al.Deep Learning for Retrospective Motion Correction in MRI: A Comprehensive Review. TMI (2024)

Not very familiar with motion correction, I would like to start by writing some general thoughts on this subject. First of all, motion can be divided into two categories: rigid motion and non-rigid motion. Regarding rigid motion, most MRI papers focused on brain motions; Regarding non-rigid motion, the most mentioned motions are those from the heart and lungs.

Matrix description of rigid motion

For any rigid/translational/affine motion, We can build a linear relationship between the motion-free/corrected k-space $v (k)$ and motion-corrupted k-space $v^{'} (k^{'})$ .

v (k) = v^{'} (k^{'}) \frac{e x p (2 π i k^{'} \cdot t)}{| d e t (A) |}

where $k^{'}$ is the motion-corrupted k-space location, and $k$ is the motion-free k-space location. $t = [t_{x}, t_{y}, t_{z}]^{T}$ is the translation in three directions. and A is an affine motion matrix, which can be written as follows:

A = R G S = R_{x} R_{y} R_{z} G S

= [\begin{matrix} 1 & 0 & 0 \\ 0 & c o s ϕ & - s i n ϕ \\ 0 & s i n ϕ & - s i n ϕ \end{matrix}] [\begin{matrix} c o s θ & 0 & s i n θ \\ 0 & 1 & 0 \\ - s i n θ & 0 & c o s θ \end{matrix}] [\begin{matrix} c o s ψ & - s i n ψ & 0 \\ s i n ψ & c o s ψ & 0 \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} 1 & G_{x y} & G_{x z} \\ 0 & 1 & G_{y z} \\ 0 & 0 & 1 \end{matrix}] [\begin{matrix} S_{x} & 0 & 0 \\ 0 & S_{y} & 0 \\ 0 & 0 & S_{z} \end{matrix}]

where $R_{x / y / z}$ are rotation matrix with angle $ϕ, θ, ψ$ . matrix $G$ is the shearing matrix, and S is the scaling matrix.

Affine motion

Applied to MRI, in the image domain, the motion corrupting process can be described as follows:

[\begin{matrix} x^{'} \\ y^{'} \\ z^{'} \\ 1 \end{matrix}] = [\begin{matrix} t_{x} \\ R & t_{y} \\ t_{z} \\ 0 & 0 & 0 & 1 \end{matrix}] G S [\begin{matrix} x \\ y \\ z \\ 1 \end{matrix}]

Rigid motion

$G = S = I$ , where $I$ is the identity matrix.

Translational motion

For this kind of motion, we will have $R = G = S = I$

Non-rigid motion

The dominant approach for non-rigid motion compensation in MRI is to estimate the motion field. For this purpose, estimation can be done from low-resolution images or directly from k-space data.

What's the motion field?

[Wikipedia]:
A camera captures 2D images from a 3D world. In other words, it maps each points $(x_{1}, x_{2}, x_{3})$ in 3D space to a point $(y_{1}, y_{2})$ in 2D image.

(\begin{matrix} y_{1} \\ y_{2} \end{matrix}) = (\begin{matrix} m_{1} (x_{1}, x_{2}, x_{3}) \\ m_{2} (x_{1}, x_{2}, x_{3}) \end{matrix})

When it comes to motion, we differentiated the above expression concerning time as follows:

(\begin{matrix} \frac{d y_{1}}{d t} \\ \frac{d y_{2}}{d t} \end{matrix}) = (\begin{matrix} \frac{d m_{1} (x_{1}, x_{2}, x_{3})}{d t} \\ \frac{d m_{2} (x_{1}, x_{2}, x_{3})}{d t} \end{matrix}) = (\begin{matrix} \frac{d m_{1}}{d x_{1}} & \frac{d m_{1}}{d x_{2}} & \frac{d m_{1}}{d x_{3}} \\ \frac{d m_{2}}{d x_{1}} & \frac{d m_{2}}{d x_{2}} & \frac{d m_{2}}{d x_{3}} \end{matrix}) (\begin{matrix} \frac{d x_{1}}{d t} \\ \frac{d x_{2}}{d t} \\ \frac{d x_{3}}{d t} \end{matrix})

Here, we define the motion field $u$ as:

u = (\begin{matrix} \frac{d y_{1}}{d t} \\ \frac{d y_{2}}{d t} \end{matrix})

However, the problem is in practice we only can estimate the motion field based on image data. This will lead to our next topic-- how to estimate the field map?

Optical flow

...

Free-form deformation

As we discuss in rigid motion, The 3D affine transformation is defined as:

T^{a f f} (r | A, t) = A r + t

The free-form deformation is defined as:

T^{F F D} (r | c^{x}, c^{y}, c^{z}) = r + (\begin{matrix} b^{x} (r) c^{x} \\ b^{y} (r) c^{y} \\ b^{z} (r) c^{z} \end{matrix})

where $b^{p} (r) \in ℝ^{1 \times N_{c}^{p}}$ are the row-vectors with $N_{c}^{p}$ 3D B-spline basis functions, evaluated at the coordinate $r$ , and $c^{p} \in ℝ^{N_{c}^{p} \times 1}$ denotes expansion coefficients. The 3D basis functions are constructed as a Kronecker product of three components of the motion-field.

How to estimate the motion field

Generally, you can estimate motion field from both $k - s p a c e$ domain and image domain (image registration). For accelerated cases, k-space-based motion field estimation will be a fit, as the undersampling-reduced artifacts will degrade the estimated motion field.

Here, we will take the MR-MOTUS, a k-space-based method, as the example. The article gives the basic signal equation at first:

s_{t} (k) = \int_{Ω} q_{t} (r) e^{- i 2 π k \cdot r} d r

where $q_{t} (r)$ denote the transverse magnetization of a deforming object at time $t$ and spatial coordinate $r = (x, y, z)$ . The $k$ -space signal from $q_{t}$ at coordinate $k = (k_{x}, k_{y}, k_{z})$ .

The motion can be described as:

U_{t} (r) = r + δ_{t} (r)

where $U_{t} : ℝ^{3} \mapsto ℝ^{3}$ denote the motion-field that deforms $q_{0}$ to $q_{t}$ . $δ_{t} : ℝ^{3} \mapsto ℝ^{3}$ denote disolacement function.

According to the local spin conservation: $ρ (r_{t}) d r_{t} = ρ_{0} (U_{t} (r_{t})) | d e t (▽ U_{t} (r)) | d r$ , the basic signal equation can be rewritten as follows:

s_{t} (k) = \int_{Ω} q_{0} (r_{0}) e^{- i 2 π k \cdot T_{t} (r_{0})} d r_{0}

Note that $T_{t}$ is defined as the inverse of $U_{t}$ , such that:

T_{t} (r) = r + η (r), U_{t} \cdot T_{t} = I,

If you want to know more details about the proof, please go to MR-MOTUS paper.

Currently, the $k$ -space signal is related to a reference object (assumed it is known) through non-rigid/non-linear motion-fields $T_{t}$ . In addition, the authors represent the motion-field in a low-dimensional basis using coefficients $θ_{t} \in ℝ^{N_{c}}$ , where $N_{c} ≪ 3 N$ , and N is the number of voxels per motion-field. Now we have:

s_{t} = F (θ_{t} | q_{0}) [k] = \int_{Ω} q_{0} (r_{0}) e^{- i 2 π k \cdot T_{t} (r_{0} | θ_{t})} d r_{0}

In order to better reconstruct the motion-field, a regularization term was added. The general form of motion-field estimation can be rewritten as follows:

min_{θ_{t}} ‖ F (θ_{t}) - s_{t} ‖_{2}^{2} + λ R (T_{t} (\cdot | θ_{t}))

The regularization term is designed to smooth the motion-field by penalizing the spatial curvature of the motion-fields.

R (T_{t} (\cdot | θ_{t})) = \sum_{p \in x, y, z} \int_{Ω} | △ T_{t}^{p} (r | θ_{t}^{p}) |^{2} d r

where $△$ denotes the Laplace operator, and $T_{t}^{x, y, z} : R^{3} \mapsto ℝ$ denote the individual components of the motion-field $T_{t}$ . The final objective function will be:

min_{θ_{t}} ‖ F (θ_{t}) - s_{t} ‖_{2}^{2} + λ \sum_{p \in x, y, z} \int_{Ω} | △ T_{t}^{p} (r | θ_{t}^{p}) |^{2} d r

How to correct/compensate?

[To be continued]

#MRI