## 2 Theory  
### Joint Reconstruction of an Image and a Motion Model  
A review of motion-compensated reconstruction methods can be found in [11]. In this work, motion-compensated reconstruction was accomplished using GRICS [12]. The algorithm solves for both a motion-corrected image \(\rho\) and the parameters of a motion model \(\alpha\). The GRICS problem is mathematically formulated as follows:  
\[ \begin{cases}\rho,\alpha=\operatorname*{arg\,min}_{\rho,\alpha}||E(u_{\alpha })\rho-m||^{2}+\lambda||\rho||^{2}+\mu R(\alpha),\\ u_{\alpha}(x,y,z,t)=\begin{cases}[x,y,z,1]\alpha(t)\text{ (rigid model), or}\\ \alpha(x,y,z)s(t)\text{ (non-rigid model)}\end{cases}\end{cases}  \tag{1}\]  
\(m\) stands for the acquired MRI \(k\)-space data; \(m\) is a vector of size \(N_{m}\). \(\rho\) is the motion-corrected image, a vector of \(N_{\rho}\) elements. \(u\) is the displacement field at each \(k\)-space sample time; it is a vector of \(N_{m}N_{\rho}N_{d}\) elements, where \(N_{d}\) is the number of dimensions in the image (\(N_{d}\in\{2,3\}\)). In the case of a rigid motion model (e.g., for cerebral imaging), \(\alpha(t)\) is an affine transformation matrix using homogeneous coordinates. In the case of a non-rigid motion model (e.g., for chest or abdominal imaging), a separable motion model is used, with a spatial component \(\alpha(x,y,z)\), and the temporal component \(s(t)\) is a motion signal, acquired from a motion sensor (here, a pneumatic respiratory belt). A regularization term is added in the non-rigid case as \(R(\alpha)=||\
abla\alpha||^{2}\). \(E\) is an encoding operator; it is expressed using the following operators:\[E=\left[ \begin{array}{c}S_{t_{1}}FC_{1}T(u_{t_{1}})\\ \vdots\\ S_{t_{N}}FC_{1}T(u_{t_{N}})\\ \vdots\\ S_{t_{1}}FC_{N_{c}}T(u_{t_{1}})\\ \vdots\\ S_{t_{N}}FC_{N_{c}}T(u_{t_{N}})\end{array} \right]\]  
The acquisition is subdivided into \(N_{t}\) motion states (shots) and \(N_{c}\) receiver coils. For a given \(i\leq N_{t}\), each motion state \(t_{i}\) is described by a dense displacement vector field \(u_{t_{i}}\) and a sparse spatial transformation matrix \(T(u_{t_{i}})\). The size of \(T(u_{t_{i}})\) is \(N_{p}\times N_{p}\). The matrices \(S_{t_{i}}\) describe the sampling pattern of the \(k\)-space data during the \(i\)th motion state. In the Cartesian sampling case, \(S_{t_{i}}\) represents sparse matrices, having 0 and 1 as elements. The matrix \(S_{t_{i}}\) has \(N_{p}\) columns, and its number of rows is the number of \(k\)-space samples acquired at each shot \(t_{i}\). \(F\) stands for the Fourier transform operator. \(C_{1},C_{2},\ldots,C_{N_{c}}\) represent the diagonal coil sensitivity operator of the receiver coils. For each \(j\leq N_{c}\), the size of the matrices \(C_{j}\) and operator \(F\) is \(N_{p}\times N_{p}\). In GRICS, the Fast Fourier Transform of a vector is computed using the FFTW implementation [24]. The two least squares problems in Algorithm 1 are solved using the conjugate gradient solver. The Gauss-Newton scheme is introduced to linearize the cost function, line 8 around the current estimate of \(\alpha\); each iteration suggests a refinement \(\delta\alpha\) of \(\alpha\). The tolerances are \(\epsilon_{r}=10^{-3}\) for the reconstruction solver and \(\epsilon_{m}=10^{-2}\) for the motion problem. In the last resolution level, only the reconstruction step is performed.