Multi-tasking to Correct: Motion-Compensated MRI via Joint Reconstruction and Registration
Veronica Corona Angelica I. Aviles-Rivero Noemie Debroux Martin Graves Carole Le Guyader Carola-Bibiane Schoenlieb Guy Williams
























V. Corona et al.   
 Motion-Compensated MRI via Joint Reconstruction and Registration








DAMTP, University of Cambridge,  United Kingdom

DPMMS, University of Cambridge,  United Kingdom

Department of Radiology, University of Cambridge, United Kingdom


LMI, Normandie Université, INSA de Rouen,  France



Department of Clinical Neurosciences, University of Cambridge, United Kingdom












Multi-tasking to Correct: Motion-Compensated MRI via Joint Reconstruction and Registration
    Veronica Corona1Angelica I. Aviles-Rivero2Noémie Debroux1Martin Graves3 Carole Le Guyader4Carola-Bibiane Schönlieb1 Guy Williams5 
    December 30, 2023
======================================================================================================================================

va-ria-tional
quasicon-vexity
si-gnificant
analy-sis
facili-tate
di-sease
reg-is-tra-tion
diffe-rence
dissimila-rity
desi-gning
hyper-elastic
pro-per-ty
li-te-ra-tu-re
boun-da-ry
mi-ni-mi-zing
e-ve-ry
fea-tu-re
pe-na-li-zing
ne-ga-ti-ve
pro-blem
regu-larization
mo-del-ling
drop-ped
lo-ca-ting
ref-er-ence
min-i-mi-za-tion
em-pir-i-cal-ly



This work addresses a central topic in  Magnetic Resonance Imaging (MRI) which is the motion-correction problem in a joint reconstruction and registration framework. From a set of multiple MR acquisitions corrupted by motion, we aim at - jointly - reconstructing a single motion-free corrected image and retrieving the physiological dynamics through the deformation maps. To this purpose,  
we propose a novel variational model. First, we introduce an L^2 fidelity term, which  intertwines  reconstruction and registration along with  the weighted total variation. Second, we introduce an additional regulariser which is based on the hyperelasticity principles to allow large and smooth deformations. We demonstrate through numerical results that this combination creates synergies in our complex variational approach resulting in higher quality reconstructions and a good estimate of the breathing dynamics. We also show that our joint model outperforms in terms of contrast, detail and blurring artefacts, a sequential approach.





§ INTRODUCTION

Magnetic Resonance Imaging (MRI) is a well-established modality that allows to capture details of almost any organ or anatomical structure of the human body. However, the lengthy period of time needed to acquire the necessary measurements to form an MR image continues to be a major contributor to image degradation, which compromises clinical interpretation <cit.>. This prolonged time makes MRI highly sensitive to motion mainly because the timescale of physiological motion -including involuntary motion- is generally shorter than the acquisition time <cit.>. This motion is manifested as undesirable artefacts including geometric distortions and blurring, which causes a significant degradation of the image quality and affects the clinical relevance for diagnosis. Therefore, the question of how to do motion correction in MRI? is of great interest at both theoretical and clinical levels, and this is the problem addressed in this paper. 

Although motion reduction strategies such as faster imaging (e.g. <cit.>) have shown potential results, we focus on post-acquisition correction of MRI sequences, which involves explicit estimation of the motion. In this context, there have been different attempts for motion correction in MRI, whose predominant scheme is based on image registration after reconstructing the images (i.e. to estimate a mapping between a pair of images). A set of algorithmic approaches use rigid registration including <cit.>. 
However, the intrinsic nature of the problem yields those techniques to fail to capture complex transformations - e.g.  cardiac and breathing motion - which, in many clinical cases, are inevitable. To mitigate this limitation, different non-rigid image registration methods have been proposed in this domain - e.g. <cit.>. 

On the other side of the aforementioned methods that separately compute the image reconstruction and the motion estimation, there is another type of algorithmic approaches that computes those sub-tasks jointly. This philosophy was early explored in <cit.>. In that work, authors proposed a factorisation method, which uses the singular-value decomposition for recovering shape and motion under orthography. Particularly in the medical domain, and following a variational perspective <cit.>, different works have been reported including for SPECT imaging <cit.>, PET <cit.> and  MRI <cit.>- to name a few. The works with a closer aim to ours are discussed next. 



Schumacher et al. <cit.> presented an algorithmic approach that combines reconstruction and motion correction for SPECT imaging. The authors proposed a novel variational approach using a regulariser that penalises an offset of motion parameter to favour a mean location of the target object. However, the major limitation is that they only consider rigid motions.
In the same spirit, Fessler <cit.> 
proposed a generic joint reconstruction/registration framework with a model based on a penalised-likelihood functional using a weighted least square fidelity term along with a spatial and a motion regulariser. 
Authors in <cit.> propose a joint model composed of a motion-aware likelihood function and a smoothing  term for a simultaneous image reconstruction and motion estimation for PET data. More recently, Odille et al <cit.> proposed a joint model for MRI image reconstruction and motion estimation. This approach allows for an estimate of both intra and inter-image motion, meaning that, not only the misalignment problem is addressed but also it allows correcting for blurring/ghosting artefacts. 

Most recently and motivated by the deep learning (DL) performance breakthrough, a set of algorithmic approaches has been devoted to investigate the benefits of DL for image registration - e.g. <cit.>. Although, certainly, those approaches deserve attention, their review goes beyond the scope of this paper.

Contributions. In this work, we follow the joint model philosophy, in which we seek to compute - simultaneously and jointly - the MR image reconstruction representing the true underlying anatomy and the registration tasks. Whilst the use of joint models (multi-tasking) has been widely explored for modalities such as PET and SPECT, the application for MRI is a largely untouched area. This is of great interest as 
MRIs highly sensitive to motion. With this motivation in mind we propose a variational model that allows for intra-scan motion correction in MRI. The significance of our approach is that by computing the reconstruction and registration tasks jointly one can exploit their strong correlation thus reducing error propagation and resulting in a significant motion correction. Whilst this is an important part of our solution, our main contributions are:


    
  * A novel theoretically well-motivated mathematical framework for motion correction in MRI. It relies on multi-tasking to address reconstruction and registration jointly, and guarantees the preservation of anatomical structures and topology by ensuring the positivity of the Jacobian determinant.
   

 
  
   
    
  * Evaluation of our theory on a set of numerical and visual results using MRI datasets. We also show that our joint model leads to a better output than a sequential approach using the FAIR toolbox with similar characteristics. 




§ JOINT MODEL FOR MOTION CORRECTION IN MRI
  
In this section, we  describe our approach for motion correction using MRI data. 
Firstly, we theoretically motivate our joint reconstruction and registration model and secondly, we formulate our joint variational model for motion correction. 
Finally, we describe an optimisation scheme to numerically solve our proposed joint model.




 §.§ Mathematical Model and Notation

In an MRI setting, a target image u representing a part of the patient body  is acquired in spatial-frequency space.  Let us denote by x_i∈ L^2(ℝ^2) the i-th  measurement acquired by an MRI scanner for i=1,⋯,T, where T is the number of acquisitions. The measured x_i can be modelled as x_i=𝒜u_i+ε, where 𝒜 : L^2(ℝ^2) → L^2(ℝ^2) is the forward Fourier transform and by 𝒜^*:L^2(ℝ^2)→ L^2(ℝ^2) the backward inverse Fourier transform. 


𝒜 is a linear and continuous operator according to Plancherel's theorem. In this work, we assume to have real-valued images, therefore, our MRI forward operator takes only the real part of complex-valued measurements.

Let Ω be the image domain, a connected bounded open subset of ℝ^2 of class 𝒞^1, and u:Ω→ℝ be the sought single reconstructed image depicting the true underlying anatomy. For theoretical purposes, we assume that it is equal to 0 on the boundary of Ω, so that, we can extend it by 0 and apply the operator 𝒜.

We introduce  the unknown deformation as φ_i : Ω̅→ℝ^2, between the i-th acquisition and the image u. Moreover, let z_i be the associated displacements such that φ_i=+z_i, where  is the identity function. At the practical level, these deformations should be with values in Ω̅. Yet, at the theoretical level, working with such spaces of functions result in losing the structure of vector space. Nonetheless, we can show that our model retrieves deformations with values  in Ω̅ -  based on Ball's results <cit.>. A deformation is a smooth mapping that is topology-preserving and injective, except possibly on the boundary if self-contact is allowed. We consider  ∇φ_i : Ω→ M_2(ℝ) to be the gradient of the deformation, where M_2(ℝ) is the set of real square matrices of order two. After introducing the set of notations for this work, we now turn to to motivate the introduction of our mathematical model.



 §.§ Joint Variational Model


In this work, we seek to extract - jointly -  from a set of MR acquisitions x_i which are corrupted by motion, a clean static free motion-corrected reconstructed image u, along with the physiological dynamics through the deformation maps φ_i. 
With this purpose in mind, we consider the motion correction issue as a registration problem. In a unified variational framework, u and all φ_i's for i=1,⋯,T, are seen as optimal solutions of a specifically designed cost function composed of: (1) a regularisation on the deformations φ_i which prescribes the nature of the deformation, (2) a regularisation on u that allows for removing artefacts whilst keeping fine details and  (3) a fidelity term measuring the discrepancy between the deformed mean - to which the forward operator has been applied -  and the acquisitions, intertwining then the reconstruction and the registration tasks.

In order to allow for large and smooth deformations, we propose viewing the shapes to be matched as isotropic, homogeneous  and hyperelastic materials. This outlook dictates the regularisation on the deformations φ_i, which is based on the stored energy function of an Ogden material. We introduce the following notations : A : B =tr(A^TB) is the matrix inner product whilst A=√(A:A) is the associated Frobenius matrix norm. In two dimensions, the stored energy function of an Ogden material, in its general form, is given by the following expression. W_O(F)=i=1M∑a_iF^γ_i+Γ(F), with a_i>0, γ_i≥ 1 for all i=1,⋯,M and Γ  :  ]0;∞[ →ℝ, is a convex function satisfying δ→ 0^+limΓ(δ)=δ→ +∞limΓ(δ)=+∞. 

In this work, we focus on the particular energy, which reads: 

W_Op(F)={[ a_1F^4 + a_2( F - 1/F)^4                   if F>0;                       +∞               otherwise, ]. 

with a_1>0, and a_2>0. The first term penalises the changes in length, whilst the second term enforces small changes in area. We check that this function falls within the general formulation of the stored energy function of an Ogden material: W_Op(F) = W̃(ξ,δ)={[ a_1ξ^4 + a_2( δ - 1/δ)^4                  if δ>0,;                       +∞               otherwise, ]. 
W̃ : M_2(ℝ) ×ℝ→ℝ is continuous since δ→ 0^+limW̃(ξ,δ)=δ→ +∞limW̃(ξ,δ)=+∞, and is convex with respect to δ (g: | [     ℝ^+_* →ℝ; x ↦(x-1/x)^4 ]. and g” : | [                         ℝ^+_* →ℝ; x↦ 4(x-1/x)^2(3x^4+4x^2+5/x^4)≥0 ].). The design of the function Γ is driven by Ball's results <cit.> guaranteeing the deformation to be a bi-Hölder homeomorphism, and therefore, preserving the topology. It also controls that the Jacobian determinant remains close to one to avoid expansions or contractions that are too large.

The aforementioned regulariser is then applied along a discrepancy measure, which allows to interconnect the reconstruction and the registration tasks, and a regularisation of our free motion-corrected slice based on the weighted total variation.
With this purpose, let g : ℝ^+ →ℝ^+ be an edge detector function satisfying g(0)=1, g strictly decreasing and r→ +∞limg(r)=0. For the sake of simplicity, we  set g_i = g( ∇ G_σ *𝒜^*x_i), with G_σ being a Gaussian kernel of variance σ and for theoretical purposes, we assume that there exists c>0 such that 0<c≤ g_i≤ 1, and g_i is Lipschitz continuous for each i=1,⋯,T. In practice, we choose g to be the Canny edge detector. 
We follow Baldi's arguments (<cit.>) to introduce the weighted BV-space and the associated weighted total variation related to the weight g_i, for each i=1,⋯,T.

Let w be a weight function satisfying some properties (defined in <cit.> and fulfilled here by g_i). We denote by BV_w(Ω) the set of functions u∈ L^1(Ω,w), which are integrable with respect to the measure w(x)dx, such that:

    sup{∫_Ω u (φ) dx :  |φ|≤ w  everywhere , φ∈ Lip_0(Ω,ℝ^2) }<+∞,

with Lip_0(Ω,ℝ^2) the space of Lipschitz continuous functions with compact support. We denote by TV_w the previous quantity.



Given E a bounded open set in ℝ^2, with boundary of class 𝒞^2, then we have TV_g_i(ξ_E)=|∂ E|(Ω,g_i)=∫_Ω∩∂ E g_i dH^1, where ξ_E is the characteristic function of the set E. This quantity can be viewed as a new definition of the curve length- with a metric depending on the observations x_i. Minimising it is equivalent to locating the curve on the edges of 𝒜^*x_i where g_i is close to 0. We thus introduce the following fidelity term and regulariser for our single reconstruction:

    F(u,(φ_i)_i=1,⋯,T)= 1/Ti=1T∑δ TV_g_i(u ∘φ_i^-1)+1/2𝒜((u∘φ_i^-1)_e)-x_i_L^2(ℝ^2)^2,

where φ_i^-1 is the inverse deformation 
and (u∘φ_i^-1)_e is the extension by 0 of (u∘φ_i^-1) outside Ω. Indeed, we assume that u=0 and φ_i=φ_i^-1= so that u∘φ_i^-1=0 on the boundary ∂Ω. Moreover, if u∘φ_i^-1∈ BV_g_i,0(Ω), with BV_g_i,0(Ω)=BV_0(Ω) ∩ BV_g_i(Ω) (see <cit.> for more details on the trace operator on the space of functions of bounded variations), then u∘φ_i^-1 is null on the boundary ∂Ω and can be extended by 0 outside the domain Ω.
Due to the embedding theorem (<cit.>), we have that (u∘φ_i^-1)_e∈ L^2(ℝ^2) and then Plancherel's theorem gives us 𝒜((u∘φ_i^-1)_e)∈ L^2(ℝ^2) leading to the well-definedness of the fidelity term.
The first term of F aims at aligning the edges of the deformed u (u∘φ_i^-1) with the ones of the different acquisitions, whilst regularising the reconstructed image u. The second quantity targets to make the Fourier transform of the deformed u (u∘φ_i^-1) as close as possible to the acquisitions x_i (in the L^2 sense) - inspired by the classical fidelity term in MRI reconstruction models.

By combining both terms, we finally get the following minimisation problem

    inf I(u,(φ_i)_i=1,⋯,T) = F(u,(φ_i)_i=1,⋯,T)+1/Ti=1T∑∫_Ω W_Op(∇φ_i) dx. P












We now describe the numerical resolution of our joint model (<ref>).


 §.§ Numerical Method: Optimisation Scheme


In order to tackle the nonlinearity of the problem, we introduce three auxiliary variables v_i,  w_i,f_i to mimic ∇φ_i, u∘φ_i^-1 and w_i, respectively, through quadratic penalty terms. This corresponds to reformulating the equality constraints v_i=∇φ_i, w_i=u∘φ_i^-1 and w_i=f_i, into an unconstrained minimisation problem.
We then propose solving the following decoupled (discretised) problem

    min_u,φ_i,v_i,w_i,f_i    1/T∑_i=1^T x ∈Ω∑ W_Op(v_i(x)) + γ_1/2v_i - ∇φ_i_2^2 
       + γ_3/2𝒜w_i -x_i _2^2 + γ_2/2(w_i - u ∘φ_i^-1) √(det∇φ_i^-1)_2^2 
       + 1/2θf_i -w_i_2^2 +TV_g_i(f_i).

We can now minimise (<ref>) with an alternating splitting approach, where we solve the subproblem according to each individual variable whilst keeping the remaining variables fixed. 
We iterate the following updates for k=0,1,…. 

Sub-problem 1: Optimisation over v_i.  In practice, v_i simulates the gradient of the displacements z_i=(z_i,1,z_i,2) associated to the deformations φ_i. For every v_i, we have v_i=[ v_11  v_12; v_21  v_22 ].  For the sake of readability, we drop here the dependency on i. We solve the Euler-Lagrange equation with a semi-implicit finite difference scheme (for N iterations  and all i) and update v_i as: 

    v_11^k+1   =1/1+dt γ( v_11^k + dt(-4a_1 v_i^k^2_F(v_11^k+1) - 4 a_2 (1+v_22^k)  c_0 c_1
    + γ∂ z_1^k/∂ x), 
    
    v_12^k+1   =1/1+dt γ( v_12^k + dt(-4a_1 v_i^k^2_F v_12^k + 4  a_2 v_21^k c_0 c_1
    + γ∂ z_1^k/∂ y), 
    
    v_21^k+1   =1/1+dt γ( v_21^k + dt(-4a_1 v_i^k^2_F v_21^k + 4  a_2 v_12^k c_0 c_1
    + γ∂ z_2^k/∂ x), 
    
    v_22^k+1   =1/1+dt γ( v_22^k + dt(-4a_1 v_i^k^2_F(v_22^k+1) - 4 a_2 (1+v_11^k)  c_0 c_1+ γ∂ z_2^k/∂ y),

with c_0=(detv_i^k- 1/detv_i^k)^3 and c_1=1+ 1/(detv_i^k)^2.

Sub-problem 2: Optimisation over φ_i. We solve the Euler-Lagrange equation in φ_i, after making the change of variable y=φ_i^-1(x) in the L^2 penalty term, for all i, using an L^2 gradient 
flow scheme with a semi-implicit Euler time stepping (for N iterations in the same inner loop as v_i):

    0=γ_1 Δφ_i^k+1 + γ_1 [ divv_i,1^k+1; divv_i,2^k+1 ] + γ_2 (w_i^k ∘φ_i^k - u^k) ∇ w_i^k(φ_i^k).


Sub-problem 3: Optimisation over w_i. The update in w_i, for all i, reads:

    w_i^k+1 = 𝒜^*{𝒜(γ_2 det∇ (φ_i^-1)^k+1 u ∘ (φ_i^-1)^k+1 + f_i^k/θ)+ γ_3 x_i/γ_2det∇(φ_i^-1)^k+1+γ_3+1/θ}.


Sub-problem 4: Optimisation over f_i. For each  k and i, this is solved via Chambolle projection algorithm <cit.>, as in <cit.>, for an inner loop over n=0,1,…

    f_i^n+1   =w_i^k+1-θdivp_i^n, 
    
        p_i^n+1   =p_i^n+ δ_t ∇ (div p_i^n - w_i^k+1/θ)/1+δ_t/g_i | ∇(div p_i^n - w_i^k+1/θ)| .

After enough iterations, we set f_i^k+1=f_i^n+1.

Sub-problem 5: Optimisation over u. Finally, using the same change of variables as in sub-problem 2, the problem in u is simply

    u^k+1=1/T∑_i=1^T w_i^k+1∘φ_i^k+1.


In the following section, we present the results obtained with our algorithmic approach, where the registration problem in v_i and φ_i, i=1,…,T, has been solved in a multi-scale framework, from coarser to finer grids. It is worth noticing that - even though in W_Op, we have a term controlling that the Jacobian determinant remains positive, we introduce an additional step in our implementation, to monitor that the determinant does not become negative. 
To do so, we use the following regridding technique:
if at iteration N, det(∇φ^N) becomes too small or negative, we save  φ^k,N-1 and re-initialise φ^k,N=Id and v^k,N=0, and w^k = w^k ∘φ^k,N-1, and we then continue the loop on N. Finally, φ^k,final is the composition of all the saved deformations.



§ NUMERICAL EXPERIMENTS
 
This sections describes the set of experimental results that are performed to validate our proposed joint model. 

Data Description. We evaluate our framework on two datasets. The Dataset I[<http://www.vision.ee.ethz.ch/ organmot/chapter_download.shtml>] is a 4DMRI data acquired during free-breathing of the right liver lobe <cit.>. It was acquired on a 1.5T Philips Achieva system, TR = 3.1 ms, coils =4, slices =25,  matrix size = 195×166. Dataset II[<https://zenodo.org/record/55345#.XBOkvi2cbUZ>] is a 2D T1-weighted dataset <cit.> acquired during free breathing of the entire thorax. It was acquired with a  3T Philips Achieva system  with matrix size = 215 × 173,  slice thickness= 8mm, TR=3.1ms and TE=1.9ms. The experiments reported in this section were run under the same condition in a CPU-based Matlab implementation. We used an Intel core i7 with 4GHz and a 16GB RAM.  

l0.5


    0.5!  
    
    2-11
          [HTML]EFEFEF a_1     [HTML]EFEFEFa_2     [HTML]EFEFEFγ_1     [HTML]EFEFEFγ_2     [HTML]EFEFEFγ_3     [HTML]EFEFEFθ     [HTML]EFEFEFσ     [HTML]EFEFEFk     [HTML]EFEFEFN     [HTML]EFEFEFn

    [HTML]EFEFEFDataset I         1    50    5    10^5    30    5    2    2    500   500  

    [HTML]EFEFEFDataset II          1    50    5     10^5    15     5    1.5    2     500    500 

    
    
    Parameter values.
    


 
Parameter Selection. 
We present a discussion on the role of parameters in our model. The ranges of these parameters are rather stable for both experiments as shown in Table <ref>.
The parameters in the regularisation term for the registration W_Op are a_1 and a_2. The former controls the smoothness of the deformation whilst the latter is more sensitive and can be viewed as a measure of rigidity allowing for a trade-off between topology preservation and the capacity of handling large deformations. 
If a_2 is chosen too small, the deformation may lose its injectivity property resulting in the loss of topology preservation, while if chosen too big, the registration accuracy may be impaired. 
The parameters γ_1,γ_2 and θ, balance the L^2 penalty terms. 
Finally, γ_3 weights the data-fitting term for the reconstruction. For the FAIR approach - which is described next - we set the regularisation parameter for Dataset I=0.1 and DatasetII=1. 



Comparison with FAIR. We begin by evaluating our joint model against a sequential approach using the FAIR toolbox <cit.>, which is a well-established framework for image registration. For a fair comparison, we set the characteristics of the compared approach closer to ours, as follows. Firstly, we choose the sum of squared differences as a distance measure on  A^*x_i, and secondly, we set the hyperelastic regulariser - as it allows for similar modeling of the deformations.

The computations of FAIR and our approach were performed under the same conditions.
We ran the same multi-scale approach. Moreover, we set the same reference image that we selected for the initialisation of our alternating scheme. 
This is important because of the non-convexity of our joint problem. After registering each frame to the reference, we average the results to obtain a mean image. In Fig. <ref>, we show the outputs - for both datasets - uncorrected being the coarse Euclidean mean, FAIR referring to the results of the process described above; and ours being the reconstructed image u obtained with our model. 
By visual inspection, we observe that the coarse Euclidean mean displays highly blurry effects along with loss of details whilst our approach improves in this regard. 
We also observed that - although FAIR improves over the coarse mean - outputs generated with FAIR have a remarkable loss of contrast  and detail. In addition, we point out that even though our method needs to solve 5 subproblems, the final computational time is less than the sequential FAIR approach.



Overall, visual inspection of the results shows the following drawbacks in which our joint model improves: (1) blurry effects: for both datasets, our model gives clearer images compared to the other approaches in which small details and shapes are lost, and (2) contrast preservation: our method is able to better preserve the contrast than other approaches, especially the one based on FAIR in which loss of contrast is clearly visible.


Further Analysis of our Approach. In Fig. <ref>, we can observe that our method produces a reasonable estimation of the breathing dynamics, where the motion fields are visualised by a deformation grid. More precisely, it displays - for the two datasets - the estimated  motion φ_i and its inverse φ_i^-1 for the frame i=5 along with their corresponding Jacobian determinant maps, whose values are interpreted as follows: small deformations when values are closer to 1,  big expansions when values are greater than 1, and big contractions when values are smaller than 1.

Moreover, one can observe that the determinants remain positive in all cases, that is to say, our estimated deformations are physically meaningful and preserve the topology as required in a registration framework.

In Fig. <ref>, we show further results for the registration.
We display a difference map (top) between our (adaptive) reference u, and the registered data w_i ∘φ_i for sample frames for Dataset I and II. 
Moreover, we show the mutual information (MI) between the same images as a quantitative measure of performance (bottom), 
where larger MI values indicates better alignment. Since the L^2 fidelity term drives the registration process in our optimisation framework, it sounds relevant to use the difference map to illustrate the accuracy of our model. However, one should also be aware that a small change in contrast - which often happens between two frames and thus between u and w_i ∘φ_i - is clearly visible on the difference map, even if the overall anatomies are well aligned. Keeping in mind that, for the sake of tractability, the parameters are set globally and not adjusted for each frame, we can see that the intensities and the structures are well aligned for all the examples. This is corroborated by large MI values.





§ CONCLUSION

In this work we presented a new mathematical model to tackle motion correction in MRI in a joint reconstruction/registration framework. The underlying idea is to exploit redundancy in the temporal resolution in the data to compensate for motion artefacts due to breathing. The proposed multi-tasking approach solves simultaneously the problems of reconstruction and registration to achieve higher image quality and accurate estimation of the breathing dynamics. We mathematically motivated our joint model and presented our optimisation scheme. Furthermore, our experimental results show the potential of our approach for clinical applications. 



§ ACKNOWLEDGEMENTS
 VC acknowledges the financial support of the Cambridge Cancer Centre. Support from the CMIH and CCIMI University of Cambridge is greatly acknowledged.


splncs03



