Patent Publication Number: US-11645791-B2

Title: Systems and methods for joint reconstruction and segmentation of organs from magnetic resonance imaging data

Description:
RELATED APPLICATIONS 
     The present application claims the benefit of U.S. Provisional Application Ser. No. 62/916,455 filed on Oct. 17, 2019, the entire disclosure of which is expressly incorporated herein by reference. 
    
    
     BACKGROUND 
     Technical Field 
     The present disclosure relates generally to the field of medical imaging and computer vision technology. More specifically, the present invention relates to systems and methods for joint reconstruction and segmentation of organs from magnetic resonance imaging (“MRI”) data. 
     Related Art 
     Magnetic resonance imaging (“MRI”) is widely used due to its high resolution, high contrast, and no radiation, but a fully-sampled MRI scan requires a lot of time for a patient to undergo the scan. Such scans could greatly benefit from advances being made in machine learning techniques (artificial intelligence or “AI”) for MRI reconstruction of under-sampled MRI scan data that specifically benefit cardiac MRIs. In particular, there is a need to improve the performance of reconstruction (making current MRI scans faster), remove the need for manual intervention, and support reconstruction and organ segmentation where there is patient and/or organ movement, thereby improving the ability to reconstruct MRI images of the heart and other organs, as well as performing cardiac segmentation. 
     Better MRI imaging of the heart that requires less time and the ability to use AI techniques to address motion and organ segmentation is of particular importance. One of the biggest hurdles to more extensive use of cardiac MRI is how long the exam takes (over an hour) and how complicated the procedures are (patients need to hold their breath during certain sequences). The breath hold is necessary to cease any motion that might obscure the image, but not all patients are able to do this for the required length of time or repeatedly. There is also a need to reduce manual intervention/correction required by current approaches to myocardium segmentation in cardiac magnetic resonance (MR) images which play a vital role in clinical diagnosis of the cardiovascular diseases. Because of the low contrast and large variation in intensity and shapes, myocardium segmentation has been a challenging task. Finally, there is a need to reduce both scan and reconstruction time. Accordingly, the systems and methods of the present address the foregoing, and other, needs. 
     SUMMARY 
     The present disclosure relates to systems and methods for joint reconstruction and segmentation of organs from magnetic resonance imaging (MRI) data. Sparse MRI data (k-space MRI data) is received at a computer system. The computer system jointly processes the MRI data using a plurality of reconstruction and segmentation processes. For example, the MRI data is processed using a joint reconstruction and segmentation process to identify an organ from the MRI data. Additionally, the MRI data is processed using a channel-wise attention network to perform static reconstruction of the organ from the MRI data. Further, the MRI data can be processed using a motion-guided network to perform dynamic reconstruction of the organ from the MRI data. Advantageously, the joint processing performed by the systems and methods of the present disclosure allow for rapid static and dynamic reconstruction and segmentation of organs from sparse MRI data, with particular advantage in clinical settings. 
    
    
     
       BRIEF DESCRIPTION OF THE DRAWINGS 
       The foregoing features of the invention will be apparent from the following Detailed Description, taken in connection with the accompanying drawings, in which: 
         FIG.  1    is a flowchart illustrating overall processing steps carried out by the systems and methods of the present disclosure; 
         FIG.  2    is a diagram illustrating processing performed by processing step  16   a  of  FIG.  1   , in greater detail; 
         FIGS.  3 A- 3 B  are graphs illustrating the RSME of different structures ( FIG.  3 A ) as well as the training time of different structures ( FIG.  3 B ); 
         FIGS.  4 - 5    illustrate operation of the processing step  16   a  of the system; 
         FIG.  6    is a diagram illustrating processing performed by step  16   b  of  FIG.  1   , in greater detail; 
         FIG.  7    illustrates operation of the processing step of  16   b  of the system; 
         FIGS.  8 - 9    are diagrams illustrating processing performed by step  16   c  of  FIG.  1   , in greater detail; 
         FIGS.  10 A- 10 B  are graphs illustrating performance of the process  16   c  of the system; 
         FIGS.  11 - 12    are images illustrating performance of the process  16   c  of the system; and 
         FIG.  13    is a diagram illustrating the motion-guided network of the present disclosure implemented as a single, differentiable network. 
     
    
    
     DETAILED DESCRIPTION 
     The present disclosure relates to systems and methods for joint reconstruction and segmentation of organs from MRI data, as discussed in detail below in connection with  FIGS.  1 - 13   . 
       FIG.  1    is a flowchart illustrating overall processing steps carried out by the systems and methods of the present disclosure, indicated generally at  10 . In step  12 , the system obtains k-space MRI data  14  from a suitable source, such as directly from an MRI scanner, from a computer system, and/or from a database of MRI data. The k-space MRI data  14  typically is undersampled data generated by the MRI machine when a scan is performed on a patient. In the past, the k-space MRI data  14  often did not have sufficient information to allow for adequate reconstruction and segmentation of organs in clinical applications, including both static and dynamic reconstruction and segmentation. Importantly, the three joint processes  16   a - 16   c  of the systems and methods of the present disclosure, discussed in more detail below, are performed jointly (or, in series) on the k-space MRI data  14  to allow for joint reconstruction and segmentation of organs, including both static and dynamic reconstruction and segmentation. This represents a significant improvement in the speed of MRI scans as well as the ability to perform dynamic and static reconstruction and segmentation of organs using sparse k-space MRI data. The processes  16   a - 16   c  include process  16   a  for performing joint reconstruction and segmentation of the organ (such as the heart or other organ), process  16   b  for performing static reconstruction of the organ using a channel-wise attention network, and process  16   c  for performing dynamic reconstruction of the organ using a motion-guided network. Each of these processes will be described in greater detail below. 
     It is noted that one or more of the processes disclosed herein (e.g., the processing steps shown in  FIG.  1    and more fully described herein with reference to the other figures) could be performed via a processor in communication with, or embedded within, an MRI scanner. Indeed, such processor could be an existing processor of an MRI scanner that is programmed in accordance with the processes disclosed herein. Further, the processes disclosed herein could be performed by a standalone processor that processes MRI image data, such as a personal computer, server, cloud computing device, smart phone, or any other suitable computing device. Further, the processes discussed herein could access MRI image data stored in a memory (e.g., a memory of an MRI scanner or a separate memory), and could be embodied as computer-readable instructions executed by a processor in communication with the memory which process the MRI image data as disclosed herein. 
     Additionally, it is noted that the processes  16   a - 16   c  need not occur jointly, and that each process can be operated independently and with varying input data without departing from the spirit or scope of the present invention. For example, if the input data is static data, process  16   b  (independent of processes  16   a  and  16   c ) can be utilized to process the static data in order to perform reconstruction of an organ using the channel-wise attention network disclosed herein. As another example, if the input data is dynamic data, process  16   c  (independent of processes  16   a  and  16   b ) can be utilized to perform dynamic reconstruction of the organ using the motion-guided network disclosed herein. Still further, any combination of the processes  16   a - 16   c  could be performed in any desire order. For example, the process  16   a  can be used to perform segmentation of the organ using the joint reconstruction and segmentation network disclosed herein, in combination with one, or both, of the processes  16   b - 16   c.    
       FIG.  2    is a diagram illustrating processing performed by process  16   a  of  FIG.  1   , in greater detail. As shown in  FIG.  2   , process  16   a  is a joint image reconstruction and myocardium segmentation deep neural network based approach that has the following structure: an image reconstruction network and a segmentation network. 
     The general objective of a CS-MRI reconstruction problem is to find the best reconstructed image x such that it minimizes the difference between undersampled data in k-space y, and undersampled data converted from x through the Fourier transform. The problem is formalized as follows: 
                       x   *     =       arg   ⁢           ⁢       min   x     ⁢           ⁢       1   2     ⁢           ⁢              E   ⁡     (     F   ⁡     (   x   )       )       -   y          2   2           +     λ   ⁢           ⁢     g   ⁡     (   x   )             ,           (   1   )               
where F is Fourier transform, E is an undersampled matrix and λg(⋅) is a weighted regularization term. In MRI image reconstruction, the regularization function g usually takes    1 -norm or total variation norm.
 
     Let f(x)=½E(F(x))−y∥ 2   2 . Equation (1) immediately above is solvable if f is a smooth convex function with a Lipschitz continuous gradient L(f) and g is a continuous convex function and possibly nonsmooth. For any L&gt;0, consider the following quadratic approximation of Equation (1) at a given point x. 
                         Q   L     ⁡     (     x   ,     x   ′       )       =       f   ⁡     (     x   ′     )       +       (     x   -     x   ′       )     ·     ∇     f   ⁡     (     x   ′     )           +       L   2     ⁢            x   -     x   ′            2       +     λ   ⁢           ⁢     g   ⁡     (   x   )             ,           (   2   )               
which admits a unique minimizer p L . By ignoring constant terms related to x′, we can rewrite p L  as:
 
     
       
         
           
             
               
                 
                   
                     
                       p 
                       L 
                     
                     ( 
                     
                       x 
                       ′ 
                     
                     ) 
                   
                   = 
                   
                     
                       arg 
                       ⁢ 
                       min 
                       ⁢ 
                       
                         { 
                         
                           
                             Q 
                             L 
                           
                           ( 
                           
                             x 
                             , 
                             
                               x 
                               ′ 
                             
                           
                           ) 
                         
                         } 
                       
                     
                     = 
                     
                       arg 
                       ⁢ 
                       
                         
                           min 
                           x 
                         
                         
                           { 
                           
                             
                               
                                 λ 
                                 ⁢ 
                                 
                                   g 
                                   ⁡ 
                                   ( 
                                   x 
                                   ) 
                                 
                               
                               + 
                               
                                 
                                   L 
                                   2 
                                 
                                 ⁢ 
                                 
                                    
                                   
                                     x 
                                     - 
                                     
                                       
                                         ( 
                                         
                                           
                                             x 
                                             ′ 
                                           
                                           - 
                                           
                                             
                                               1 
                                               L 
                                             
                                             ⁢ 
                                             
                                               ∇ 
                                               
                                                 f 
                                                 ⁡ 
                                                 ( 
                                                 
                                                   x 
                                                   ′ 
                                                 
                                                 ) 
                                               
                                             
                                           
                                         
                                          
                                       
                                       2 
                                     
                                   
                                   } 
                                 
                               
                             
                             , 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     According to Equation (3) immediately above and reconstruction problem Equation (1) immediately above, we can compute x k  iteratively using the following minimization approach: 
                       x   k     =     arg   ⁢       min   x       {       λ   ⁢     g   ⁡   (   x   )       +         L   k     2     ⁢            x   -     (       v   k     -       1     L   k       ⁢     ∇     f   ⁡   (     v   k     )           )            2         }           ,           (   4   )               
where v k  is a linear combination of x k−1  and x k−2  from the previous two iterations and L k  is a constant chosen by specific rule in the backtracking FISTA. This can generate a sequence of {x 1 , x 2 , . . . , x k }. In theory, if k goes to infinity, then x k  can approximate the optimal solution of Equation (1) immediately above. In order to speed up the converge, x k  is replaced by v k , a linear combination of x k−1  and x k−2 .
 
     Note that in our MRI reconstruction problem, the gradient of f(v k ) is defined as: ∇f(v k )=F T (E(F(v k ))−y), where F T  represents inverse Fourier transform. For specific forms of g(⋅) such as  1 norm, we can derive a closed form of x k  and then obtain the solution of x k  through gradient descent. More generally, we use gradient descent to generate a sequence of {z k,n } that gradually converges to x k . Then we can generalize the problem to any form of g(⋅). Basically, z k,n  is computed as follows:
 
 z   k,n   =z   k,n−1   −μ{λ∇g ( z   k,n−1 )+ L   k ( z   k,n−1   −v   k )+ F   T ( E ( F ( v   k ))− y )},  (5)
 
where μ denotes the gradient step and ∇g(⋅) denotes the gradient of g(⋅).  FIG.  2    shows the relationship between z and x. Suppose the input is x i . Using Equation (5) immediately above iteratively n times, we can estimate x i+1 . Hence, there are n steps between x i  and x i+1 . z i,j  denotes the j th  step between x i  and x i+1 .
 
     However, there are several drawbacks in this FISTA-based computational framework. First, it is very hard to derive a closed form of x k  when g(⋅) is one of such joint norms. Second, it requires a huge amount of computation to find a good solution using the iterative process. Furthermore, it is time-consuming to find a good step size μ, since it varies given different undersampling rates. Additionally, L k  is chosen according to rigorous conditions in backtracking FISTA, so it is sometimes hard to satisfy all conditions. 
     In order to overcome the aforementioned difficulties, we first propose a deep neural network model called FR-Net for MRI reconstruction. Specifically, we create a corresponding neural network for x k  that can be solved via gradient descent. One advantage is that the neural network can be adapted to different forms of g(⋅) and hence closed form derivation is no longer required. Suppose at step k, the optimal value of Equation (4) immediately above is x k . At each iteration n in the gradient descent, a value z k−1,n  is generated to approximate x k . As n becomes larger, z k−1,n  converges to x k .  FIG.  2    shows how to iteratively compute z and x. Process  16   a  includes the use of a joint reconstruction and segmentation deep neural network  20  which processes undersampled k-space data  22 , an FR-Net  24  which generates a reconstructed image  26 , and a U-net (discussed below)  28  which generates an image  30  indicating myocardium contours. Process  16   a  also includes an iteration block  32  and processing chain  34 , both discussed below. 
     In order to solve Equation (5) immediately above, we need to come out the close form of ∇g(z i,j ). However, the regularization term g is always complicated. For example, it can be a joint norm whose gradient is hard to derive. Therefore, we replace g(⋅) with a convolution neural network G k . The intuition is that the convolutional layer has powerful ability of feature representation that it can generalize to any regularization norms. Specifically, G k  takes as input k-space data of a shape 256×256 and outputs feature vectors of the same shape. It consists of three convolution layers and each uses 3×3 kernels with 32 channels, followed by a Batch Normalization layer and a PReLU activation layer. The advantage of G k  over g(⋅) is that convolution neural networks can be generalized to any form of regularization term and hence derivation of closed form of g(⋅) is no longer needed. The second modification is that we substitute L k  with a fully connected layer called difference layer D k  since L k  is not easy to compute under some rigorous constraints in backtracking FISTA. Furthermore, F T  (E(Fv k )−y) is a constant value and we name it as C k . The outputs of G k , D k  and C k  are combined through element-wise addition and connected to another fully connected layer F k . Finally, equation (5) is transformed to:
 
 x   k   =z   n−1   −F   k {( G   k−1 ( z   n−1 )+ D   k ( z   n−1   −v   k )+ C   k ( v   k )}  (6)
 
     Table 1, below, shows the average reconstruction performance of 3-fold cross validation performed by the process  16   a : 
     
       
         
           
               
               
               
               
               
             
               
                   
                 TABLE 1 
               
             
            
               
                   
                   
               
               
                   
                 RMSE↓ 
                 PSNR↑ 
                 RMSE↓ 
                 PSNR↑ 
               
            
           
           
               
               
               
            
               
                 Method 
                 50% 
                 25% 
               
               
                   
               
            
           
           
               
               
               
               
               
            
               
                 L1-wavelet ┌14┐ 
                 0.0680□0.0 
                 23.5844□2.0 
                 0.0975□0.0 
                 20.6242□2.7 
               
               
                 Total variation 
                 0.0713□0.0 
                 23.1699□2.0 
                 0.0908□0.0 
                 21.0009□1.6 
               
               
                 Low-rank ┌10┐ 
                 0.0508□0.0 
                 26.2022□2.4 
                 0.1014□0.0 
                 20.2634□2.7 
               
               
                 RU-Net 
                 0.0331□0.0 
                 29.1018□1.9 
                 0.0546□0.0 
                 24.0784□1.4 
               
               
                 FR-Net-2D 
                 0.0225□0.0 
                 32.8431□2.3 
                 0.0430□0.0 
                 26.5066□1.6 
               
               
                 FR-Net- 
                   0.0220 □ 0.0   
                   33.0129 □ 2.2   
                   0.0411 □ 0.0   
                   26.9550 □ 1.7   
               
               
                 Joint-FR-2D 
                 0.0236□0.0 
                 32.3520□2.3 
                 — 
                 — 
               
               
                 Joint-FR- 
                 0.0241□0.0 
                 32.3455□2.1 
                 — 
                 — 
               
               
                   
               
            
           
         
       
     
       FIGS.  3 A- 3 B  are graphs illustrating reconstruction results with k-space data at undersampled rate 25%. 
     Conventional MRI segmentation approaches consist of two separate steps. Basically, the first step is to reconstruct an image from undersampled k-space data and the second step is to feed the reconstructed image into an established segmentation model using automatic methods. For simplicity, we call this conventional approach Two-step Model. We propose a joint model called Joint-FR-Net. As shown in  FIG.  1   , Joint-FR-Net is an end-to-end deep neural network, which takes k-space data as input and directly learns segmentation mask as output. One significant difference from Two-step Model is that the “reconstructed image” is now a set of differentiable parameters that connect FR-Net  24  and U-Net  28 . The motivation behind this change is to bridge the gap between two originally-isolated models to become capable of learning parameters through backpropagation of shared gradients in the whole pipeline. From the perspective of model effectiveness and training efficiency, this mechanism allows two models to mutually benefit from each other. FR-Net  26  passes to the input of U-Net  28  the on-the-fly features of the “reconstructed image” instead of fixed noisy images so that features are shared and learned by both models. Meanwhile, U-Net  28  directly backpropagates gradients to FR-Net  24  to make the whole model more segmentation-oriented. 
     Technically speaking, let Re be the reconstruction network parameterized by θ and S Φ  be the segmentation network parameterized by Φ. In our application, R θ  and S Φ  respectively represent FR-Net  24  and U-Net  28 . Let y denote the input of k-space data. The reconstructed image is R θ (y) and the segmentation result is S Φ (R θ (y)). Note that the input of the segmentation network is the output of the reconstruction network and hence the parameters and θ are updated simultaneously. 
     We also define a combination loss    com  which is based on both the reconstruction loss function I R  and the segmentation loss function l S .    com  is defined as: £ com =I R (R θ (y), {tilde over (x)})+β·l S (S Φ (R θ (y)), {tilde over (x)}), where {tilde over (x)} is the ground truth of the reconstruction image and {tilde over (x)} is the ground truth of the reconstruction mask. B is the hyper-parameter that balances between reconstruction and segmentation performances. Reconstruction loss l R  is defined as l 1  loss and segmentation loss l s  is defined as dice loss. During the training procedure, the key issue is to find an optimal equilibrium between reconstruction and segmentation such that they are mutually beneficial. 
     The process  16   a  was experimentally tested. All experiments are conducted on a dynamic cardiac MRI dataset, which is collected from 25 volunteers and three of them are diagnosed with cardiac dyssynchrony disease. It contains 4, 000 2D slices of SAX images with manually annotated LV contours over various spatial locations and cardiac phases. We randomly sample ⅔ of the 25 subjects as the training set and the rest ⅓ as the test set. We apply Fourier transform to convert SAX image to k-space data and undersample it. Particularly, we keep eight lowest spatial frequencies and adopt Cartesian un-dersampling on k-space data along phase direction. Then undersampled k-space data can be simulated. Our model is assessed on the data at undersampled rates of 25% and 50%. Since the dataset is a dynamic sequence of SAX images, we consider both the 2D case and 2D+time case in the following experiments. In the reconstruction task, we quantitatively evaluate our method by root-mean-square error (RMSE) and peak signal-to-noise ratio (PSNR). In the segmentation task, we choose the commonly-used Dice&#39;s score and average perpendicular distance (APD) as our measure. 
     As shown in  FIG.  2   , we initialize x 0  by Fourier transform on the k-space data. The number of steps (k) and the number of blocks between each step (n) are unknown and there could be infinite combinations. Therefore, we run an empirical experiment to decide the structure of the FR-Net.  FIGS.  3 A- 3 B  show the corresponding RMSE and training time of different structures. We train different structures with undersampled data at rate of 50% and only run 15 epochs. Specifically, the line k=1 means there is only one step, x 0  and x 1 . The line k=2 means there are two steps x 0 , x i  and x 2 . We increase the number of blocks between each step and see how their RMSE changes. Balancing between time efficiency and reconstruction quality, we choose five blocks respectively between x 0  and x 1 , x 1  and x 2  (k=2). X 2  is the output of reconstructed image. 
     We compare our model with both CS-based methods and the state-of-the-art deep learning approach. Specifically, we consider following three classic CS-based methods with different regularization terms:    1 -wavelet, the total variation approach, and the low-rank approach. We also consider one deep learning approach called RU-Net that is a variation of U-Net. Both FR-Net and RU-Net are implemented in PyTorch and trained with ADAM optimizer on Tesla K80. In order to further improve the performance, we also adopt a k-space correction strategy by replacing zeros of the original k-space data with the values of the reconstructed image that in k-space. The average results are reported in Table 1, above. As we can see, our FR-Net-2D model achieves the best performance. And it also outperforms RU-Net in two cases with the lowest RMSE 0.0225 and 0.0430 respectively. For 2D+time data, we use 5 successive frames as input. The results show that the FR-Net-2D+time model which performs on dynamic data, achieves the best performance in terms of RMSE and PSNR in both cases. This validates that our model performs very well at different undersampled rates with and without temporal information. 
       FIG.  4    shows the reconstruction results at different undersampled rates. As we can see, the CS-based methods miss lots of details and the reconstructed images are blurry and fuzzy. In contrast, deep learning methods are more precise and efficient. Specifically, our FR-Net 2D and FR-Net 2D+time models yield much clearer and high-quality images. 
     Table 2, below, shows average endo- and epicardium segmentation performance of test set. The input k-space data is at undersampled rate 50%. 
     
       
         
           
               
               
               
               
               
             
               
                   
                 TABLE 2 
               
             
            
               
                   
                   
               
               
                   
                 Dice↑ 
                 ADP(mm)↓ 
                 Dice↑ 
                 ADP(mm)↓ 
               
            
           
           
               
               
               
            
               
                 Method 
                 Endo. 
                 Epi. 
               
               
                   
               
            
           
           
               
               
               
               
               
            
               
                 Fully-sampled 
                 0.7661 ± 0.2658 
                 4.1774 ± 5.0946 
                 0.7754 ± 0.2677 
                 5.7527 ± 10.8612 
               
               
                 Zero-filling 
                 0.6357 ± 0.3009 
                  7.7575 ± 10.9018 
                 0.6570 ± 0.3040 
                 9.8359 ± 15.0050 
               
               
                 Two-step 
                 0.6635 ± 0.3081 
                  9.1255 ± 21.1796 
                 0.7494 ± 0.2734 
                 8.2817 ± 18.4909 
               
               
                 Joint-FR-2D 
                 0.7260 ± 0.2787 
                 4.9771 ± 5.7992 
                 0.7503 ± 0.2922 
                 6.9593 ± 12.3803 
               
               
                 Joint-FR-2D + time 
                 
                   0.7310 ± 0.2603 
                 
                 
                   4.2884 ± 4.9119 
                 
                 
                   0.7573 ± 0.2763 
                 
                 
                   1.2628 ± 1.6988  
                 
               
               
                   
               
            
           
         
       
     
       FIG.  5    shows segmentation results of input data is at undersampled rate 50%. In this experiment, we train a U-Net with fully-sampled images and use it as the best achievable model. We also train another model with zero-filling images as the lower bound of segmentation performance. In order to show the mutual benefits of the Joint-FR-Net on solving the myocardium segmentation problem. We compare our method with the Two-step model which takes as input the images reconstructed from the FR-Net. Our method respectively takes as input the 2D k-space data and the 2D+time dynamic data and we call them Joint-FR-2D and Joint-FR-2D+time in this experiment. The training and test data are based on the undersampled rate of 50%. Experiments are conducted on both endo- and epicardium segmentation. The average segmentation results are reported in Table 2 above. We observe that our two joint models outperform the Two-step model in terms of Dice and APD for both endo- and epicardium cases. The Two-step model only gains 0.6635 and 0.7494 Dice score in endo. and epi. case, respectively. We also find that the Joint-FR-2D+time model achieves better results than the Joint-FR-2D model. This highlights the generalization ability of our model to dynamic data. Moreover, we can see the Dice index of Joint-FR-2D+time model is very close to that of fully-sampled model. The reconstruction performance of the Joint-FR-2D and Joint-FR-2D+time models are also reported in Table 1. It achieves comparable results with the FR-Net-2D and FR-Net-2D+time models and outperforms other CS-Based methods. Thus, we conclude that our joint model not only benefits the segmentation task but also achieves promising reconstruction results compared with the other CS-based models. We also plot the segmentation results in  FIG.  5   . Note that for better visualization, we use fully-sampled image as reconstructed image here. As the Joint-FR-Net models are trained with the raw k-space data, they compute similar contours to the ground truth. 
       FIG.  6    is a diagram illustrating processing performed by step  16   b  of  FIG.  1   , in greater detail. By way of background, in Compressed Sensing (CS)-based MRI reconstruction problems, the goal is to find a reconstructed image {tilde over (x)} such that it minimizes the reconstruction error between original k-space data y and Fourier transform of x: 
     
       
         
           
             
               
                 
                   
                     x 
                     ~ 
                   
                   = 
                   
                     
                       arg 
                       ⁢ 
                       
                         
                           min 
                           x 
                         
                         
                           
                              
                             
                               
                                 
                                   F 
                                   u 
                                 
                                 ( 
                                 x 
                                 ) 
                               
                               - 
                               y 
                             
                              
                           
                           2 
                           2 
                         
                       
                     
                     + 
                     
                       
                         R 
                         ⁡ 
                         ( 
                         x 
                         ) 
                       
                       . 
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
     Here F u  is an operator that transforms x into Fourier domain with undersampling. R(⋅) is a regularization term that depends on the prior information of data and general choices are    1  or    0 . 
     Traditionally, the objective function of Equation 1 immediately above is solved in an iterative manner that requires thousands of iteration steps until convergence. A more efficient way is to approximate the optimal reconstructed image from undersampled k-space data via deep neural network. One issue is that as the network goes deeper, the original information may be degraded. A Data Consistency (DC) layer can be used to avoid loss of information. Basically, DC layer takes as input the reconstructed image x n  from the n-th reconstruction block and outputs an updated “reconstructed” image x dc   n . Formally, DC layer is defined as
 
 x   n   dc   =DC ( x   n )= F   H (τ( F ( x   n ))).  (2)
 
     Here F(x) is Fourier transform that takes image x as input and outputs {circumflex over (x)} in Fourier domain and F H (x) is inverse Fourier transform. τ(x) is the data fidelity operation whose output has the same dimension as {circumflex over (x)}: 
                       τ   ⁡   (     x   ^     )     [     i   ,   j     ]     =     {               x   ^     [     i   ,   j     ]     ,             (     i   ,   j     )     ∉   Ω                       x   ^     [     i   ,   j     ]     +     vy   [     i   ,   j     ]         1   +   v       ,             (     i   ,   j     )     ∈   Ω                     (   3   )               
where [i,j] is the matrix indexing operation, Ω is the set of sampled positions of k-space data and v∈[0, ∞) is the noise level. In the noiseless case (v→∞), we have {circumflex over (x)}[i,j]=y[i,j] if (i,j)∈Ω, i.e. filling {circumflex over (x)} with the original values of k-space data y at position (i,j).
 
     The process  16   b  is shown in greater detail in  FIG.  6    and is referred to herein as MR Cascaded Channel-wise Attention Network (MICCAN). As shown in  FIG.  6   , MICCAN mainly consists of two parts: U-net with Channel-wise Attention (UCA) modules  40 ,  42  and DC layer  50 . These components are cascadedly coupled together and repeat for N times. Formally, denote the n-th UCA module  40  and the n-th DC layer  50  respectively by UCA n  and DC n . The starting point of MICCAN is undersampled k-space data y (elements  46   a ,  46   b ) which is later converted into a zero-filling image (element  48 ) x 0 =F H (y) through inverse Fourier transform F H  and fed to a UCA module  40 . The reconstructed image  52  is produced by the module  40 . Our cascaded model can be simply formalized as 
                   {                   x   n     =       U   ⁢   C   ⁢       A   n     (     x     n   -   1     dc     )       +     x     n   -   1     dc                     x   n   dc     =       DC   n         (     x   n     )                   (       n   =   1     ,   …         ,   N     )                   (   4   )               
where x 0   dc  is initialized as x 0 . The final reconstructed image of MICCAN, namely x N   dc , is produced by the last DC layer (DC N ).
 
     In the previous work on reconstruction problem, deep learning based methods have two major issues. First, they treat each channel-wise feature equally, but contributions to the reconstruction task vary from different feature maps. Second, receptive field in convolutional layer may cause to lose contextual information from original images, especially high-frequency components that contain valuable detailed information such as edges and texture. We developed the UCA module by introducing an attention mechanism that filters the useless features and enhance the informative ones. The attention technique is only applied in the decoder part. The intuition is that features of the decoder are extracted from coarse to fine feature-maps of multiple scales via skip connection. The attention module filters salient and prunes irrelevant and noisy features such that allows model parameters in shallower layers to be updated mostly that are relevant to a given task. 
     Specifically, we use global average pooling to extract the channel-wise global spatial information to vector z∈   C , whose c dimension is defined as 
                     z   c     =       1     H   ×   W       ⁢       ∑     i   =   1     H         ∑       j   =   1       W           f   c     [     i   ,   j     ]     ⁢     {       c   =   1     ,   …         ,   C     )                     (   5   )               
where f c ∈   W×H  is the feature map in the c-th channel. Such operation squeezes the spatial information of the whole image into a vector length of C. To further extract feature related to the final task, we introduce another gating mechanism as follows:
 
 {circumflex over (x)}   c =σ(δ( z*W   1 )* W   2 )⊙ f   c ,  (6)
 
where “*” is convolution operator and δ(⋅) is ReLU activation function to encode the channel-wise dependencies. W 1  is a kernel in the first convolutional layer that reduces the C-dimensional feature vector into C/r. On the contrary, kernel W 2  increases feature size back to C. Sigmoid function σ(⋅) is used to compute weighted attention map, which is later applied to rescaling the input feature f c . Based on this attention mechanism in the UCA module, our model MICCAN achieves very promising results and outperforms several state-of-the-art methods.
 
     The end-to-end network consists of UCA blocks  40 , Data Consistency blocks  50 , and Channel-wise Attention units  42 . It takes as input the undersampled k-space data  46   a  and the undersampled mask  46   b  (leftmost) and outputs the reconstructed MRI image  52  (right-most image). The zero-filled reconstruction image  48  (second-left image) works as the start point for reconstruction. The Data Consistency unit  50  employs the original k-space data  46   a  for further refinement. 
     To address the problem of vanishing low frequency in deep learning based MRI reconstruction, we utilize a long skip connection from the zero-filling image to the final reconstruction block. Specifically, we replace the residual in the last UCA module, namely X N =UCA N  (X dc   N−1 )+x dc   N−1 , with a long skip connection:
 
 x   N   =UCA   N ( x   N−1   dc )+ x   0   (7)
 
This simple modification is used to learn the global residual and to stabilize the gradient flow in deep residual network.
 
     A common choice of loss function for reconstruction problems is  2, but the resulting reconstructed image is of low quality and lacks high frequency detail. Therefore, we propose to use a combination of loss functions including    1  loss and perceptual loss    p . Given target image x s  and reconstructed image x=x dc   N  of MICCAN parameterized by θ, the combined loss is defined as 
                         ℓ   θ     (     x   ,     x   s       )     =         λ   1     ⁢       ℓ   1   θ     (     x   ,     x   s       )       +       λ   p     ⁢       ℓ   p   θ     (     x   ,     x   s       )           ⁢   
   where           (   8   )                               ℓ   1   θ     (     x   ,     x   s       )     =            x   -     x   s            1             (   9   )                               ℓ   p   θ     (     x   ,     x   s       )     =       ∑     k   =   1     K                  ϕ   VGG   k     (   x   )     -       ϕ   VGG   k     (     x   s     )            2   2               (   10   )               
where λ 1  and λ p  are weighing factors for two losses, ΦkV GG(⋅) represents features of the k-th activation layer in VGG network. Note that perceptual loss    p  minimizes the    2  distance between reconstruction image and target image in K different feature spaces, or equivalently it encourages the predicted image to be perceptually similar to the target image.
 
     To evaluate the effectiveness of process  16   b , we compare it with several state-of-the-art approaches on a simulated cardiac k-space dataset. We evaluate following methods in the experiment. Two traditional CS-MRI methods include  1-wavelet and TV norm, which are implemented in the BART toolbox. Two deep learning based models include DC-CNN and MRN5. The only difference between our proposed MICCAN and MRN5 is that MRN5 does not used attention module and long skip connection. We also consider three variants of our methods: MICCAN with    2  loss (MICCAN-A), MICCAN with combined loss (MICCAN-B), and MICCAN with both long skip connection and combined loss (MICCAN-C). We set the reduction ratio r as 8 for all our MICCAN models. For the combined loss, we set λ 1  as 10 and λ p  as 0.5. As the shallow features encode details of images, we compute the perceptual loss of layer relu1-2, relu2-1, relu2-2 and relu3-1 of the VGG-16 network. All deep learning models are implemented using PyTorch and trained on NVIDIA K80. Learning rate is initialized as 10 −4  with decreasing rate of 0.5 for every epochs. The training batch is 8 and the maximum number of epochs is 50. For fair comparison, we set the number of reconstruction blocks N as 5. 
       FIG.  7    illustrates wisualization of reconstructed images and reconstruction errors with k-space data at undersampled rate 12.5%. Also, Table 3, below, shows reconstruction results of the proposed MICCAN models and other methods. 
     
       
         
           
               
               
               
               
               
             
               
                   
                 TABLE 3 
               
               
                   
                   
               
               
                   
                 Methods 
                 NRMSE 
                 PSNR 
                 SSIM 
               
               
                   
                   
               
             
            
               
                   
                 TV 
                 0.1087 ± 0.0 
                 19.7388 ± 2.8 
                 0.5923 ± 0.0 
               
               
                   
                 i1-wavelet 
                 0.0753 ± 0.0 
                 22.7054 ± 2.0 
                 0.6333 ± 0.0 
               
               
                   
                 DC-CNN 
                 0.0587 ± 0.0 
                 24.7993 ± 1.7 
                 0.6612 ± 0.0 
               
               
                   
                 MRN5 
                 0.0427 ± 0.0 
                 27.5373 ± 1.5 
                 0.7851 ± 0.0 
               
               
                   
                 MICCAN 
                 0.0402 ± 0.0 
                 28.0664 ± 1.6 
                 0.8005 ± 0.0 
               
               
                   
                 MICCAN 
                 0.0391 ± 0.0 
                 28.3283 ± 1.7 
                 0.8214 ± 0.0 
               
               
                   
                 MICCAN 
                 0.0385 ± 0.0 
                 28.4489 ± 1.6 
                 0.8198 ± 0.0 
               
               
                   
                   
               
            
           
         
       
     
     A cardiac MRI dataset with 15 subjects is adopted in our experiments. We randomly choose 10 subjects as training set, 2 subjects as validation set and the rest 3 subjects as test set. We follow the k-space data simulation method. It assumes the sampled mask follows a zero-mean Gaussian distribution and the Cartesian undersampling method is adopted, also keeps the eight lowest spatial frequencies. Our model is evaluated on data with undersam-pled rate at 12.5% (acceleration rate 8×). 
     We also quantitatively evaluate all models with three widely used measurements: normalized root square mean error (NRMSE), peak signal-to-noise ratio (PSNR) and the structural similarity index measure (SSIM). The results are shown in Table 3 and we basically observe following three trends. First, all our MICCAN variants outperform other baseline models. This mainly attributes to the attention mechanism. Even trained with £2 loss, MICCAN-A achieves NRMSE of 0.0402, PSNR of 28.0664 and SSIM of 0.8005 that beats MRN5 model. Second, MICCAN-B model that is trained with combined loss gains better result than MICCAN-A. For example, MICCAN-B has 2.7% decrease in NRMSE compared with MICCAN-A and also achieves the best SSIM value of 0.8214. This indicates that our combined loss is better than the £ 2  loss, making image less blur and keeping the perceptual details. Third, with long skip connection, MICCAN-C further improves reconstruction performance with the lowest NRMSE of 0.0385 and the highest PSNR value of 28.4489. Overall, all these results demonstrate the effectiveness of the channel-wise attention modules and the proposed combined loss. 
     We further visualize the reconstructed image and reconstruction errors on a test sample and qualitatively analyze the results of all methods. As shown in  FIG.  7   , most of the baseline methods cannot completely recover detail of the image and suffer from severe blurring artifacts. In contrast, our three MICCAN methods eliminate most blurring artifacts and recover more details from low undersampled k-space data. Furthermore, for the transitional methods such as    1 -wavelet and TV, the reconstructed images are similar to the zero-filling image and suffer from heavy blurring artifacts. We highlight that our MICCAN achieves best results and with much fewer reconstruction errors. 
       FIGS.  8 - 9    are diagrams illustrating processing performed by step  16   c  of  FIG.  1   , in greater detail. The following framework is provided in connection with step  16   c . Given a sequence of under-sampled k-space data {yt} t∈[T]  of T frames, the dynamic MRI reconstruction problem is to predict reconstructed images (zt) t∈[T]  from {y t }, which can be formalized as an optimization problem: argmin {zt}  L rec ({z t }, where: 
     
       
         
           
             
               
                 
                   
                     
                       ℒ 
                       rec 
                     
                     ( 
                     
                       { 
                       
                         z 
                         t 
                       
                       } 
                     
                     ) 
                   
                   = 
                   
                     
                       
                         ∑ 
                         
                           t 
                           = 
                           1 
                         
                         T 
                       
                       
                         
                           1 
                           2 
                         
                         ⁢ 
                         
                           
                              
                             
                               
                                 
                                   F 
                                   u 
                                 
                                 ( 
                                 
                                   z 
                                   t 
                                 
                                 ) 
                               
                               - 
                               
                                 y 
                                 t 
                               
                             
                              
                           
                           2 
                           2 
                         
                       
                     
                     + 
                     
                       λ 
                       ⁢ 
                       
                         
                           R 
                           ⁡ 
                           ( 
                           
                             z 
                             t 
                           
                           ) 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   1 
                   ) 
                 
               
             
           
         
       
     
     The term ∥F u (z t )−y t ∥ 2   2  is used to guarantee data consistency by restricting the reconstructed image z, to be close to the input measurement y t . F u (⋅) is an operator that transforms image-domain zt into Fourier domain followed by undersampling. R(⋅) is a regularization function that depends on the prior knowledge of the input {y t }. Common choices include sparsity in transformed domain, total variation (TV) penalties, and low-rank property. λ is a weighting factor. 
     In order to capture anatomical motion in the dynamic MRI acquisition, it is natural to incorporate motion estimation/motion compensation (ME/MC) technique in the reconstruction process. Specifically, based on the brightness constancy assumption, for a temporal 2D image z t (x, y, t) with small movement (Δx, Δy, Δt) with respect to the next frame, we add the following motion estimation constraint to the objective function (1): 
                         ℒ   me     (     {     v   t     }     )     =         ∑     t   =   1       T   -   1                    ∇     z   t   T       ⁢     v   t       +       ∂     z   t         ∂   t              1       +     δ   ⁢            v   t          1           ,           (   2   )               
where
 
               ∇       z   t     (     x   ,   y     )       =     (         ∂     z   t         ∂   x       ,       ∂     z   t         ∂   y         )           
are the derivatives of image z t  at position (x,y), and
 
                 v   t     (     x   ,   y     )     =     (         Δ   ⁢   x       Δ   ⁢   t       ,       Δ   ⁢   y       Δ   ⁢   t         )           
is the estimated displacement motion fields or optical flow.
 
     Furthermore, given the estimated motion field v t , the reconstructed image z t  can be re-fined through MC process, i.e. c t =MC(z t , z 1 , z T )+r t , where c t  is the motion-compensated reconstructed image and r t  is a residual term for better exploiting temporal redundancy. Therefore, we can derive the motion compensation constraint as follows. 
     
       
         
           
             
               
                 
                   
                     
                       ℒ 
                       mc 
                     
                     ( 
                     
                       { 
                       
                         r 
                         t 
                       
                       } 
                     
                     ) 
                   
                   = 
                   
                     
                       ∑ 
                       
                         t 
                         = 
                         1 
                       
                       
                         T 
                         - 
                         1 
                       
                     
                     
                       
                         1 
                         2 
                       
                       ⁢ 
                       
                         
                           
                              
                             
                               
                                 
                                   F 
                                   u 
                                 
                                 ( 
                                 
                                   c 
                                   t 
                                 
                                 ) 
                               
                               - 
                               
                                 y 
                                 t 
                               
                             
                              
                           
                           2 
                           2 
                         
                         . 
                       
                     
                   
                 
               
               
                 
                   ( 
                   3 
                   ) 
                 
               
             
           
         
       
     
     By combining with two motion-based constraints of Equations (2) and (3) immediately above, the motion-guided dynamic MRI reconstruction problem is defined as: 
     
       
         
           
             
               
                 
                   
                     
                       
                         arg 
                         ⁢ 
                         min 
                       
                       
                         { 
                         
                           
                             z 
                             t 
                           
                           , 
                           
                             v 
                             t 
                           
                           , 
                           
                             r 
                             t 
                           
                         
                         } 
                       
                     
                     ⁢ 
                     
                       
                         ℒ 
                         rec 
                       
                       ( 
                       
                         { 
                         
                           z 
                           t 
                         
                         } 
                       
                       ) 
                     
                   
                   + 
                   
                     
                       ηℒ 
                       me 
                     
                     ( 
                     
                       { 
                       
                         v 
                         t 
                       
                       } 
                     
                     ) 
                   
                   + 
                   
                     
                       
                         ζℒ 
                         mc 
                       
                       ( 
                       
                         { 
                         
                           r 
                           t 
                         
                         } 
                       
                       ) 
                     
                     . 
                   
                 
               
               
                 
                   ( 
                   4 
                   ) 
                 
               
             
           
         
       
     
     The solution to Equation (4) immediately above is non-trivial and traditional CS-based algorithms are usually computationally expensive and require long running time for hyper-parameter tuning. Recent advances in deep learning provide an alternative way for efficient MRI reconstruction, but very few works focused on the dynamic reconstruction problem and they only targeted for the simpler Equation (1) immediately above without considering motion information. To this end, we propose a deep learning based method called Motion-guided Dynamic Reconstruction Network (“MODRN”) to solve Equation (4) immediately above. 
     Process  16   c  separates the motion-guided dynamic reconstruction problem into three closely-connected parts: (i) Dynamic Reconstruction (DR) component for estimating initial reconstructed image from Equation (1) immediately above; (ii) Motion Estimation (ME) component for generating motion information through Equation (2) immediately above; and (iii) Motion Compensation (MC) component for refining reconstructed image guided by learned motion based on Equation (3) immediately above. 
     Instead of directly solving Equation (1) immediately above, an iterative process is adopted through DR component to approximate reconstructed images z t . Formally, given under-sampled k-space data {y t}t∈[T]  with sampled mask Ω, DR component learns to reconstruct images in N iterations: 
     
       
         
           
             
               
                 
                   
                     z 
                     t 
                     
                       ( 
                       n 
                       ) 
                     
                   
                   = 
                   
                     
                       DR 
                       ⁡ 
                       ( 
                       
                         
                           z 
                           t 
                           
                             ( 
                             
                               n 
                               - 
                               1 
                             
                             ) 
                           
                         
                         , 
                         
                           y 
                           t 
                         
                         , 
                         Ω 
                       
                       ) 
                     
                     ⇔ 
                     
                       { 
                       
                         
                           
                             
                               
                                 
                                   x 
                                   bt 
                                   
                                     ( 
                                     n 
                                     ) 
                                   
                                 
                                 , 
                                 
                                   
                                     b 
                                     t 
                                     
                                       ( 
                                       n 
                                       ) 
                                     
                                   
                                   = 
                                   
                                     B 
                                     ⁡ 
                                     ( 
                                     
                                       
                                         
                                           f 
                                           enc 
                                         
                                         ( 
                                         
                                           z 
                                           t 
                                           
                                             ( 
                                             
                                               n 
                                               - 
                                               1 
                                             
                                             ) 
                                           
                                         
                                         ) 
                                       
                                       , 
                                       
                                         b 
                                         t 
                                         
                                           ( 
                                           
                                             n 
                                             - 
                                             1 
                                           
                                           ) 
                                         
                                       
                                     
                                     ) 
                                   
                                 
                               
                             
                           
                           
                             
                               
                                 
                                   x 
                                   st 
                                   
                                     ( 
                                     n 
                                     ) 
                                   
                                 
                                 , 
                                 
                                   
                                     s 
                                     t 
                                     
                                       ( 
                                       n 
                                       ) 
                                     
                                   
                                   = 
                                   
                                     S 
                                     ⁡ 
                                     ( 
                                     
                                       
                                         
                                           f 
                                           
                                             dec 
                                             ⁢ 
                                             1 
                                           
                                         
                                         ( 
                                         
                                           x 
                                           bt 
                                           
                                             ( 
                                             n 
                                             ) 
                                           
                                         
                                         ) 
                                       
                                       , 
                                       
                                         s 
                                         t 
                                         
                                           ( 
                                           
                                             n 
                                             - 
                                             1 
                                           
                                           ) 
                                         
                                       
                                     
                                     ) 
                                   
                                 
                               
                             
                           
                           
                             
                               
                                 
                                   z 
                                   t 
                                   
                                     ( 
                                     n 
                                     ) 
                                   
                                 
                                 = 
                                 
                                   DC 
                                   ( 
                                   
                                     
                                       
                                         
                                           f 
                                           
                                             dec 
                                             ⁢ 
                                             2 
                                           
                                         
                                         ( 
                                         
                                           x 
                                           st 
                                           
                                             ( 
                                             n 
                                             ) 
                                           
                                         
                                         ) 
                                       
                                       + 
                                       
                                         z 
                                         t 
                                         
                                           ( 
                                           
                                             n 
                                             - 
                                             1 
                                           
                                           ) 
                                         
                                       
                                     
                                     , 
                                     
                                       y 
                                       t 
                                     
                                     , 
                                     Ω 
                                   
                                   ) 
                                 
                               
                             
                           
                         
                         , 
                         
 
                         	 
                         
                           n 
                           ∈ 
                           
                             
                               [ 
                               N 
                               ] 
                             
                             . 
                           
                         
                       
                     
                   
                 
               
               
                 
                   ( 
                   5 
                   ) 
                 
               
             
           
         
       
     
     where z t   (0)  is zero-filling image and z t   (n)  is the reconstructed image of y t  after iteration n. B and S are two ConvGRU units that respectively output features x bt   (n)  and x st   (n)  together with hidden states b t   (n)  and s t   (n) . f enc  and f dec1 , f dec2  are convolutional encoder and decoders in the U-Net, which is used as the backbone of the DR component to capture course-to-fine features of reconstructed images. Equation (5) immediately above is visualized in  FIG.  8    for better understanding. 
     As can be seen in  FIG.  8   , process  16   c  includes input data  60 , recurrent processing steps  62 , under-sampled k-space data  64 , intermediate data  66 , and final processing steps  68 . One benefit here is that regularization function R(⋅) in Equation (1) immediately above is now built upon the convolutional network for automated feature learning and hence avoid the requirements of prior knowledge on the selection of R. DC(⋅) is the differentiable DC layer that takes the same effect as the data consistency term ∥Fu(z t )−y t ∥ 2   2  in Equation (1) immediately above to force the reconstructed image to be consistent with the input data. It fills the reconstructed image z t , with the original values of input data y t  in the Fourier domain by the sampled mask Ω. 
     In order to capture dynamic information of image sequence during each iteration, two kinds of ConvGRU units are introduced in recurrent processing steps  62 , namely, B and S, inspired by Equation (5) immediately above. The difference between B and S is that GRU unit S is used to improve the performance of image z t  over N iterations while the role of B is to connect dynamic information of neighboring images z t−1  and z t , which is implemented by initializing hidden state b t   (0)  as b t−1   (N) . Finally, we impose l 1  on the reconstructed images {z t   N } with respect to ground truth for penalizing 
     In analogy to Equation (2) immediately above, the Motion Estimation (ME) component takes as input the sequence of reconstructed images {zt} t∈[T]  and learn to predict displacement motion fields {vt} t∈[T] . As shown in  FIG.  9   , our proposed ME component embraces two parts: a FlowNet backboned by convolutional U-Net (U-FlowNet)  70  for motion field estimation. The other is a differentiable sampling layer  72  based on Spatial Transformer Network (STN), which endows convolutional network with the ability to warp the spatial deformation between images. Unlike traditional optimization algorithms for motion estimation that depend on a strong assumption that the brightness of two fames should be consistent and the movement of the foreground object is small, our method does not succumb to any assumption and hence is more applicable in practical dynamic MRI reconstruction. The performance of ME is heavily affected by noisy input, therefore it is pre-trained with two fully sampled images z t−1  and z t . The image pair is first fed to the U-FowNet  70 , which produces two-channel displacement v t−1  along the x and y directions. Then, the sampling layer warps z t−1  towards z t  by using v t−1  and yields a warping image denoted by z t   0  through differentiable bilinear interpolation. This leads to a natural re-formalization of motion estimation (ME) loss    me  between z t−1  and z t  from Equation (2):
 
   me ( z   t−1   ,z   t )=∥ z′   t   −z   t ∥ 1   +β∥v   t−1 ∥ 1   +γ∥v   t−1 ∥ TV .  (6)
 
The first term is an image reconstruction loss used to keep the majority of high-frequency parts on images. Two additional regularization terms reinforce constraints on the motion field v t−1 , where l 1  regularization is to suppress unreal large magnitude of displacement and total-variation (TV) regularization is to make the displacement locally smooth.
 
     In addition, the above loss only enforces temporal consistency between consecutive frames, but there is no guarantee for long-term coherence. Therefore, we consider to train the U-FlowNet  70  with three sets of ME losses to capture long-term motion information, as illustrated in  FIG.  9   . 
                         ℒ   me     (     {     z   t     }     )     =           ∑     T   -   1         t   =   2           ℓ   me     (       z   1     ,     z   2       )       +       ∑     t   =   2       T   -   1           ℓ   me     (       z   t     ,     z   T       )       +       ∑     t   =   2       T   -   2           ℓ   me     (       z   t     ,     z     t   +   1         )           ,           (   7   )               
where three terms on the right-hand-side are respectively forward ME loss, backward ME loss and neighboring ME loss.
 
     Motion Compensated (MC) component is used to refine reconstructed images {zt} t∈[T]  through motion information and to generate motion compensated image {ct} t∈[T] . During the MC stage  74  shown in  FIG.  9   , two additional fully sampled reference frames are added to learn more accurate displacement motion fields. The pre-trained U-FlowNet  70  is fixed and directly used as an operator in the MC stage  74 . As shown in  FIG.  9   , the MC stage  74  takes a reconstructed image z t  from the DR component and two reference frame z 1  and z T  as input. It first retrieves two warping images z′ t  and z″ t  from the ME stage  74  by feeding z 1 , z t  and z t , z T  respectively. These two images represent forward and backward motion information, which is then concatenated and fed to a residual network to generate residual information r t , as described in Equation (3) immediately above. Finally, the reconstructed image z t  together with the residual r t  are summed up to generate the motion-guided refined image c t , which is penalized by l 1  loss with respect to the ground truth image. 
     An experiment was conducted with a short-axis (SAX) cardiac dataset composed of 15 patients. Each subject contains around 12 SAX planes and each plane includes 24 phases (2D images) that form a whole cardiac cycle. The image resolution is normalized to 1.25 mm and image size is cropped to 152×152 pixels. In order to simulate k-space data, we adopt a Cartesian under-sampling method which assumes that sampled mask Ω follows a zero-mean Gaussian distribution and keeps 8 center spatial frequencies. We consider two different settings on the dataset respectively with under-sampling rates of 20% (or acceleration rate 5×) and 12.5% (8×). For convenience, we refer to these two cases as Rate 5× and Rate 8×. We perform 3-fold cross-validation in the following experiments that each fold contains 10 training subjects and 5 test subjects. Table 4, below, shows the average performance of dynamic MRI reconstruction on the test subjects in both cases of Rate 5× and Rate 8×. The best results are highlighted in bold font. 
     
       
         
           
               
               
               
               
               
               
               
             
               
                   
                 TABLE 4 
               
             
            
               
                   
                   
               
               
                   
                 NRMSE↓ 
                 PSNR↑ 
                 SSIM↑ 
                 NRMSE↓ 
                 PSNR↑ 
                 SSIM↑ 
               
            
           
           
               
               
               
            
               
                 Method 
                 5× 
                 8× 
               
               
                   
               
            
           
           
               
               
               
               
               
               
               
            
               
                 k-t SLR 
                 0.0934 
                 21.0858 
                 0.6794 
                 0.1054 
                 19.9504 
                 0.6193 
               
               
                 k-t FOCUSS 
                 0.0766 
                 22.7471 
                 0.6581 
                 0.0879 
                 21.4063 
                 0.5920 
               
               
                 k-t FOCUSS + ME/MC 
                 0.0758 
                 22.8139 
                 0.6701 
                 0.0854 
                 21.6547 
                 0.6131 
               
               
                 DC-CNN(3D) 
                 0.0360 
                 29.1292 
                 0.8449 
                 0.0513 
                 25.9709 
                 0.7441 
               
               
                 DRN w/o GRU 
                 0.0381 
                 28.7187 
                 0.8286 
                 0.0519 
                 25.9120 
                 0.7448 
               
               
                 DRN 
                 0.0349 
                 29.5394 
                 0.8502 
                 0.0485 
                 26.5275 
                 0.7687 
               
               
                 MODRN 
                 
                   0.0274 
                 
                 
                   32.0403 
                 
                 
                   0.9104 
                 
                 
                   0.0364 
                 
                 
                   29.4774 
                 
                 
                   0.8702 
                 
               
               
                   
               
            
           
         
       
     
     All of the deep learning models with were implemented with PyTorch and trained on an NVIDIA K80 system. All models are trained for total 80 epochs using Adam optimizer, with initialized learning rate of 5×10 −4  and decreasing rate of 0.5 for every 20 epochs. Due to hardware limitations, the number of iterations is set to be N=3 and the length of image sequence is T=12. 
     The dynamic reconstruction performance of the proposed methods was evaluated quantitatively and qualitatively in both cases of Rate 5× and Rate 8×. We consider three variants of our models: DRN w/o GRU (the one without GRU hidden unit), DRN (the one with DR component only) and MODRN (the complete version). We compare with four state-of-the-art approaches including k−t SLR, k−t FOCUSS, k−t FOCUSS+ME/MC, and DC-CNN (3D). The first three are traditional CS-based methods and only k−t FOCUSS+ME/MC includes ME/MC procedures. The last one is also a deep learning based method that explores spatio-temporal information using 3D convolution. Three common quantitative metrics are used: root square mean error (NRMSE), peak signal-to-noise ration (PSNR) and structural similarity index measure (SSIM). 
     The results of all methods are reported in Table 4 above. We observe that all our methods consistently outperform four state-of-the-art approaches in both Rate 5× and Rate 8× cases. In particular, MODRN achieves the best performance for all metrics, mainly attributing to the motion information exploited by ME/MC components. We also find that DRN outperforms DRN w/o GRU by a large margin, which indicates the importance of utilizing dynamic sequence of image. 
     To further investigate the performance of four deep learning methods, we plotted NRMSE values within a complete cardiac cycle of one example in  FIGS.  10 A- 10 B . It shows that our method MODRN consistently achieves the smallest error of dynamic reconstruction for the sequence of images. In contrast, the models without ME/MC are unstable along the temporal dimension, especially in the case of DC-CNN(3D). For example, in the case of Rate 8×, the gap between DRN and MODRN model become larger, which implies the significance of using motion information. 
     We visualize the reconstructed images and error with respect to ground truth of all methods in  FIG.  11   . It is obvious that all CS-based methods have streaking artifacts and larger reconstruction error while our MODRN model eliminates the most blurring artifacts and recovers more high-frequency details. 
     
       
         
           
               
               
               
               
             
               
                   
                 TABLE 5 
               
               
                   
                   
               
               
                   
                 Method 
                 Dice↑ 
                 HD↓ 
               
               
                   
                   
               
             
            
               
                   
                 Reference 
                 0.8130 
                 1.9254 
               
               
                   
                 Lucas-Kanade 
                 0.8125 
                 1.9577 
               
               
                   
                 U-FlowNet-A 
                 0.8297 
                 1.8755 
               
               
                   
                 U-FlowNet-B 
                 
                   0.8306 
                 
                 
                   1.8584 
                 
               
               
                   
                   
               
            
           
         
       
     
     The motion estimation results generated by the U-FlowNet were estimated from the ME component. Two baseline methods, Reference and Lucas-Kanade, are compared with our U-FlowNet-A (trained with only neighboring loss) and U-FowNet-B (trained with combined loss). Reference method directly calculates metrics using the segmentation of the target phase and the reference phase. Since it is impractical to obtain the ground truth of optical flow from cardiac MR, we compute the overlapped area of the myocardium between the targeting image and the warping image. In particular, we calculate the average Dice&#39;s score and Hausdorff Distance between Z 1  and other frames, Z T  and other frames and also neighboring frames. The results of 3-fold cross-validation are reported in Table 5 above. We observe that U-FowNet-B method achieves the best performance, which indicates that compared with neighboring loss, our combined loss contributes more to accurate motion estimation with large movement between frames. 
     Second, we compare the quality of motion-guided refined image by MC component of MODRN with that of reconstructed image by DRN alone. The results of three consecutive frames are visualized in  FIG.  12   . We can observe clear improvements of MODRN that its reconstruction error is reduced around cardiac region and no noticeable artifact is generated. 
       FIG.  13    is a diagram illustrating the motion-guided network discussed above in connection with process  16   c , implemented as a single, differentiable network. As can be seen, the network includes the DRN discussed above as the first layer of the differentiable network, which receives as input a zero-filling image sequence (a Fourier transform of under-sampled k-space sequences). The DRN outputs an initial reconstruction series of images. Two ME components are provided in a second layer of the differentiable network, and each processes both the initial reconstruction series of images and a reference frame and outputs a warped image. The warped images are then processed by the MC component (which forms the third layer of the differentiable network) to output a final reconstruction. The differentiable network shown in  FIG.  13    provides an end-to-end model that improves the performance of the MODRN network discussed above. Training in such an end-to-end manner is beneficial for allowing gradients from the MC component to be backpropagated to the ME components as well as to the DRN. Because of this, the ME components and the DRN are not only optimized for motion estimation and initial reconstruction, but also for generating final, motion-guided frames. 
     Having thus described the system and method in detail, it is to be understood that the foregoing description is not intended to limit the spirit or scope thereof. It will be understood that the embodiments of the present disclosure described herein are merely exemplary and that a person skilled in the art may make any variations and modification without departing from the spirit and scope of the disclosure. All such variations and modifications, including those discussed above, are intended to be included within the scope of the disclosure. What is desired to be protected by Letters Patent is set forth in the following claims.