4496
Blind Sparsity Based Motion Estimation and Correction Model for Arbitrary MRI Sampling Trajectories1Institute for Signal Processing, Universität zu Lübeck, Lübeck, Germany
Synopsis
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.
Introduction
Several motion compensation methods exist for each sampling trajectory, measuring modality, and imaging object, separately. Often they are based on extra navigator measurements, reduction of patient comfort or extension of measurement time like in gating. But at all, free non-periodic patient motion like swallowing is still not compensatable. There is a lack of a general blind motion estimation algorithm retrospectively compensating free patient motion that is applicable to arbitrary sampling trajectories. A model for such an algorithm is designed based on the sparsity of MR images, gradient descent and fast gridding. Simulations were driven on PROPELLER1 and radial trajectories.Methods
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 space3. Rotation is described by applying a rotation matrix to the sampling trajectory coordinates followed by a barycentric interpolation4 combined with Delaunay triangulation5 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 ADMM6 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 gradients7 using the inverse NDFT2. 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 method7. 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 gridding8 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 trees9.
Results
Representatively, PROPELLER and radial trajectories were simulated. The algorithm was tested on BrainWeb10 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.
Discussion
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.
Acknowledgements
This work has been supported by the German Research Foundation under Grant No. ME 1170/11-1.References
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.
Figures
