In a conventional MR-MOTUS framework, we define target time-series images\(H = [h_ ,...,h_ ]\) a motion static reference image in its vectorized form \(}\) of length \(N\) (total number of voxels) and motion fields \(D = [D_ ,...,D_ ]\) for each dynamic index \(t\), up to \(M\) total number of dynamics. These quantities are mutually related as follows,

$$h_ (r) = q(D_ (r)) \cdot \det (\nabla D_ )$$

(1)

where \(}=\left(x,y,z\right)\) are spatial coordinates and \(\nabla\) is the Jacobian operator. From Eq. (1), \(}}_\) is approximated by warping the reference image \(}\), which is the fundamental requirement for the MR-MOTUS framework. This approximation is incorporated into the MR signal model leading to the following MR-MOTUS signal model \(}}\) [18].

$$S_ (k) = \hat(D_ |q) = \int_ \,e^^ (r)}} dr$$

(2)

The notation \(\hat(D_ |q)\) denotes that the model is conditional to the specific value of \(}\). Here, the reference image \(}\) is related to the motion fields \(}\) at time point \(t\) for the measured k-space data \(s_\), and \(}=\left(_,_,_\right)\) are k-space coordinates.

One limitation in the previous MR-MOTUS implementations, the coil dimension of the k-space has been compressed into a single homogeneous coil to reduce computation time. In this paper, we assume that the coil maps are static throughout the acquisition. Meaning that a single reference image needs to be warped before coil maps are applied.

Another limitation is that it does not consider contrast variation in the reference image, which can occur in cardiac examinations due to blood flow or contrast administration. Additionally, \(}\) must be motion static, meaning that breath-hold or gated images need to be acquired in a preparation phase before the examination to use the MR-MOTUS framework. To tackle these challenges, we made a small but crucial extension to the MR-MOTUS framework where the reference image is allowed to change over time, resulting in a time series of images \(Q = [q_ ,...,q_ ]\). Although this modification enhances the validity and flexibility of the MR-MOTUS model, it also drastically augments the number of unknowns in the reconstruction problem, making it even more ill-posed and challenging to solve for the motion fields. To implement this, we use the low-rank plus sparse (L + S) decomposition [12] as a regularizer for \(}\). Additionally, there is no longer a need for a preparation phase to acquire a motion static reference image, instead contrast-varying motion static reference images are jointly estimated during alternating iterations using the entire time-resolved data. We call this joint approach CMR-MOTUS.

According to this approach, H is approximated by

$$h_ (r) = q_ (D_ (r)) \cdot \det (\nabla D_ )$$

(3)

Hence, a slightly modified signal model F is derived to include this time-varying reference image Q

$$S_ (k) = \hat(D_ |q) = \int_ e^^} (r)}} dr dr$$

(4)

The goal is to gradually disentangle the motion and contrast variations by alternating between the L + S decompositions and the MR-MOTUS framework.

In summary, CMR-MOTUS relies on the alternating solution of two reconstruction problems, namely:

1.

Image reconstruction to determine the time-dependent image contrast changes \(Q\).

2.

Motion estimation to determine the motion fields \(D\).

We briefly describe both alternating steps.

Image reconstruction

In the literature of L + S decomposition in time-resolved MRI [12], the entire target time-series images \(}\) is approximated by a low-rank term \(}\) plus a sparse term \(}\). The \(}\) term captures the larger scale slow changes in the time series, while the \(}\) term captures the finer details in the series, including the contrast enhancement [12].

In this work, we instead leverage motion fields using MR-MOTUS to motion correct the L+S reconstruction process. This is achieved by using the CMR-MOTUS model presented in eq (4) in the data consistency term of the original L+S inverse problem [12].

$$[L^ ,S^ ] = \arg \min _} \sum\limits_^ |D_ ) - s_ } \right\|_^ } + \lambda _ \left\| L \right\|_ + \lambda _ \left\| \right\|_ } Q = L + S$$

(5)

where \(}+}\) now approximates the time-varying reference image \(}\), and the Fourier transform now contains the time-dependent motion fields as described by \(\mathbf\). Note that the variables in the notation \(F(q_ |D_ )\) are switched, now meaning that the model is conditional on the specific value of \(}\). \(\left\| }} \right\|_\) denotes the nuclear norm, \(\left\| }} \right\|_\) the L1-norm, and \(_}\) and \(_}\) are the weights for the regularization terms involving \(}\) and \(}\), respectively. Additionally, we take advantage of sparsity in the temporal Fourier domain using the temporal Fourier operator \(\mathcal\), as it has shown promising results in first-pass myocardial perfusion applications[12].

Technically, we have the flexibility to designate any motion state at time \(t\) as a reference motion static state by subtracting its influence on all the other motion fields. For instance, the reference static phase could correspond to the diastolic one. This adjustment is done before solving Eq. 6 and it steers the L + S reconstruction process towards a diastolic reference motion static state. In this work, we opt to reconstruct a mid-position \(}\) by subtracting the temporal mean of the motion fields from \(}\).

Since we use a non-differentiable regularizer, we use a custom proximal gradient method based on the Fast Ite (rative/Shrinkage Thresholding Algorithm (FISTA) [21, 22] to solve equation (6).

Motion estimation

As described in Huttinga et al. [17], non-rigid motion fields in MRI can be decomposed into a low-rank model with principal components\(\Phi\) and \(\Psi\). For an explicit Rank \(R\),\(\Phi\) is a matrix with size \(2N \times R\) (\(3N \times R\)for 3D), describing the spatial components and \(\Psi\)is a matrix with size \(M\times R\) describing the temporal components. The full motion fields are then approximated by \(D \approx \Phi \Psi^}}\) . Additionally, we parameterize the low-rank motion fields using a B-spline basis in space and time as described by Huttinga et al. [17]. Hence, we can solve the following inverse problem for an explicit low-rank representation of the motion fields, similar to the low-rank MR-MOTUS method, but using a time-varying contrast image qt rather than a fixed reference:

$$[\Phi^ ,\Psi^ ] = \arg \,\min_ }} \sum\limits_^ \left| ) - s_ \left\| ^ } \right.} \right.} \right.} + \left. } \right\|TV\left. \right\|_$$

(6)

The 2, 1 norm refers to the “mixed norm” of the isotropic total variation TV as defined in [19, 23]

$$TV: = \sqrt ^ ^ ]_ } \right\|_ } } \right)^ + \left( ^ ^ ]_ } \right\|_ } } \right)^ }$$

(7)

We use the L-BFGS [24] algorithm to solve equation (67).

Figure 1 summarizes the CMR-MOTUS framework. The iterative process alternates between steps 1 and 2 until the time-varying reference image of the motion state is accurately captured in step 1 and the motion dynamics are effectively captured by the motion fields in step 2. The optimal regularization weights \(_}}\), \(_}}\), \(_}}}\), and the maximum number of alternations and inner iterations are determined empirically (see Methods section).

Fig. 1

Schematic representation of the proposed CMR-MOTUS framework. The framework takes raw time-resolved \(t\) MRI data \(_\). In step 1, the contrast-varying reference images \(_\) are solved while the motion fields \(D\) are fixed, using an L + S reconstruction. In step 2, the motion fields \(D\) are solved while the contrast-varying reference images \(q\) are fixed. The output of step 1, shown on the right, includes the contrast component \(S\) and structure component \(L\). Step 1 and step 2 are solved iteratively in an alternating manner until the motion and contrast are disentangled. The output of step 2, shown on the left, includes the motion fields \(D\)