In this section, we present the basic mathematical formulation underpinning the proposed motion correction method (further details can be found in [29]). The contrast volume, in the remainder of this section, will be denoted by \(}} \in }}^\), where \(n_x\) is the number of voxels contained in a rectangular field of view. The 3D image undergoes a time-dependent rigid motion

where \(t\) is a time-related label. In practice, \(t\) corresponds to the index of the k-space readout line in the phase-encoding plane. The corresponding rigid transformation is given by \(T_}}_t }\), and is parameterized by a time-dependent motion parameter \(}}_t \in }}^6\), which includes translations and rotations in 3D:

$$}} = \left( }} ,}} \right),\quad }} = \left( \right),\quad } = \left( \right).$$

(2)

The rigid motion consists of a 3D rotation (defined by the 2D rotation angles \(\varphi_x , \varphi_y , \varphi_z\), performed in this order around the corresponding axes) followed by a translation (governed by the translation parameters \(\tau_x , \tau_y , \tau_z\).

Without loss of generality, we are assuming a Cartesian acquisition. At each given time \(t\), the MR acquisition process corresponds to the evaluation of the Fourier transform \(}}\) of \(u_t\) in a particular subset \(K_t\) of the k-space. In practice, the acquisition is structured in such a way that all the subsets \(K_t\) consist of parallel lines in the k-space (the common direction being the readout direction). We refer to the Fourier transform of a rigidly moving object \(u_}} : = T_}} u\) as the perturbed Fourier transform \(}}_}} u: = }}u_}}\), and can be directly characterized as

$$}}_}} u \left( }} \right): = }\left( }k \cdot \tau } \right)}}u \left( }^ u} \right)$$

(3)

where the rotational operator with respect to the 3D angle φ is indicated by \(R_}^\). This definition is motivated by classical Fourier identities that describe the action of rigid motion under the Fourier transform. Due to rotational effects, one must resort to the non-uniform discrete Fourier transform (NUFFT) to evaluate Eq. (3) [33, 34].

Note that we implicitly assumed that no motion occurs while sampling the elements of \(K_t\), since the state of the object at the time \(t\) is associated to a single motion parameter \(}}_t\). The assumption is motivated by the fact that \(K_t\) will correspond, in practice, to a single Cartesian readout line, which lasts few milliseconds. Hence, the data acquisition at time \(t\) is symbolized by the application of the selection operator \(S_t\) to the Fourier-transformed volume:

$$d_t = S_t }}_}}_t } \mathbf = \left( }}_}}_t } \mathbf\left( _1 } \right),...,}}_}}_t } \mathbf\left( _ } \right)} \right),\quad \mathbf_1 ,...,\mathbf_ \in K_t .$$

(4)

Here, \(n_r\) is the number of k-space samples in a single readout. The resulting inverse problem can be cast as an optimization problem over the reconstruction unknowns u and the motion parameters \(}}_t\), that is:

$$\limits_, }}_ }} f\left( , }}_ } \right) + \mu\, g_\theta \left( }}_ } \right)\quad }}.\quad g_u \left( \mathbf \right) \le \varepsilon$$

(5)

where \(}}_ = ( }}_1 , \ldots , }}_ )\), and \(n_t\) is the number of time steps. The parameters \(\varepsilon , \mu\) (both positive numbers) set the strength of the corresponding regularization terms. The first term of the objective functional in Eq. (5) corresponds to the data misfit:

$$f\left( , }}_ } \right) = \mathop \sum \limits_^ \frac\|}}_}}_t } } - \mathbf_t \|^2 .$$

(6)

The least-squares norm is indicated here by ∥·∥. The regularization terms \(g_u\) and \(g_\theta\) are crucial in ensuring the well-posedness of the problem. Indeed, the objective in Eq. (6) will be sensitive to the relatively low signal-to-noise ratios of the high-frequency components of the data. Moreover, the objective is highly non-convex as a function of \(}}_\). The motion-parameter regularization is designed to ensure some form of regularity in time (e.g., smoothness), this can be achieved for example by setting

$$g_\theta \left( }}_ } \right) = \mathop \sum \limits_^ \frac\|}}_ - }}_t \|^2 .$$

(7)

Alternatively, higher-order derivatives may be used. Another strategy, adopted in this paper, is to impose smoothness by setting hard constraints for the motion parameters, rather than via an additive penalty term as in Eq. (7) [29].

The crux of the proposed method is related to the choice of the regularization term \(g_u\) in Eq. (5). We adopt the structure-guided total variation scheme proposed in [30] in the context of multi-contrast imaging, that is:

$$g_u \left( \mathbf \right) = \sum \limits_}\| _} |_} \nabla \mathbf|_}\| ,\quad\quad _} |_} = I_3 - _} |_} \upxi_} |_}^} ,$$

(8)

where \(I_3\) is the 3 × 3 identity matrix, \(\nabla \cdot |_}\) is the discretized gradient operator evaluated at the voxel with center x, and \(_} |_}\) is the projection operator on the linear space that is orthogonal to the vector \(_} |_} \in }}^3\). The symbol H indicates the adjoint operation. The vector \(_} |_}\) corresponds to the normalized gradient of a given motion-free contrast v, e.g., \(\upxi_} |_} \approx \nabla \mathbf|_} / \|\nabla \mathbf|_}\|\). The actual definition is

$$_} |_} = \frac|_} }|_}\|^2 + \eta^2 } }} ,$$

(9)

for some constant η > 0. The regularization term in Eq. (8) enforces the gradient structure of v onto u, when v and u are anatomically compatible. It is important to observe that v is not required to be registered with the target contrast \(\mathbf\), since the estimation of the motion parameters in Eq. (5) will automatically compensate for the initial misalignment (see also [31]). In this work, we actually adopt a constrained formulation, meaning that structural similarity is imposed by forcing the solution to belong to the constraint set \(C_u = \ : g_u \left( \mathbf \right) \le \varepsilon \}\), where \(\varepsilon > 0\) is a prescribed regularization level (see [29, 35], for more details).

To solve Eq. (5), we adopt an alternating update scheme based on the proximal alternating minimization algorithm (PALM) described in [36]. The algorithm of the optimization strategy is exemplified in Algorithm 1. Each update requires the linearization of the smooth objective \(f\) and the application of the proximal operators associated to \(g_u\) and \(g_\theta\). We will make use of multi-scale methods to ease the ill-posedness of the problem. Two types of scale are considered, here:

spatial/temporal grid: this scale is associated to the spatial and temporal grid sizes of the reconstructed image \(}}\) and motion parameters \(}}_\), respectively, by considering a sequence of optimization problems defined on progressively finer grids. The pragmatic approach considered in this work actually consists in fixing a relatively coarse temporal grid for the motion unknowns, along with a corresponding time-interpolation operator (this procedure effectively acts as an additional implicit regularizer for \(}}_\)). Therefore, only the spatial grid of \(}}\) is scaled up. Note that the spatial scale considered at a certain multi-scale stage poses a limit on how well the motion parameters can be resolved at that stage, due to the Nyquist criterion, since they are associated to specific coordinates in k-space;

regularization strength: this scale is related to the regularization level \(\varepsilon\), as defined in Eq. (5). Hence, strongly regularized problems are solved first, and the regularization is gradually relaxed as in a continuation strategy. As mentioned above, this regularization term is explicitly implemented by forcing the solution to lay in the set \(C_u = \ : g_u \left( \mathbf \right) \le \varepsilon \}\), where \(\varepsilon > 0\) and \(\varepsilon = \alpha\, g_u \left( _}} } \right)\) and \(\alpha = 0.01, 0.1, 0.5, 0.8\). This choice was preliminarily fine-tuned on earlier results. See for more details [29].

Overall, this results in two nested sequences of optimization problems (see Algorithm 1).

Algorithm 1 Joint motion correction and reconstruction with alternating proximal operator evaluation

The proposed motion correction algorithm estimates the rigid motion \(}} = (}} ,})\) that the object of interest undergoes at some point during the scan, to undo its effect on the reconstructed 3D image. We parameterize the rigid motion in terms of three translation parameters \(}} = (\tau_x , \tau_y ,\tau_z )\), each one corresponding to motion in the spatial coordinate direction \(x, y, } z,\) and three rotation angles \(} = \left( \right)\), which describe 2D rotations around the axes \(x, y, z\), respectively. We conventionally assume that the motion is performed by an initial translation, followed by three plane rotations. The order of these rotations is implicitly defined by \(}\).

For a consistent display of our results, we are assuming that: the \(x\) direction corresponds to the left–right direction, \(y\) to the posterior-anterior direction, and \(z\) to the inferior-superior direction, the \(xy\) plane corresponds to the axial plane, \(xz\) to the coronal plane, and \(yz\) to the sagittal plane. Left/right, anterior/posterior, and inferior/superior are meant from the patient perspective. The orientation of the rotation planes is determined by the lexicographic order of their labels, according to the right-hand rule.