Introduction

Retrospective cardio-respiratory motion correction during dynamic MRI scans relies on grouping the acquired data into motion states or gates. The data from each gate are then reconstructed and the resulting images co-registered to a “reference” gate to produce an image incorporating all the acquired data [1]. This approach is of limited effectiveness when highly accelerated scans are used due to spatial transformation errors introduced by under sampling artifacts. We demonstrate that by enforcing data consistency on the transformation fields ([2], i.e., transformed images must match k-space data) the performance of retrospective motion correction can be greatly improved.

Methods

Joint reconstruction of a single motion-corrected MR image $$$y$$$ and motion fields $$$S_k$$$ using $$$l$$$ gates can be formulated as the solution of:

$$\underset{S_k,y}{\text{argmin}} \ \alpha \, R_1(S) + \sigma \, R_2(y) + \sum_{k=1}^l D_k(S_k y) $$

where $$$D_k(\cdot)$$$ is the $$$L_2$$$ differences of the acquired MR k-space data from gate $$$k$$$ and the MR forward model of the transformed image $$$y$$$ using the motion fields $$$S_k$$$. $$$R_1$$$ is a motion field's regularizer $$$S_k$$$ and $$$R_2$$$ a (total variation) penalty on the image $$$y$$$. This optimization problem can be reformulated as a constrained minimization problem
$$S,y = \underset{\underset{s.t. z_k=S_k y}{y,z_k,S_k, }}{\text{argmin}} R_1(S) + R_2(y) + \sum_k D_k(z_k)$$
which has the advantage of decoupling $$$S_k$$$ and $$$y$$$ in the data fidelity terms. This problem can be solved using the augmented Lagrangian
$$\mathcal{L}_{\rho}(z_k,y,S_k,u_k)=R_1(S)+R_2(y)+ \sum_k \Bigl( D_k(z_k) \\ + \frac{\rho}{2} ||z_k - S_k y +u_k||_2^2 -\frac{\rho}{2}||u_k||_2^2 \Bigr) (1)$$
where $$$u_k$$$ are the scaled dual variables and $$$\rho$$$ is an internal penalty parameter of $$$\mathcal{L}_{\rho}$$$. Minimizing Eq. 1 using the well-known Alternating Direction Method of Multipliers (ADMM) [3], yields an iterative scheme consisting of three sub-problems and a dual variable update
$$z_k^{(n+1)} = \underset{z_k}{\text{argmin}} \, \mathcal{L}_{\rho} (z_k,y^{(n)},S_k^{(n)},u_k^{(n)} )$$
$$y^{(n+1)} = \underset{y}{\text{argmin}} \, \mathcal{L}_{\rho} ( z_k^{(n+1)},y,S_k^{(n)},u_k^{(n)} ) $$
$$S_k^{(n+1)} = \underset{S_k}{\text{argmin}} \, \mathcal{L}_{\rho} ( z_k^{(n+1)},y^{(n+1)},S_k,u_k^{(n)} )$$
$$u_k^{(n+1)} = u_k^{(n)} + z_k^{(n+1)} - S_k^{(n+1)}y^{(n+1)}$$

The first sub-problem (data-fidelity plus quadratic penalty),
$$z_k^{(n+1)} = \underset{z_k}{\text{argmin}} \, D(z_k) + \frac{\rho}{2}||z_k - S_k^{(n)}y^{(n)} + u_k^{(n)} ||_2^2$$
is solved by applying the conjugate gradient method. Solution of the second sub-problem,
$$y^{(n+1)} = \underset{y}{\text{argmin}} \, \sigma ||\nabla_3 y||_{2,1} +\frac{\rho}{2} \sum_k ||S_k^{(n)} y - z_k^{(n+1)} -u_k^{(n)} ||_2^2$$
is accelerated using the approximation
$$y^{(n+1)} \approx \underset{y}{\text{argmin}} \, \sigma ||\nabla_3 y||_{2,1} +\frac{\rho}{2} || y - v ||_2^2 (2)$$
with $$$v=\frac{1}{l}\sum_k S_k^{(n)^{-1}}\left(z_k^{(n+1)}+u_k^{(n)}\right)$$$, with $$$S_k^{(n)^{-1}}$$$ being the inverse of motion field $$$S_k^{(n)}$$$. Eq. (2) is a simple denoising-like problem solved using the first-order primal-dual hybrid gradient algorithm. The third sub-problem
$$S_k^{(n+1)}=\underset{S}{\text{argmin}} \, \sum_k \left(||S_k y - z_k^{(n+1)} -u_k^{(n)}||_2^2 \right) $$
is solved using the diffeomorphic demons algorithm with Gaussian smoothing [4].

The proposed joint reconstruction algorithm was evaluated using simulated data with known ground truth motion fields derived from the XCAT phantom [5] as well as experimental data. K-space data were generated from XCAT images with a voxel size of $$$1\,mm^3$$$ for 5 different respiratory gates with a varying number of golden angle k-space "spokes'' distributed randomly across the gates, mimicking a free-breathing radial stack-of-stars acquisition. The reconstructions in the proposed method were performed using a voxel size of 4.5 x 2 x 2 $$$mm^3$$$ (80 x 256 x 256 voxels). Data sets of varying under-sampling levels were created: 400 (full sampling), 200, 100, 50, and 25 spokes in-plane per gate, respectively. Joint reconstruction was applied to all data sets using $$$\sigma=(0.01, 0.1, 1.0)$$$ for regularization and $$$(3,5,7,9)$$$ for demons' smoothing. Iterative, TV-regularized reconstructions using the same values for $$$\sigma$$$ were also performed, followed by image-based spatial co-registration (demons) of the resulting images using the same levels of smoothing as above. Pixel-wise spatial transformation errors were calculated for all methods using the ground-truth simulated values.

Results

Figure 1 shows calculated spatial deformation errors for reconstructions using different levels of acceleration. Spatial deformation errors increase significantly for the image-based method with increased acceleration, especially towards the edges of the abdomen (figure 2). Images corresponding to the same reconstruction conditions for both methods are presented in figure 3 and demonstrate significant spatial blur for image-based motion correction at high acceleration levels. By contrast, the behavior of the data-consistent approach remains stable under the same conditions. The same behavior is observed when experimental data from a healthy human control is used (figure 4). Notably, most of the image-based motion correction blur appears in areas with significant artifact levels (arrows). These errors are not present in the data-consistent approach.

Conclusions

We have shown a data-consistent spatial deformation approach for motion correction of dynamic MRI scans. This approach exhibits significant improvements in performance relative to the the traditional image-based spatial co-registration approach.

Acknowledgements

Supported in part by PHS grant R01 EB029306.