Anita Möller^{1}, Marco Maass^{1}, Tim Jeldrik Parbs^{1}, and Alfred Mertins^{1}

A blind retrospective MRI motion estimation and compensation algorithm is designed for arbitrary sampling trajectories. Using the idea of natural images being sparsely representable, the algorithm is based on motion estimation between a motion corrupted image and it’s sparse representative. Therefore, rigid motion models are designed and used in gradient descent methods for image quality optimization. As the motion estimation and compensation work on arbitrary real valued sampling coordinates, the algorithm is capable for all trajectories. Image reconstruction is performed by computationally efficient gridding. The exact motion estimation results are shown for PROPELLER and radial trajectory simulation.

One MRI k-space measurement is composed of several partial measurements $$$y_n$$$ along specified trajectories within one readout per different time intervals $$$n$$$. Patient motion within one readout is assumed to be constant due to small intervals. A suitable image reconstruction problem is given by $$\hat{\boldsymbol{x}}=\underset{\boldsymbol{x}}{\arg\min}\sum_{n=1}^N\left\|\mathcal{S}_n\mathcal{F}\mathcal{D}_n\boldsymbol{x}-\boldsymbol{y}_n\right\|_2^2+\frac{1}{\sigma^2}\Phi(\boldsymbol{x})$$ with the sampling operator $$$\mathcal{S}_n$$$ representing the trajectory at time $$$n$$$, $$$\mathcal{F}$$$ denoting the Fourier transform of the image $$$\boldsymbol{x}$$$ which is motion corrupted by operator $$$\mathcal{D}_n$$$.

Sampling is
described by $$$\mathcal{S}\mathcal{F}\boldsymbol{x}$$$ on arbitrary k-space
trajectories and needs the computation with real-valued nonequidistant
frequency coordinates. It is calculated by the nonequidistant discrete Foruier
transform (NDFT)^{2}.

Two motion
models are combined in $$$\mathcal{D}_n$$$ for translation and rotation
separately. Full and sub pixel translation shifts are mathematically modelled
by convolution matrices in image space^{3}. Rotation is described by
applying a rotation matrix to the sampling trajectory coordinates followed by a
barycentric interpolation^{4} combined with Delaunay triangulation^{5}
of the coordinates. This delivers a differentiable motion model on arbitrary nonequidistant
frequency coordinates.

A meaningful regularization on motion corrupted MR images is enforcing sparsity in the wavelet domain by $$$\Phi=\|W\boldsymbol{x}\|_1$$$ to give advantage to clear natural structures.

For motion
estimation, a three-step iteration algorithm is performed. A sparsifying ADMM^{6}
is evaluated in the first step. It solves the minimization for fixed motion
parameters by splitting it into a measurement fitting and a sparsifying problem
which are iteratively solved and updated by
$$\hat{\boldsymbol{x}}_{d+1}=
\underset{\boldsymbol{x}}{\arg\min}\sum_{n=1}^N\|\mathcal{S}_n\mathcal{F}\mathcal{T}_n\boldsymbol{x}-\boldsymbol{y}_n\|_2^2+\frac{1}{\lambda^2}\left\|\boldsymbol{x}-\bar{\boldsymbol{x}}_{d+1}\right\|_2^2,
\quad\bar{\boldsymbol{x}}_{d+1}=\hat{\boldsymbol{v}}_d-\boldsymbol{u}_d,$$
$$\hat{\boldsymbol{v}}_{d+1}=
\underset{\boldsymbol{v}}{\arg\min}\left\|W\boldsymbol{v}\right\|_1+\frac{\sigma^2}{\lambda^2}\left\|\boldsymbol{v}-\bar{\boldsymbol{v}}_{d+1}\right\|^2_2,\quad
\bar{\boldsymbol{v}}_{d+1}=\hat{\boldsymbol{x}}_{d+1}+\boldsymbol{u}_d,$$
$$\boldsymbol{u}_{d+1}=\boldsymbol{u}_d+\hat{\boldsymbol{x}}_{d+1}-\hat{\boldsymbol{v}}_{d+1}.$$
The
measurement fitting problem is solved by conjugate gradients^{7} using
the inverse NDFT^{2}. The sparsity of the solution for the second
problem is enforced by soft thresholding in the wavelet domain. In the second
step, motion estimation is achieved by using the sparsified reconstruction of
the ADMM in the derivatives of the motion models with Newton’s gradient method^{7}.
Thereby, the parametrization of translation is improved with a Gaussian model.
The last step updates the globally estimated motion on the k-space frequency
coordinates and coefficients.

Finally, image reconstruction is performed by
gridding^{8} to overcome blurring of the inverse NDFT. The frequency
coefficients are resampled onto a Cartesian grid. As the arbitrary trajectory
coordinates were not equally dense over the whole k-space, the coefficients are
convolved with an area density compensation function. This weighting is
computationally efficient calculated by partitioning the k-space into
rectangles around the sampling points and saving the structure
in k-d trees^{9}.

Representatively,
PROPELLER and radial trajectories were simulated. The algorithm was tested on
BrainWeb^{10} and Shepp-Logan phantom data in a FOV of 455 and 160
pixels, respectively. Smooth patient motion was modeled as an autoregressive
moving average process with maximum translation amplitudes up to 30 pixels and maximum
rotation up to 45°.

Figure 1 shows motion corrupted images and the motion compensated results of the algorithm calculated with PROPELLER trajectories. In the corrupted images, no structure or even no image is identifiable. The reconstructed images do not show motion artefacts, but clear contours and details are visible. Only a few gridding artefacts appear due to low resolution. In Figure 2, the progresses of the corresponding originally applied and by the algorithm estimated motion are shown. Their courses follow very closely. The table in Figure 3 shows the mean percentage improvement of the PSNR and mutual information (MI) between motion corrupted images and motion compensated reconstructions for different maximum translation and rotation amplitudes.

The results prove the ability of the algorithm to estimate rigid motion very exactly. Even from images corrupted with motion shifts about one fifth of the FOV very detailed images are reconstructed. Small differences between the original and estimated motion just shift the position of the image centroid but do not effect the image quality. The high improvements of the image quality measures emphasize the impression given by the reconstructed images. The method is easily adaptable to other motion models and with it expandable to physical descriptions of elastic motion. As sampling is mathematically described for real-valued nonequispaced frequency coordinates the algorithm quality is largely independent of the design of MRI trajectories.

1. J.G. Pipe. Motion Correction With PROPELLER MRI: Application to Head Motion and Free-Breathing Cardiac Imaging. Magn Reson Med, 42(5):963-969, 1999.

2. J. Keiner, S. Kunis, D. Potts. Using NFFT 3-A Software Library for Various Nonequispaced Fast Fourier Transforms. ACM Trans Math Softw, 36(4):19:1-19:30, 2009.

3. A. Möller, M. Maass, A. Mertins. Blind Sparse Motion MRI with Linear Subpixel Interpolation. Proc BVM, 510-515, 2015.

4. K. Hormann. Barycentric Interpolation. In: G.E. Fasshauer, editor. Approximation Theory XIV: San Antonio 2013. Springer, Cham, 197-218, 2014.

5. F. Aurenhammer, R. Klein, D.T. Lee. Voronoi Diagrams and Delaunay Triangulations. WORLD Scientific, 2013.

6. N. Parikh, S. Boyd, Proximal Algorithms. Found Trends Optim, 1(3):127-239, 2014.

7. J. Nocedal, S.J. Wright. Numerical Optimization-2nd ed. Springer, New York, 2006.

8. O.K. Johnson, J.G. Pipe. Convolution Kernel Design and Efficient Algorithm for Sampling Density Correction. Magn Reson Med, 61(2):439-447, 2009.

9. J.L. Bentley. Multidimensional Binary Search Trees Used for Associative Searching. Commun ACM, 18(9):509-517, 1975.

10. C.A. Cocosco, V. Kollokian, R.K.S. Kwan, E.C. Evans. BrainWeb: Online Interface to a 3D MRI Simulated Brain Database. Neuroimage, 5(4):425, 1997.

Figure 1: Top line: Motion corrupted images of the BrainWeb and Shepp-Logan phantom. Bottom line: Corresponding motion compensated reconstructions. Original motion with maximum shift of 30 pixels per dimension and maximum rotation of 45° was applied. It is correspondingly shown with the estimated motion used for motion compensation in Figure 2.

Figure 2: Original and estimated absolute motion for the corresponding images in Figure 1 per measurement time $$$n$$$.

Figure 3: Mean percentage improvement of the image quality measures MI and PSNR for different images and maximum translation and rotation amplitudes over five runs each.