Compressed Sensing: From Research to Clinical Practice with Data-Driven Learning

Compressed Sensing: F rom Researc h to Clinical Practice with Data-Driv en Learning Joseph Y. Cheng ∗ , F eiyu Chen, Christopher Sandino, Morteza Mardani, John M. P auly , Senior Memb er, IEEE, Shreyas S. V asanaw ala † Abstract Compressed sensing in MRI enables high subsampling factors while maintaining diagnostic image qualit y . This technique enables shortened scan durations and/or im- pro ved image resolution. F urther, compressed sensing can increase the diagnostic infor- mation and v alue from eac h scan p erformed. Ov erall, compressed sensing has significan t clinical impact in impro ving the diagnostic quality and patient exp erience for imaging exams. Ho wev er, a n umber of challenges exist when moving compressed sensing from researc h to the clinic. These c hallenges include hand-crafted image priors, sensitive tun- ing parameters, and long reconstruction times. Data-driv en learning pro vides a solution to address these challenges. As a result, compressed sensing can hav e greater clinical impact. In this tutorial, w e will review the compressed sensing formulation and outline steps needed to transform this form ulation to a deep learning framework. Supplemen- tary op en source co de in python will b e used to demonstrate this approach with op en databases. F urther, we will discuss considerations in applying data-driv en compressed sensing in the clinical setting. Keyw ords: compressed sensing, deep learning, clinical translation 1 In tro duction Compressed sensing is a p ow erful to ol in magnetic resonance imaging (MRI). As the scan duration is directly related to the num ber of data samples measured, collecting few er mea- suremen ts enables faster imaging. Multi-c hannel imaging, also kno wn as “parallel imaging,” lev erages the lo calized sensitivit y profiles of eac h receiv er element in a coil array to enable subsampling factors on the order of 2–6 [1 – 3]. By exploiting prop erties of the reconstructed images (suc h as the sparsit y of the W a velet co efficients of an image), compressed sensing can further increase the subsampling factors by 2–3 fold [4]. The additional factor is extremely p o werful in enabling a broad range of clinical applications. F or example, high-resolution v olumetric imaging that would take minutes to acquire can b e ac hieved in a single breath- hold. This strategy minimizes the sensitivity of patien t motion and increases the diagnostic qualit y of the resulting images. Ov er the past decade since compressed sensing was in tro duced to MRI [4], man y dev elop- men ts hav e b een made to extend this idea and bring it into clinical practice. One area with significan t research and clinical activit y is in m ulti-dimensional imaging. With more dimen- sions to exploit sparsit y , the subsampling factor can b e increased substantially (o ver 10 fold). ∗ J. Y. Cheng is with the Departmen t of Radiology , Stanford Universit y , Stanford, CA, 94022 USA e-mail: jyc heng@stanford.edu. † F. Chen, C. Sandino, M. Mardani, J. M. Pauly , and S. S. V asanaw ala are with Stanford Universit y . 1 As a result, multi-dimensional scans can b e completed in clinically feasible scan times. F or example, a v olumetric cardiac-resolved flow imaging sequence (4D flo w) can b e performed in a 5–10 min scan instead of an hour long scan needed to satisfy the Nyquist rate [5]. This single 4D flow scan enables a comprehensive cardiac ev aluation with flo w quantification, functional assessment, and anatomical information [6]. Instead of an hour long cardiac exam for congenital heart defect patients with complex cardiac anomalies, the exam can b e completed in a simple-to-execute single 4D flo w scan (Fig. 1A). Other examples of multi- dimensional compressed sensing include cardiac imaging [7, 8], dynamic-con trast-enhanced imaging [9, 10] (Fig. 1B), and “extra-dimensional” imaging with 4+ dimensions [11, 12]. t = 0 s 14. 5 s 29. 0 s 43. 5 s 58. 0 s A B Figure 1: Example clinical applications of compressed sensing. In A , tw o cardiac phases of the cardiac-resolved v olumetric v elo city MRI (4D flow) are shown for the purp ose of congenital heart defect ev aluation [6]. T o enable a clinically feasible scan duration of 5–15 min utes while main taining high spatial (0 . 9 × 0 . 9 × 1 . 6 mm 3 ) and high temp oral resolutions (22.0 ms), a subsampling factor of 15 w as used [5]. In B , high spatial (1 × 1 × 2 mm 3 ) with a 14.5-s temp oral resolution w as achiev ed for a dynamic con trast enhanced MRI using an acceleration factor of 6.5 [10]. High spatial-temp oral resolutions are required for capturing the rapid hemo dynamics of p ediatric patients. Additionally , rapid imaging has significant impact on the clinical w orkflow. Exam times can b e significantly shortened to reduce patient burden and discomfort. F or p ediatric imaging, the shortened exam time enables the reduction of the depth and length of anes- thesia [13]. F or extremely short scan times (less than 15 min utes), anesthesia can entirely eliminated. 1.1 Remaining Challenges Muc h success has b een observed b y applying compressed sensing to specific clinical appli- cations such as for cardiac cine imaging or for MR angiography . Ho wev er, the p otential that compressed sensing brings is not fully realized in terms of maximal acceleration and breadth of applications. A num b er of c hallenges limit the use of compressed sensing for 2 clinical practice. First, compressed sensing is sensitiv e to tuning parameters. One exam- ple of these parameters is the regularization parameter. These parameters determine the w eight of the regularization term with respect to the data consistency term in the cost function. Increasing the regularization parameter will improv e the p erceived signal-to-noise ratio (SNR) of reconstructed images. Ho wev er, while the p erceived SNR impro ves, fine structures in the reconstructed images may b e ov er-smo othed, or new image textures ma y b e in tro duced, resulting in an artificial appearance of images. The optimal v alue of the regularization parameter usually v aries among scans and patien ts, making parameter tun- ing clinically infeasible. The impact of the regularization parameter on the reconstruction p erformance is illustrated in Fig. 2. Figure 2: Example results of data-driv en compressed sensing (CS) with v ariational net works (VN) [14, 15]. Con ven tional compressed sensing requires tuning of the regularization param- eter λ . The optimal v alue of this parameter ma y v ary with differen t scans. Compared to con ven tional CS using L1-ESPIRiT [3] (left and right columns), learning a VN (middle col- umn) ac hieves prop er regularization without the need of manual tuning the regularization parameter. Second, compressed sensing may result in textural artifacts. Compressed sensing theory assumes that the unkno wn signal is sparse in some transform domain. This domain is often man ually chosen based on the characteristic of the image. F or MRI applications, common c hoices of this domain transformation are w av elet [4], total v ariation [8], and locally lo w- rank transforms [9]. How ever, these transforms ma y promote some textural artifacts in the reconstructed images. F or example, locally lo w-rank reconstruction ma y cause blo c k artifacts in reconstructed images. Lastly , compressed sensing usually has unknown reconstruction times b ecause of the use of iterativ e conv ex optimization algorithms. T o achiev e clinically acceptable image quality , the n umber of iterations for these algor ithms may v ary b et ween 50 to o ver 1000, resulting in v arying reconstruction times. F or v olumetric scans, these num b ers of iterations corresp ond to several minutes to o v er an hour long reconstruction times. This uncertaint y of long 3 reconstruction times may further lead to dela ys and queues in clinical scanning, as well as require exp ensiv e dedicated computational hardware. 1.2 Curren t Solutions Muc h researc h has attempted to solve these c hallenges. The performance of compressed sensing is highly dep enden t on the regularization function and regularization parameters used. T o reduce textural artifacts, domain transformations used in the regularization terms (i.e. W av elet transform, total v ariation) are carefully c hosen and tuned for eac h application. Impact of these different regularization functions must b e considered as differen t functions in tro duce different types of bias to the reconstructed images. Ideally , the regularization parameters should b e set based on the system noise. Ho w- ev er, man y other factors suc h as patient motion contribute to measuremen t error. Th us, regularization parameters are often manually tuned for eac h clinical application to achiev e acceptable p erformance. Automated parameter selection metho ds, such as Stein’s un biased risk estimate (SURE) [16], hav e b een developed to select the regularization parameters for compressed sensing MRI. How ever, this approach comes at the cost of additional reconstruc- tion times; as a result, this approach is uncommon in practice. In regards to reconstruction time, a trade-off betw een reconstruction time and p erformance can also be ac hiev ed b y early truncation of the iterative reconstruction [17]. Most of these solutions require manual adjustment of the original compressed sensing problem to sp ecific applications. This tuning and redesigning pro cess requires extra time and effort, and constant attention as hardw are and clinical proto cols evolv e. Previous solutions to automate some of this pro cess in volv e additional memory and more complex computations. 1.3 Ov erview T o ov ercome the obstacles with compressed sensing, data-driven learning has recently b e- come a compelling and practical approach [14, 18 – 21]. The purpose of this tutorial is to build the basic framew ork for extending compressed sensing to a data-driv en learning ap- proac h and to describe the considerations for clinical deplo yment. W e also discuss a num b er of new c hallenges when using this framew ork. W e hav e released supplementary python code on GitHub 1 to demonstrate an example of this data-driv en compressed sensing framework. 2 Data-Driv en Compressed Sensing 2.1 Bac kground The MRI reconstruction problem can be form ulated as a minimization problem [4]. The optimization consists of solving the following equation: ˆ m = arg min m 1 2 k Am − y k 2 2 + λR ( m ) (1) where m is the reconstructed image set, A describ es the imaging model, and y is the mea- sured data in the k-space domain. The imaging mo del for MRI consists of signal mo dulation 1 https://github.com/MRSRL/dl- cs 4 b y coil sensitivity profile maps S , F ourier transform op eration, and data subsampling. These sensitivit y profile maps S are sp ecific for each dataset y . The goal of this optimization is to reconstruct an image set m that b est matches the measured data y in the least squares sense. F or highly subsampled datasets, this problem is ill posed – man y solutions satisfy Eq. 1; th us, the regularization function R ( m ) and corresp onding regularization parameter λ incorp orate image priors to help constrain the problem. Man y optimization algorithms hav e been developed to solv e the minimization problem in Eq. 1. F or simplicit y , we base our discussion on the proximal gradient metho d. W e refer the reader to similar approac hes based on other optimization algorithms [14, 18, 19, 21]. T o solve Eq. 1, we split the problem in to tw o alternating steps that are rep eated. F or the k -th iteration, a gradien t up date is performed as m ( k +) = m ( k ) − 2 tA H ( Am ( k ) − y ) , (2) where A H is the transpose of the imaging mo del, and t is a scalar sp ecifying the size of the gradient step. The current guess of image m is denoted here as m ( k +) . Afterw ards, the pro ximal problem with regularization function R is solved: m ( k +1) = pro x λR  m ( k +)  = arg min u R ( u ) + 1 2 λ k u − m ( k +) k 2 2 , (3) where u is a help er v ariable that transforms the regularization into a con vex problem that can b e more readily solved. The up dated guess of image m at the end of this k -th iteration is denoted as m ( k +1) , and m ( k +1) is then used for the next iteration in Eq. 2. This pro xi- mal problem is a simple soft-thresholding step for specific regularization functions such as R ( m ) = k Ψ m k 1 where Ψ is a W a velet transform [22]. Previously , this regularization function has b een hand-crafted and hand-tuned for every sp ecific application. F or example, spatial wa velets is a p opular choice for general t wo and three dimensional images [4], spatial total v ariation (or finite differences) for angiography [4], and temp oral total v ariation or sparsity in the temp oral F ourier space for cardiac motion [8]. Unfortunately , the design and testing of differen t regularization functions require significan t engineering effort, and its practicality is hampered b y the need to empirically tune the asso ciated regularization parameters. 2.2 Data-Driv en Learning T o b e able to dev elop fast and robust reconstruction algorithms for different MRI sequences and scans, a compelling alternative is to tak e a data-driven approac h to learn the opti- mal regularization functions and parameters. Though it ma y be possible to directly learn this regularization function, a simpler and more straigh t-forw ard approac h is to learn the pro ximal step in Eq. 3 whic h will implicitly learn b oth the regularization function and regularization parameter. In this setup, the pro ximal step is replaced with a deep neural net work to be learned: E θ k where θ k are the learned parameters. Eq. 3 b ecomes m ( k +1) = E θ k ( m ( k +) ) . (4) The steps of Eqs. 2 and 4 can b e unrolled with a fixed num b er of iterations and be denoted as mo del G θ ( y , A ) with inputs of measuremen ts y and imaging mo del A . The training of such a model can then b e p erformed using the follo wing loss function: min θ X i k G θ ( y i , A i ) − m i k 2 2 , (5) 5 where m i is the i -th ground truth example that is retrosp ectiv ely subsampled b y a sampling mask M in the k-space domain to generate y i = M A H i m i . F or deplo yment, new scans are acquired according to the sampling mask M , and the measured data y i can b e used to reconstruct images ˆ m i as ˆ m i = G θ ( y i , A i ) , (6) where A i con tains sensitivity profile maps S i that can b e estimated using algorithms like JSENSE [2] or ESPIRiT [3]. An example netw ork is shown in Fig. 3AB. C m i Compl ex to Channel s Compl ex to Channel s Featu re extract ion for adversa rial loss R esBlock (32 ) Batch Norm ReLU 1x1 Conv, 64 stride 2 R esBlock (64 ) Batch Norm ReLU 1x1 Conv, 128 stride 2 R esBlock (12 8) Batch Norm ReLU 1x1 Conv, 256 stride 2 Batch Norm Tanh Compl ex to Channel s Channe ls to Comple x A i y i m i ( k ) Updat e block Proximal block - A i H + x t ( k ) A i H m i ( k +1) Compl ex to Channel s R esBlock (12 8) R esBlock (12 8) R esBlock (12 8) Batch Norm ReLU 3x3 Conv, 2 N c Channe ls to Comple x + B ResBlock ( f ) Batch Norm ReLU 3x3 Conv, f Batch Norm ReLU 3x3 Conv, f + A Figure 3: Neural net work arc hitecture for MRI reconstruction. One residual blo c k (Res- Blo c k) [23] is illustrated in A with f feature maps. The ResBlo ck is used as a building blo c k for the different netw orks. In B , one iteration of the reconstruction netw ork is de- picted where the i -th dataset is passed through the k -th iteration. Matrix A i represen ts the imaging mo del, m i represen ts the dataset in the image domain, and y i represen ts the dataset in the k-space domain. The final output can b e passed through the net work in C to extract feature maps that can b e compared to the feature maps extracted from the ground truth data using the same netw ork. The tanh activ ation function in C is used to ensure that the v alues in the outputted feature maps are within ± 1. 2.3 Neural Netw ork Design for MRI Reconstruction Imp ortan t considerations must b e made when applying deep neural net works to MRI re- construction. This includes handling of complex data [24], circular con volutions in image domain for Cartesian data [21], incorporating the acquisition model, flexibilit y in acquisition parameters and geometry , data normalization strategies, and plausible data augmen tation approac hes. Since the data measured during an MRI scan is complex, complex data is used for the MRI reconstruction pro cess. As a result, the net works used for MRI reconstruction need to handle complex data types. Two approaches can be used to solv e this issue. The first 6 Z er o P adding Cir cular P adding A B 5x5 K ernel Phase Enc ode Figure 4: Circular con volutional lay er for MRI reconstruction. Current deep learning frame- w orks support zero-padding for conv olutional lay ers (shown in A ) to maintain the original dimensions. Ho wev er, data are measured in the frequency space, and the FFT algorithm is used to con vert this data in to the image domain. Thus, for Cartesian imaging, circular con volutions should be p erformed. This circular conv olution can b e p erformed b y first cir- cularly padding the input data as sho wn in B and cropping the result to the original input dimension. approac h is to con vert the complex data in to tw o channels of data. This conv ersion can b e p erformed as concatenating the m agnitude of the data with the phase of the data in the c hannels dimension, or concatenating the real comp onent with the imaginary comp onen t of the data in the channels dimension. This con v ersion can be reversed without loss of data. The h yp othesis is that the complex data prop erties may not need to b e fully mo deled with sufficiently deep net w orks. This first approac h requires no modifications to current deep learning framew orks that do not ha ve complex num b er supp ort; how ever, the known structure of complex n umbers are not fully exploited. The second approach is to build the neural net work with op erations that supp ort complex data. Sev eral efforts hav e b een made to enable complex op erations in the conv olutional la yers (with complex back-propagation) and the activ ation functions [24]. F or demonstration purp oses, w e construct the netw orks in Fig. 3BC with the first approach b ecause of its simplicit y . Another consideration in building neural net works for MRI reconstruction is handling the con volutional operation at the image edges. Assuming zeros b ey ond the image edges for image-domain con volutions is sufficient for most cases, esp ecially when the imaging volume is surrounded by air and when no high-intensit y signal exists near the edges of the field-of- view (FO V). Ho w ever, when these conditions are not satisfied, this assumption may result in residual aliasing artifacts. These artifacts are usually observed when the imaging ob jec t is larger than the F OV. Data are measured in the F ourier space and transformed using the fast F ourier transform algorithm (FFT) to the image domain; as a result, signals are circularly wrapp ed. Thus, circular conv olutions should be p erformed when applying conv olutions in the image domain. F or this purp ose, w e first circularly pad the images b efore applying the con volutional lay ers (Fig. 4B). The padding width is set such that the “v alid” p ortion of the output has the same image dimensions as the input. Lev eraging a more accurate MRI acquisition mo del in the neural net work will impro ve the reconstruction performance and simplify the netw ork arc hitecture. More sp ecifically , data are measured using m ulti-channel coil receiv ers. The sensitivity profile maps of the differen t coil elemen ts can b e leveraged as a strong prior to help constrain the reconstruction problem. F urther, by exploiting the (soft) SENSE mo del [1, 3], the m ulti-channel complex data ( y ) are reduced to a single-c hannel complex image ( m ) b efore each conv olutional 7 net work blo ck denoted as the proximal blo ck in Fig. 3B. As a further benefit, the learned net work can be trained and applied to datasets with different n umbers of coil channels, b ecause the input to the learned conv olutional netw ork blo ck only requires a single-c hannel complex image. Acquisition parameters ma y c hange for different patients, including spatial resolution, matrix size, and image contrast. Therefore, when applied to clinical scans, the neural net works need to b e able to deal with these differen t scan parameters that are determined b y eac h sp ecific clinical application. One common parameter is the acquisition matrix size. Con volutional lay ers are flexible to different input sizes. F or lay ers that are strict in sizes, the input k-space can b e zero-padded to the same size to mak e the learned netw orks flexible to different dimensions of the acquisition matrix. Appropriate data normalization helps impro ve the training and the final performance of the learned mo del. The input data should ideally b e pre-whitened and normalized by an estimate of the noise statistics among the different coil arra y receivers. In the training loss function of Eq. 5, this data whitening helps balance the training examples based on the signal-to-noise ratio of each training example. These statistics can be measured during a fast calibration scan during which data are measured with no RF excitation. Alternativ ely , if this noise information is not av ailable such as for already collected datasets, the noise can b e estimated from the background signal for a fully sampled acquisition. F or simplicity and with some p ossible loss of p erformance, the raw measurement data can also b e normalized according to the L2 norm of the cen tral k-space signals. In our demonstration, w e normalize the input data b y the L2 norm of the central 5 × 5 region of k-space. F or in termediate lay ers, batc h normalization can b e used to minimize sensitivit y to data scaling. As training examples are difficult to obtain, the num b er of av ailable training examples can b e limited. T o address this concern, data augmen tation can b e applied to train the reconstruction net work. Care must b e tak en when applying data augmentation transforma- tions. F or example, data interpolation are needed for random rotations which may degrade the quality of the input data. Other image domain op erations ma y in tro duce unrealistic errors in the measuremen t k-space domain, such as aliasing in the k-space domain. Flipping and transp osing the dataset will preserv e the original data quality; thus, these op erations can b e included in the training. 2.4 Loss F unction The p erformance of the data-driv en approach is highly dependent on the loss function used. The easiest loss functions to use for training are L1 and L2 losses. The L2 loss is describ ed in Eq. 5 and can b e conv erted to an L1 loss by using the L1 norm instead of the L2 norm. The L1 and L2 losses do not adequately capture the idea of structure or perception. More sophisticated loss functions can b e used suc h as using a netw ork pre-trained on nat- ural images to extract “p erceptual” features [25]. Though general, this feature extraction net work should b e trained for the sp ecific problem domain and the sp ecific task at hand. Generativ e adv ersarial net works (GANs) [26] can b e used to mo del the prop erties of the ground truth images and to exploit that information for improving the reconstruction qual- it y . In Fig. 3C, w e constructed a feature extraction net work D ω ( m ) that is trained to extract the necessary features to compare the reconstruction output with the gr ound truth [27]. The 8 Input Output T ruth Princip al c omponents Figure 5: Example feature maps extracted b y the learned net w ork for adversarial loss. A feature extraction net work w as join tly trained with the reconstruction netw ork. This feature extraction net work was used to pro duce feature maps to compute the adv ersarial loss. Extracted features for the subsampled input (first ro w), reconstruction output (second ro w), and fully-sampled ground truth (last ro w) are display ed. The net work extracted 128 differen t feature maps. Here, the dominant 10 principal components of the resulting feature maps are sho wn where all principal comp onents are rotated to b e aligned to the principal comp onen ts of the truth data. training loss function in Eq. 5 b ecomes: min θ max ω X i k D ω ( G θ ( y i , A i )) − D ω ( m i ) k 2 , (7) where parameters ω are optimized to maximize the difference b etw een the reconstructed image G θ ( y i , A i ) and the ground truth data m i . At the same time, parameters θ are optimized to minimize the difference b et ween the reconstruction and the ground truth data after passing b oth these images through net w ork D ω to extract feature maps. The optimization of Eq. 7 consists of alternating betw een the training of parameters ω with θ constan t and the training of parameters θ with ω constan t. W e refer to this training approac h as training with an “adversarial loss.” This min-max loss function in Eq. 7 can b e unstable and difficult to train. The man y tric ks used to train GANs, such as the training describ ed in Ref. [28], can be leveraged here. W e implement tw o main comp onen ts to help stabilize the training [26]. First, G θ is pre-trained using either an L1 or L2 loss so that the parameters in this net work are prop erly initialized. The netw ork is then fine tuned with the adversarial net work using a reduced training rate. Second, a pixel-wise cost function is added to Eq. 7: min θ max ω X i λ k D ω ( G θ ( y i , A i )) − D ω ( m i ) k 2 2 + k G θ ( y i , A i ) − m i k 2 2 . (8) The images b efore passing through D ω can b e considered as additional feature maps to help stabilize the training pro cess. Hyp erparameter λ is used to w eigh b et ween the t w o comp onen ts in the loss function. 3 Demonstration The data should be collected at the p oin t in the imaging pipeline where the mo del will b e deplo yed. F or the purpose of image reconstruction of subsampled datasets, w e would need to collect the ra w measurement k-space data. This data can b e already filtered with 9 Input L1 L o ss L2 L o ss Adv er sarial P SNR: 36 .6 NRMSE: 0.36 S SIM: 0.82 P SNR: 37 .4 NRMSE: 0.33 S SIM: 0.88 P SNR: 37 .1 NRMSE: 0.34 S SIM: 0.87 P SNR: 37 .4 NRMSE: 0.33 S SIM: 0.88 y T ruth Sampling R = 9.4 kz ky x Figure 6: Demonstration results of a volumetric knee dataset that is subsampled with an acceleration factor R of 9.4 with corner cutting (effective R or 12). In the b ottom ro w, the v olume was reconstructed slice by slice using three net works trained with different loss functions. The reconstruction using the netw ork trained using the L1 and L2 losses yielded the b est results in terms of PSNR, NRMSE, and SSIM. How ever, the reconstruction using the netw ork trained with the adversarial loss yielded results with most realistic texture. the anti-aliasing readout filter and pre-whitened. Typically , this raw imaging data are not readily av ailable as only magnitude images are sa ved as DICOM images in hospital imag- ing database. F urthermore, these stored magnitude images are often processed with image filters, and accurate simulation of raw imaging data from DICOM images is difficult to p erform if not impossible especially in sim ulating realistic phase information. T o facili- tate developmen t, a n umber of differen t op en data initiativ es hav e b een recently launched, including mridata.org [29] and fastMRI [30]. These resources pro vide an initial starting p oin t for developmen t, but more datasets of v arying contrasts, field strengths, and vendors are needed. F or demonstration of the data-driven reconstruction and for enabling repro ducibilty of the results, w e do wnloaded 20 fully sampled volumetric knee datasets [31] that are freely a v ailable in the database of MRI ra w data, mridata.org [29]. Each v olume was collected with 320 slices in the readout direction in x , and eac h of these x -slices w ere treated as separate examples during training and v alidation. The datasets w ere divided by sub ject: 15 sub jects for training (4800 x -slices), 2 sub jects for v alidation (640 x -slices), and 3 sub jects for testing. Sensitivit y maps were estimated using JSENSE [2]. P oisson-disc sampling masks [4] were generated using an acceleration factor R of 9.4 with corner cutting (effectiv e R of 12) and a fully-sampled calibration region of 20 × 20. The netw orks in Fig. 3 were implemen ted in p ython using the T ensorFlo w frame- w ork [32]. Additional reconstruction comp onents w ere performed using the SigPy p ython 10 Clinical practice Data-driven learning Data collection: Measurement of raw Simulation: k-Space subsampling Normalization and augmentation Model design, implementation, and training Pre-clinical evaluation Clinical deployment Clinical evaluation Se t up s c an p ar ame t er s Machine le arning module Simpl e & r apid r ec onst ructi on SCAN OPER A T OR MRI SCANNER C OMPUTER SER VER Init ial f eedb ack windo w Upda t ed f eedb ack windo w ... ... Machine le arning updat e module Compl e t e XD flow r ec onst ructi on ... ... Pr e vious dat a Ti me, t Contr ast in jec ti on S T AR T SC AN S T OP SCAN b a c d e f g h i k z k y Data preparation A B C Figure 7: Example developmen t cycle for the data-driven learning framework. In A , fully sampled raw k-space datasets are first collected and used to p erform realistic simulation of subsampling. The reconstruction algorithm is then built and trained in B . Lastly , the reconstruction mo del is deploy ed and ev aluated in C . pac k age 2 [33]. The reconstruction netw ork was built using 4 iterations; this setup allo wed for relatively faster training for demonstration purp oses. Performance can b e improv ed by increasing the num b er of iterations. T raining and exp erimen ts were p erformed on a single NVIDIA Titan Xp graphics card with 12GB of memory which supp orted a batc h size of 2. The netw ork w as trained multiple times using different loss functions: an L1 loss for 20k steps and L2 loss for 20k steps. F or adv ersarial loss, the net work w as first pre-trained using an L1 loss for 10k steps and then join tly trained with the adv ersarial loss and an L1 loss for 60k steps (10k steps were for the reconstruction net work and 50k steps w ere for the adv ersarial feature extraction netw ork). The data preparation, training, v alidation, and testing python scripts are a v ailable on GitHub. Example results are sho wn in Fig. 6. 4 Clinical In tegration 4.1 Example Setup Man y of these algorithms are computational in tensive, but these algorithms can lev erage off- the-shelf consumer hardware. The bulk of the computational burden is now in the net work, and graphic pro cessing units (GPUs) can b e leveraged for this purp ose. The compute system only needs the minimal central pro cessing unit (CPU) and memory requiremen ts to supp ort the GPUs. Using an NVIDIA Titan Xp graphics card, one 320 × 256 slice to ok on a verage 0.1 seconds to reconstruct, 1.8 seconds for a batc h of 16 slices, and a total of 36.0 seconds for the entire v olume with 320 slices. These benchmarks include the time to transfer the data to and from the CPU to the GPU. The computational sp eed can b e 2 https://github.com/mikgroup/sigpy 11 Compr essed Sensing 3D R ec on Ne tw ork R e f er enc e R e f er enc e (Un- Z oomed) Magnitude (x -y) y-t Figure 8: Clinical deplo yment of a 3D spatiotemp oral conv olutional reconstruction net- w ork for tw o-dimensional cardiac-resolved imaging. Ab ov e, a prosp ectively acquired short- axis dataset with 12.6 fold subsampling w as reconstructed using compressed sensing (L1- ESPIRiT) [3] and a trained 3D reconstruction netw ork. These reconstructions w ere com- pared with the fully sampled reference images (acquired in a separate scan). The 3D net- w ork was able to reconstruct fine structures such as papillary m uscles and trab eculae with less blurring than L1-ESPIRiT (white arrows). F urthermore, the images from the net work depict more natural cardiac motion compared to L1-ESPIRiT with total v ariation, which in tro duced blo ck artifacts (yello w arrows). further improv ed using new er GPU hardware and/or more cards. Alternativ ely , inference on the reconstruction netw ork can b e performed on the CPU to b e able to lev erage more memory in a more cost effectiv e manner, but this setup comes with slow er inference sp eed. 4.2 Clinical Cases Through initial developmen ts, w e hav e deplo y ed a n um b er of different mo dels at our clin- ical site. Example results of data-driven compressed sensing with v ariational net works (VN) [14, 15] are shown in Fig. 2. Con ven tional compressed sensing using L1-ESPIRiT [3] required tuning of the regularization parameter λ . Con ven tional compressed sensing re- construction ac hieved high noise when the regularization parameter was too low (0.0005, left column in Fig. 2), and high residual and blurring artifacts when the regularization pa- rameter was to o high (0.05, right column in Fig. 2). Optimal v alue of this parameter ma y v ary with different scans. Compared to con ven tional compressed sensing, learning a VN (middle column in Fig. 2) achiev ed proper regularization without the need of tuning the regularization parameter. The b est image prior for multi-dimensional imaging is more difficult to engineer, but these larger datasets b enefits tremendously from subsampling since these datasets take longer to acquire. Th us, for m ulti-dimensional space, a comp elling approach is to use data-driv en learning. In Fig. 8, a reconstruction net work consisting of 3D spatiotemp oral con volutions was trained on 12 fully-sampled, breath-held, multi-slice cardiac cine datasets acquired with a balanced steady-state free precession (bSSFP) readout at m ultiple scan orien tations. The reconstruction netw ork enabled the acquisition of a cine slice in a single heartb eat and a full stack of cine slices in a single breathhold [34]. 12 4.3 Ev aluation F or dev elopment and protot yping, image metrics provide a straightforw ard metho d to ev al- uate p erformance. These metrics include p eak signal-to-noise ratio (PSNR), normalized ro ot-mean-square error (NRMSE), and structural similarity (SSIM) [35]. These metrics are describ ed in the App endix. How ever, what is more imp ortan t in the clinical setting is the abilit y for clinicians to mak e an informed decision. Therefore, we recommend performing clinical studies. F or comparing metho ds, the ev aluation of clinical studies can b e p erformed with blinded grading of image quality by multiple radiologists. This ev aluation can b e based on set criteria including ov erall image quality , perceived signal-to-noise ratio, image con trast, image sharpness, residual artifacts, and confidence of diagnosis [15]. Scores from 1 to 5 or -2 to 2 are given for each clinical patient scan with a sp ecific definition for eac h score for ob jectivit y . An example of scoring for signal-to-noise ratio and image sharpness is describ ed in T able 1. A total num b er of more than 20 consecutiv e scans are usually collected to ac hieve a comprehensiv e and statistically meaningful ev aluation for an initial study . The difference b et ween t wo reconstruction approac hes can b e statistically tested with Wilco xon tests on the null hypothesis that there w as no significan t difference b et ween tw o approac hes. Score Ov erall Image Quality Signal-to-Noise Ratio Sharpness 1 Nondiagnostic All structures app ear to b e to o noisy . Some structures are not sharp on most images. 2 Limited Most structures app ear to b e to o noisy . Most structures are sharp on some images. 3 Diagnostic F ew structures app ear to b e to o noisy on most images. Most structures are sharp on most images. 4 Go od F ew structures app ear to b e to o noisy on a few images. All structures are sharp on most images. 5 Excellen t There is no noticeable noise on any of the images. All structures are sharp on all images. T able 1: Example scoring criteria to ev aluate metho ds. 5 Discussion This data-driv en approach to accelerated imaging has the p oten tial to eliminate many of the c hallenges asso ciated with compressed sensing, including the need to design hand- crafted priors and to hand tune the regularization functions. In addition, the describ ed framew ork inspired by iterativ e inference algorithms provides a principled approach for designing reconstruction netw orks. F urthermore, this approach is useful to in terpret the trained net work components. Ho wev er, the data-driv en reconstruction framew ork faces new challenges, including data scarcity , generalizability , reliability , and meaningful metrics. High qualit y training lab els, or fully sampled datasets, are scarce esp ecially in the clin- ical settings where patient motion impacts image quality and where length y fully-sampled acquisitions are impractical to p erform if a faster solution exists. T o address data scarcit y , an imp ortant future direction p ertains to designing compact net work architectures that are effectiv e with small training labels and are p ossibly trained in an unsup ervised fashion. Early attempts are made in [36] where recurrent neural netw orks are lev eraged for learn- 13 ing proximal operators using only a couple of residual blo cks that p erforms well for small training sample sizes. The net work can be trained or re-trained for differen t anatomies and different t yp es of scan con trasts. This strategy can b e implemen ted if sufficien t training datasets for all differen t settings are readily a v ailable. Giv en the scarcity of training examples and the cost to collect these examples, another strategy is to design and train netw orks that are highly generalizable. This generalizabilit y can b e ac hieved with smaller net works as discussed b efore or by training with a larger image manifold such as using natural images. With this approac h, a loss in p erformance may b e observed since a larger than necessary manifold is learned. In the worst case, the images reconstructed using this data-driv en approac h should not b e worse than images reconstructed using a more conserv ative approach such as zero-filling or parallel imaging [1, 3]. Regarding generalizability , an imp ortan t question concerns reconstruction p erformance for sub jects with unseen abnormalities. Patien t abnormalities can b e quite heterogeneous, and these abnormalities are rare and unlik ely to b e included in the training dataset. If not designed carefully , generative netw orks hav e the p ossibility of removing or creating critical features that will result in misdiagnosis. The optimization-based net work architecture uti- lizes data consistency , as exemplified b y the gradient up date step, for the image reco very problem. How ever, the inheren t ambiguit y of ill-p osed problems do es not guarantee faith- ful reco v ery . Therefore, a systematic study is required to analyze the reco v ery of images. Also, effective regularization techniques (p ossibly through adversarial training) are needed to av oid missing important diagnostic information. More efforts are needed to dev elop a holistic quality score capturing the uncertaint y in the acquisition sc heme and training data. Dev eloping a standard un biased metric for medical images that assesses the authen ticity of medical images is extremely imp ortant. Here, w e discuss the use of differen t loss functions to train the net work and common imaging metrics to ev aluate the images. Additionally , we discuss an example of a p ossible clinical ev aluation that can b e performed to assess the algo- rithm in the clinical setting. Ho wev er, the reconstruction task should ideally be optimized and ev aluated for the end task whic h can b e consisted of detection and quan tification. With significan t effort in automating image in terpretation, this data-driven framew ork pro vides an opp ortunity to pursue the ability to train the reconstruction end-to-end for the ultimate goal of improving patient care. In conclusion, deep learning has the potential to increase the accessibility and general- izabilit y of fast imaging through data subsampling. Previous c hallenges with compressed sensing can b e approached with a data-driven framework to create a solution that is more readily translated to clinical practice. App endix: Imaging Metrics Imaging metrics are commonly used to ev aluate results. This includes p eak signal-to-noise ratio (PSNR) and normalized ro ot-mean-square error (NRMSE). F or these quantities, w e 14 use the following equations: MSE( x, x r ) = 1 N M N X i M X j | x [ i, j ] − x r [ i, j ] | 2 , (9) PSNR( x, x r ) = 10 log 10  max( | x r | 2 ) / MSE  , (10) NRMSE( x, x r ) = √ MSE / v u u t 1 N M N X i M X j | x r [ i, j ] | 2 , (11) where x denotes the test image, x [ i, j ] denotes the v alue of the test image at pixel ( i, j ), and x r denotes the reference ground truth image. Ac kno wledgmen t The authors w ould lik e to thank GE Healthcare, NIH R01-EB009690, NIH R01-EB026136, and NIH R01-EB019241 for the research supp ort. References [1] K. P . Pruessmann, M. W eiger, M. B. Scheidegger, and P . Bo esiger, “SENSE: sensitivit y enco ding for fast MRI.” Magnetic R esonanc e in Me dicine , v ol. 42, no. 5, pp. 952–62, 11 1999. [Online]. Av ailable: h ttp://www.ncbi.nlm.nih.gov/pubmed/10542355 [2] L. Ying and J. Sheng, “Join t image reconstruction and sensitivit y estimation in SENSE (JSENSE),” Magnetic R esonanc e in Me dicine , vol. 57, no. 6, pp. 1196–1202, 6 2007. [Online]. Av ailable: http://doi.wiley .com/10.1002/mrm.21245 [3] M. Uec ker, P . Lai, M. J. Murphy , P . Virtue, M. Elad, J. M. Pauly , S. S. V asana wala, and M. Lustig, “ESPIRiT-an eigen v alue approac h to auto calibrating parallel MRI: Where SENSE meets GRAPP A,” Magnetic R esonanc e in Me dicine , vol. 71, no. 3, pp. 990–1001, 3 2014. [Online]. Av ailable: http://doi.wiley .com/10.1002/mrm.24751 [4] M. Lustig, D. Donoho, and J. M. Pauly , “Sparse MRI: The application of compressed sensing for rapid MR imaging.” Magnetic R esonanc e in Me dicine , vol. 58, no. 6, pp. 1182–1195, 12 2007. [Online]. Av ailable: h ttp://www.ncbi.nlm.nih.gov/pubmed/ 17969013 [5] J. Y. Cheng, K. Hanneman, T. Zhang, M. T. Alley , P . Lai, J. I. T amir, M. Uec ker, J. M. Pauly , M. Lustig, and S. S. V asanaw ala, “Comprehensiv e motion-comp ensated highly accelerated 4D flo w MRI with ferumo xytol enhancemen t for p ediatric congenital heart disease.” Journal of Magnetic R esonanc e Imaging , vol. 43, no. 6, pp. 1355–1368, 6 2016. [Online]. Av ailable: http://www.ncbi.nlm.nih.go v/pubmed/26646061 [6] S. S. V asanaw ala, K. Hanneman, M. T. Alley , and A. Hsiao, “Congenital heart disease assessmen t with 4D flo w MRI.” Magnetic R esonanc e in Me dicine , vol. 42, no. 4, pp. 870–86, 10 2015. [Online]. Av ailable: h ttp://www.ncbi.nlm.nih.gov/pubmed/25708923 [7] H. Jung, K. Sung, K. S. Nay ak, E. Y. Kim, and J. C. Y e, “k-t fo cuss: a general compressed sensing framework for high resolution dynamic mri,” Magnetic r esonanc e in me dicine , vol. 61, no. 1, pp. 103–116, 2009. 15 [8] L. F eng, M. B. Srichai, R. P . Lim, A. Harrison, W. King, G. Adluru, E. V. R. Dib ella, D. K. So dickson, R. Otazo, and D. Kim, “Highly accelerated real-time cardiac cine MRI using k-t SP ARSE-SENSE.” Magnetic R esonanc e in Me dicine , v ol. 70, no. 1, pp. 64–74, 7 2013. [Online]. Av ailable: http://www.ncbi.nlm.nih.go v/pubmed/22887290 [9] S. G. Lingala, Y. Hu, E. DiBella, and M. Jacob, “Accelerated dynamic MRI exploiting sparsit y and low-rank structure: k-t SLR.” IEEE T r ansactions on Me dic al Imaging , v ol. 30, no. 5, pp. 1042–1054, 5 2011. [Online]. Av ailable: h ttp://www.ncbi.nlm.nih.gov/pubmed/21292593 [10] T. Zhang, J. Y. Cheng, A. G. P otnick, R. A. Barth, M. T. Alley , M. Uec ker, M. Lustig, J. M. P auly , and S. S. V asanaw ala, “F ast p ediatric 3D free-breathing ab dominal dynamic con trast enhanced MRI with high spatiotemp oral resolution,” Journal of Magnetic R esonanc e Imaging , v ol. 41, no. 2, pp. 460–473, 2 2015. [Online]. Av ailable: http://www.ncbi.nlm.nih.go v/pubmed/24375859 [11] L. F eng, L. Axel, H. Chandarana, K. T. Blo ck, D. K. Sodickson, and R. Otazo, “XD-GRASP: Golden-angle radial MRI with reconstruction of extra motion-state dimensions using compressed sensing.” Magnetic R esonanc e in Me dicine , vol. 75, no. 2, pp. 775–88, 2 2016. [Online]. Av ailable: h ttp: //www.ncbi.nlm.nih.go v/pubmed/25809847 [12] J. Y. Cheng, T. Zhang, M. T. Alley , M. Uec ker, M. Lustig, J. M. P auly , and S. S. V asanaw ala, “Comprehensiv e Multi-Dimensional MRI for the Sim ultaneous Assessmen t of Cardiopulmonary Anatom y and Ph ysiology,” Scientific R ep orts , vol. 7, no. 1, p. 5330, 2017. [Online]. Av ailable: http://dx.doi.org/10.1038/s41598- 017- 04676- 8 [13] S. S. V asanaw ala, M. T. Alley , B. A. Hargreav es, R. A. Barth, J. M. P auly , and M. Lustig, “Impro ved Pediatric MR Imaging with Compressed Sensing,” R adiolo gy , vol. 256, no. 2, pp. 607–616, 8 2010. [Online]. Av ailable: h ttp://www.ncbi.nlm.nih.gov/pubmed/20529991 [14] K. Hammernik, T. Klatzer, E. Kobler, M. P . Rec ht, D. K. So dic kson, T. Pock, and F. Knoll, “Learning a v ariational netw ork for reconstruction of accelerated MRI data,” Magnetic R esonanc e in Me dicine , v ol. 79, no. 6, pp. 3055–3071, 6 2018. [15] F. Chen, V. T a viani, I. Malkiel, J. Y. Cheng, J. I. T amir, J. Shaikh, S. T. Chang, C. J. Hardy , J. M. Pauly , and S. S. V asanaw ala, “V ariable-Density Single-Shot F ast Spin-Ec ho MRI with Deep Learning Reconstruction by Using V ariational Net works,” R adiolo gy , v ol. 289, no. 2, pp. 366–373, 11 2018. [Online]. Av ailable: h ttp://pubs.rsna.org/doi/10.1148/radiol.2018180445 [16] S. Ramani, Z. Liu, J. Rosen, J.-F. Nielsen, and J. A. F essler, “Regularization parameter selection for nonlinear iterative image restoration and mri reconstruction using gcv and sure-based metho ds,” IEEE T r ansactions on Image Pr o c essing , v ol. 21, no. 8, pp. 3659– 3672, 2012. [17] A. Beck and M. T eb oulle, “A fast iterative shrink age-thresholding algorithm for linear in verse problems,” SIAM J Imaging Sci , v ol. 2, no. 1, pp. 183–202, 2009. [Online]. Av ailable: http://epubs.siam.org/doi/pdf/10.1137/080716542 16 [18] Y. Y ang, J. Sun, H. Li, and Z. Xu, “ADMM-Net: A Deep Learning Approach for Compressive Sensing MRI,” NIPS , pp. 10–18, 5 2017. [Online]. Av ailable: h [19] J. Adler and O. Oktem, “Learned Primal-Dual Reconstruction,” IEEE T r ansactions on Me dic al Imaging , vol. 37, no. 6, pp. 1322–1332, 6 2018. [Online]. Av ailable: h ttps://ieeexplore.ieee.org/do cument/8271999/ [20] H. K. Aggarwal, M. P . Mani, and M. Jacob, “Mo dl: Mo del based deep learning arc hi- tecture for inv erse problems,” IEEE tr ansactions on me dic al imaging , 2018. [21] J. Y. Cheng, F. Chen, M. T. Alley , J. M. Pauly , and S. S. V asanaw ala, “Highly Scalable Image Reconstruction using Deep Neural Net works with Bandpass Filtering,” arXiv:1805.03300 [cs.CV] , 5 2018. [Online]. Av ailable: h [22] I. Daubechies, M. Defrise, and C. De Mol, “An iterativ e thresholding algorithm for linear inv erse problems with a sparsity constraint,” Communic ations on Pur e and Applie d Mathematics , v ol. 57, no. 11, pp. 1413–1457, 11 2004. [Online]. Av ailable: h ttp://doi.wiley .com/10.1002/cpa.20042 [23] K. He, X. Zhang, S. Ren, and J. Sun, “Iden tity Mappings in Deep Residual Netw orks,” arXiv:1603.05027 [cs.CV] , 3 2016. [Online]. Av ailable: h [24] C. T rab elsi, O. Bilaniuk, Y. Zhang, D. Serdyuk, S. Subramanian, J. F. San tos, S. Mehri, N. Rostamzadeh, Y. Bengio, and C. J. Pal, “Deep Complex Net works,” arXiv:1705.09792 [cs.NE] , pp. 1–19, 5 2017. [Online]. Av ailable: h [25] J. Johnson, A. Alahi, and L. F ei-F ei, “P erceptual Losses for Real-Time St yle T ransfer and Sup er-Resolution,” in ECCV , 3 2016. [Online]. Av ailable: h [26] M. Mardani, E. Gong, J. Y. Cheng, S. S. V asanaw ala, G. Zaharc huk, L. Xing, and J. M. Pauly , “Deep Generative Adversarial Neural Net works for Compressiv e Sensing MRI,” IEEE T r ansactions on Me dic al Imaging , vol. 38, no. 1, pp. 167–179, 1 2019. [Online]. Av ailable: h ttps://ieeexplore.ieee.org/do cument/8417964/ [27] K. Hammernik, E. Kobler, T. Pock, M. P . Rech t, D. K. So dic kson, and F. Knoll, “V ariational Adv ersarial Netw orks for Accelerated MR Image Reconstruction,” in Joint A nnual Me eting ISMRM-ESMRMB , Paris, F rance, 2018, p. 1091. [28] I. Gulra jani, F. Ahmed, M. Arjovsky , V. Dumoulin, and A. Courville, “Improv ed T raining of W asserstein GANs,” , 2017. [Online]. Av ailable: h [29] F. Ong, S. Amin, S. S. V asana wala, and M. Lustig, “An Open Archiv e for Sharing MRI Ra w Data,” in ISMRM & ESMRMB Joint Annual Me eting , Paris, F rance, 2018, p. 3425. 17 [30] J. Zb ontar, F. Knoll, A. Sriram, M. J. Muckley , M. Bruno, A. Defazio, M. Paren te, K. J. Geras, J. Katsnelson, H. Chandarana, Z. Zhang, M. Drozdzal, A. Romero, M. Rabbat, P . Vincen t, J. Pink erton, D. W ang, N. Y akub ov a, E. Owens, C. L. Zitnic k, M. P . Rec ht, D. K. So dic kson, and Y. W. Lui, “fastMRI : An Op en Dataset and Benchmarks for Accelerated MRI,” arXiv: 11811.08839 [cs.CV] , pp. 1–29, 2018. [31] K. Epp erson, A. M. Sawy er, M. Lustig, M. T. Alley , M. Ueck er, P . Virtue, P . Lai, and S. S. V asana wala, “Creation of F ully Sampled MR Data Rep ository for Compressed Sensing of the Knee,” in SMR T 22nd A nnual Me eting , Salt Lak e Cit y , Utah, USA, 2013. [32] M. Abadi, A. Agarwal, P . Barham, E. Brevdo, Z. Chen, C. Citro, G. S. Corrado, A. Davis, J. Dean, M. Devin, S. Ghema wat, I. Goo dfellow, A. Harp, G. Irving, M. Isard, Y. Jia, R. Jozefo wicz, L. Kaiser, M. Kudlur, J. Lev enberg, D. Mane, R. Monga, S. Mo ore, D. Murra y , C. Olah, M. Sc huster, J. Shlens, B. Steiner, I. Sutsk ever, K. T alw ar, P . T uc ker, V. V anhouck e, V. V asudev an, F. Viegas, O. Vin y als, P . W arden, M. W attenberg, M. Wic ke, Y. Y u, and X. Zheng, “T ensorFlo w: Large-Scale Mac hine Learning on Heterogeneous Distributed Systems,” arXiv:1603.04467 [cs.DC] , 3 2016. [Online]. Av ailable: [33] F. Ong and M. Lustig, “SigPy: A Python Pac k age for High P erformance Iterative Reconstruction,” in ISMRM Annual Me eting & Exhibition , Montreal, Canada, 2019. [34] C. Sandino, P . Lai, M. A. Janich, A. C. S. Brau, S. S. V asanaw ala, and J. Y. Cheng, “ESPIRiT with deep priors: Accelerating 2D cardiac CINE MRI b ey ond compressed sensing,” in SCMR/ISMRM Workshop , Seattle, W ashington, USA, 2019. [35] Z. W ang, A. C. Bo vik, H. R. Sheikh, and E. P . Simoncelli, “Image qualit y assessment: from error visibilit y to structural similarit y .” IEEE T r ans Image Pr o c essing , vol. 13, no. 4, pp. 600–612, 4 2004. [Online]. Av ailable: h ttp://www.ncbi.nlm.nih.gov/pubmed/15376593 [36] M. Mardani, Q. Sun, S. V asa wanala, V. Pap yan, H. Mona jemi, J. Pauly , and D. Donoho, “Neural Proximal Gradient Descen t for Compressiv e Imaging,” Confer enc e on Neur al Information Pr o c essing Systems , v ol. 8166, p. 111, 2018. [Online]. Av ailable: h 18