Deep-learning-based Optimization of the Under-sampling Pattern in MRI
Cagla D. Bahadir Alan Q. Wang Adrian V. Dalca Mert R. Sabuncu
eess.IV cs.LG stat.ML 68T01






Deep-learning-based Optimization of the 
 Under-sampling Pattern in MRI
    Cagla D. Bahadir*, Alan Q. Wang*, Adrian V. Dalca, and Mert R. Sabuncu Cagla D. Bahadir and Mert R. Sabuncu are with the Meinig School of Biomedical Engineering, Cornell University. Alan Q. Wang and Mert R. Sabuncu are with the School of Electrical and Computer Engineering, Cornell University. Adrian V. Dalca is with the Computer Science and Artificial Intelligence Lab at the Massachusetts Institute of Technology and the A.A. Martinos Center for Biomedical Imaging at the Massachusetts General Hospital.
* indicates equal contribution. Cagla D. Bahadir and Alan Q. Wang are co-first authors.



    December 30, 2023
==========================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================================




In compressed sensing MRI (CS-MRI), k-space measurements are under-sampled to achieve accelerated scan times. 
CS-MRI presents two fundamental problems: (1) where to sample and (2) how to reconstruct an under-sampled scan.
In this paper, we tackle both problems simultaneously for the specific case of 2D Cartesian sampling, using a novel end-to-end learning framework that we call LOUPE (Learning-based Optimization of the Under-sampling PattErn).
Our method trains a neural network model on a set of full-resolution MRI scans, which are retrospectively under-sampled on a 2D Cartesian grid and forwarded to an anti-aliasing (a.k.a. reconstruction) model that computes a reconstruction, which is in turn compared with the input.
This formulation enables a data-driven optimized under-sampling pattern at a given sparsity level.
In our experiments, we demonstrate that LOUPE-optimized under-sampling masks are data-dependent, varying significantly with the imaged anatomy, and perform well with different reconstruction methods.
We present empirical results obtained with a large-scale, publicly available knee MRI dataset, where LOUPE offered superior reconstruction quality across different conditions.
Even with an aggressive 8-fold acceleration rate, LOUPE's reconstructions contained much of the anatomical detail that was missed by alternative masks and reconstruction methods.
Our experiments also show how LOUPE yielded optimal under-sampling patterns that were significantly different for brain vs knee MRI scans.
Our code is made freely available at <https://github.com/cagladbahadir/LOUPE/>.



Compressed Sensing, Magnetic Resonance Imaging, Deep Learning




§ INTRODUCTION


Magnetic Resonance Imaging (MRI) is a ubiquitous, non-invasive, and versatile biomedical imaging technology. 
A central challenge in MRI is long scan times, which constrains accessibility and leads to high costs. 

One remedy is to accelerate MRI via compressed sensing <cit.>.
In compressed sensing MRI, k-space data (i.e., the Fourier transform of the image) is sampled below the Nyquist-Shannon rate <cit.>, which is often referred to as “under-sampling.”
Given an under-sampled set of measurements, the objective is to “reconstruct” the full-resolution MRI.
This is the main problem that most of the compressed sensing literature is focused on and is conventionally formulated as an optimization problem that trades off two objectives: one that quantifies the fit between the measurements and the reconstruction (sometimes referred to as data consistency), and another that captures prior knowledge on the distribution of MRI data. 
This latter objective is often achieved via the incorporation of regularization terms, such as the total variation penalty and/or a sparsity-inducing norm on transformation coefficients, like wavelets or a dictionary decomposition <cit.>. 
Such regularization functions aim to capture different properties of real-world MR images, which in turn help the ill-posed  reconstruction problem by guiding to more realistic solutions. 
One approach to develop a data-driven regularization function is to construct a sparsifying dictionary, for example, based on image patches <cit.>.


Once the optimization problem is set up, the reconstruction algorithm often iteratively minimizes the regularization term(s) and enforces data consistency in k-space. This is classically solved for each acquired dataset, independently, and from scratch - a process that can be computationally demanding. 








As we describe in the following section, there has been a recent surge in machine learning based methods that take a different, and often computationally more efficient approach to solving the reconstruction problem.
These techniques are gradually becoming more widespread and expected to complement existing regularized optimization based approaches.




Another critical component of compressed sensing MRI is the under-sampling pattern. 
For a given acceleration rate, there is an exponentially large number of possible patterns one can implement for under-sampling. 
Each of these under-sampling patterns will in general lead to different reconstruction performance that will depend on the statistics of the data and the utilized reconstruction method. 


One way to view this is to regard the reconstruction model as imputing the parts of the Fourier spectrum that was not sampled.
For example, a frequency component that is constant for all possible datasets we might observe, does not need to be sampled when we use a reconstruction model that can leverage this information.
On the other hand, highly variable parts of the spectrum will likely need to be measured to achieve accurate reconstructions.
This simple viewpoint ignores the potentially multi-variate nature of the distribution of k-space measurements.
For instance, missing parts of the Fourier spectrum might be reliably imputed from other acquired parts, due to strong statistical dependency.
Widely used under-sampling strategies in compressed sensing MRI include Random Uniform <cit.>, Variable Density <cit.> and equi-spaced Cartesian <cit.> with skipped lines. 
These techniques, however, are often implemented heuristically, not in a data-driven adaptive fashion, and their popularity is largely due to their ease of implementation and good performance when coupled with popular reconstruction methods.




As we discuss below, some recent efforts compute optimized under-sampling patterns - a literature that is closely related to our primary objective.

In this paper, we propose to merge the two core problems of compressed sensing MRI: (1) optimizing the under-sampling pattern and (2) reconstruction.
We consider these two problems simultaneously because they are closely inter-related.
Specifically, the optimal under-sampling pattern should, in general, depend on the reconstruction method and vice versa.
Inspired by recent developments in machine learning, we propose a novel end-to-end deep learning strategy to solve the combined problem.
We call our method LOUPE, which stands for Learning-based Optimization of the Under-sampling PattErn.
A preliminary version of LOUPE was published as a conference paper <cit.>. 
In this journal paper, we present an extended treatment of the literature, an important modification to the LOUPE objective that enables us to directly set the desired sparsity level while obviating the need to identify the value of an extra hyper-parameter, more details on the methods, and new experimental results.


The rest of the paper is organized as follows.
In the following two sections, we review two closely related bodies of work: a rapidly growing list of recent papers that use machine learning for efficient and accurate reconstruction; and several proposed approaches to optimize the under-sampling pattern.
Section <ref> then presents the mathematical and implementation details of the proposed method, LOUPE. 
Section <ref> presents our experiments, where we compare the performance obtained with several benchmark reconstruction methods and under-sampling patterns. 
Section <ref> concludes with a discussion.




§ MACHINE LEARNING FOR UNDER-SAMPLED IMAGE RECONSTRUCTION



Over the last few years, machine learning methods have been increasingly used for medical image reconstruction <cit.>, including compressed sensing MRI. 
An earlier example is the Bayesian non-parametric dictionary learning approach <cit.>. 
In this framework, dictionary learning is used in the reconstruction problem to obtain a regularization objective that is customized to the image at hand. The reconstruction problem is solved via a traditional optimization approach.

More recently, fueled by the success of deep learning, several machine learning techniques have been proposed to implement efficient and accurate reconstruction models.
For instance, in a high profile paper, Zhu et al. demonstrated a versatile supervised learning approach, called AUTOMAP <cit.>. This method uses a neural network model to map sensor measurements directly to reconstruction images, via a manifold representation. 

The experimental results show that AUTOMAP can produce high quality reconstructions that are robust to noise and other artifacts.
Another machine learning based reconstruction approach involves using a neural network to efficiently execute the iterations of the Alternating Direction Method of Multipliers (ADMM) method that solves a conventional optimization problem <cit.>. This technique, called ADDM-Net, uses an “unrolled” neural network architecture to implement the iterative optimization procedure.
In a similar approach, an unrolled neural network is used to implement the iterations of the Landweber method <cit.>.





Another class of methods rely on the U-Net architecture <cit.> or its variants.
For example U-Net-like models have been trained in a supervised fashion to remove aliasing artifacts in the imaging domain <cit.>.
The U-Net architecture has also been used by appending with a forward physics based model to learn phase masks for depth estimation in other computer vision applications  <cit.>. 
In these methods, the neural network takes as input a poor quality reconstruction (e.g., obtained via a simple inverse Fourier transform applied to zero-filled k-space data) to compute a high quality output, where the objective is to minimize the loss function  (e.g., squared difference) between the output and the ground truth provided during training.

Inspired by the success of Generative Adversarial Networks (GANs) <cit.>, several groups have proposed the use of an adversarial loss, in addition to the more conventional loss functions, to obtain better quality reconstructions <cit.>. 




The method we propose in this paper builds on the neural net-based reconstruction framework, which offers a computationally efficient and differentiable anti-aliasing model.
We combine the anti-aliasing model with a retrospective under-sampling module, and learn both blocks simultaneously. 



§ DATA-DRIVEN UNDER-SAMPLING IN COMPRESSED SENSING MRI



The under-sampling pattern in k-space is closely related to reconstruction performance. 
Several under-sampling patterns are widely used, including Random Uniform <cit.>, Variable Density <cit.> and equi-spaced Cartesian <cit.> with skipped lines. 
Since the early work by Lustig et al. <cit.>, random or stochastic under-sampling strategies are commonly used, as they induce noise-like artifacts that are often easier to remove during reconstruction. 
However, as several papers have pointed out, the quality of reconstruction can also be improved by adapting the under-sampling pattern to the data and application.
We underscore that there is a related, but more general problem of optimizing projections in compressed sensing, which has received considerable attention <cit.>. 
Below, we review adaptive under-sampling strategies in compressed sensing MRI, which is the focus of our paper. 

Knoll et al. proposed an under-sampling strategy that follows the power spectrum of a provided example (template) image <cit.>. 
This method collects denser samples in parts of k-space that has most of the power concentrated, such as the low-frequency components of real-world MRIs. 
The authors argue that most of the anatomical variability is reflected in the phase and not magnitude of the Fourier spectrum, justifying the decision to rely on the magnitude to construct the under-sampling pattern.
In experiments, they apply their method to both knee and brain MRI scans, demonstrating better performance than parametric variable density under-sampling patterns <cit.>.




Kumar Anand et al. presented a second-order cone optimization technique for finding the optimal under-sampling trajectory in volumetric MRI scans <cit.>. 
Their heuristic strategy accounts for coverage, hardware limitations and signal generation in a single convex optimization problem. 
Building on this approach, Curtis et al. employed a genetic algorithm to optimize sampling trajectories in k-space, while accounting for multi-coil configurations <cit.>.

Seeger et al. employed a Bayesian inference approach to optimize Cartesian and spiral trajectories <cit.>. 
Their iterative, greedy technique seeks the under-sampling pattern in k-space that minimizes the uncertainty in the computed reconstruction, conditioned on the acquired measurements. 
In a recent paper that explores a similar direction, Haldar et al. proposed a new framework called Oracle-based Experiment Design for Imaging Parsimoniously Under Sparsity constraints (OEDIPUS) <cit.>.
OEDIPUS solves an integer programming problem to minimize a lower (Cramer-Rao) bound on the variance of the reconstruction. 







Roman et al. offered a novel theoretical perspective on compressed sensing, which underscores the importance of adapting the under-sampling strategy to the structure in the data <cit.>.
The authors also presented an innovative multilevel sampling approach that depends on the resolution and the structure of the data that empirically outperforms competing under-sampling schemes. 




An alternative class of data-driven approaches aims to find the optimal under-sampling pattern that yields the best reconstruction or imputation quality, using some full-resolution training data that is retrospectively under-sampled <cit.>.
The central challenge in this framework is to efficiently solve two computationally challenging nested optimization problems.
The first, outer problem is to identify the optimal under-sampling mask or sampling trajectory that achieves best reconstruction quality, which is in turn the result of an inner optimization step.

Several authors have proposed to solve this nested pair of problems via heuristic, greedy algorithms <cit.>.
In an empirical study, Zijlstra et al. demonstrated that a data-driven optimization approach <cit.> can yield better reconstructions than those obtained with more conventional methods <cit.>.
 However, prior data-driven methods mostly lack the computational efficiency to handle large-scale full-resolution (training) data. 
Furthermore, they are often not flexible enough to deal with different types of under-sampling schemes. 

In this paper, we present a novel,  flexible, and computationally efficient data-driven approach to optimize the under-sampling pattern.
Instead of formulating the problem as two nested optimization problems, we leverage modern machine learning techniques and take an end-to-end learning perspective.
Similar to recent deep learning based reconstruction techniques <cit.>, our implementation employs the U-Net <cit.> architecture, which we append with a probabilistic under-sampling step that is also learned. 
We assume we are provided with full-resolution MRI data, which we retrospectively under-sample. 
We note that there have been contemporaneous efforts <cit.> that build on our prior work <cit.> to take a deep learning-based approach similar to ours for optimizing the under-sampling pattern in compressed sensing.












§ METHOD





 §.§ Learning-based Optimization of Under-sampling Pattern


LOUPE solves the two problems of compressed sensing simultaneously: (1) identifying the optimal under-sampling pattern; and (2) reconstructing from the under-sampled measurements. 

Building on stochastic strategies of compressed sensing, we consider a probabilistic mask P, which describes an independent Bernoulli random variable at each k-space point. 
The probabilistic mask
P is defined on the full-resolution k-space grid, and at each point takes on non-negative continuous probability values, i.e., P∈ [0,1]^d, where d is the total number of grid points.
For example, for a 100 × 100 image, d = 10,000.
We parameterize  P with an unconstrained image O∈ℝ^d, such that P = σ_t(O).
Here, σ_t denotes an element-wise sigmoid, with slope  t, which is treated as a hyper-parameter.
P_i = 1/1+e^-t O_i, where the sub-script i indexes the grid points.
Binary realizations drawn from P represent an under-sampling mask M∈{0,1}^d: i.e., M∼∏_i=1^d ℬ(P_i), where ℬ(p) denotes a Bernoulli random variable with parameter p.
The binary mask M has value 1 for grid points that are acquired and 0 for points that were not acquired. 

Suppose we are provided with a collection of complex-valued full-resolution images, denoted as {x_j ∈ℂ^d}_j=1^n,  where j corresponds to the scan number and n is the total number of scans.
In our following treatment, without loss of generality, we will assume the provided data are in image domain. 
This can be modified to accept raw k-space measurements, by simply removing a forward Fourier transform.
In the LOUPE framework, we solve the following problem:




    min_O, θ𝔼_M∼∏_i=1ℬ(σ_t(O_i))∑_j=1^n  A_θ (F^Hdiag(M) F x_j) - x_j _2^2, 
    such that 1/dσ_t(O)_1 = α


where F ∈ℂ^d × d stands for the (forward) Fourier transform matrix, F^H is the inverse Fourier transform matrix,  A_θ denotes an anti-aliasing reconstruction function parameterized with θ, ·_2 denotes the L2 norm, ·_1 denotes the L1 norm, α∈ [0,1] is the desired sparsity level, and diag(·) denotes a diagonal matrix obtained by setting the diagonal to the argument vector.
The loss function of Equation (<ref>) quantifies the average quality of the reconstruction.
The constraint 1/dσ_t(O)_1 = 1/dP_1 = α ensures that the probabilistic under-sampling mask has an average value of α, which will correspond to an acceleration factor of R= 1/α.
We emphasize that the problem formulation of Equation (<ref>) is slightly different than the formulation in our conference paper <cit.>, where we employed a sparsity inducing regularization instead of a hard constraint.
That formulation required fine-tuning a hyper-parameter in order to obtain the desired sparsity level, which was computationally inefficient.
In this paper we implement a practical strategy, described below, to solve the problem with the hard sparsity constraint directly.


The loss function involves an expectation over the random binary mask M.
We approximate this expectation using a Monte Carlo-based sample averaging strategy:



    min_O, θ∑_j=1^n 1/K∑_k=1^K
     A_θ (F^Hdiag(m^(k)) F x_j) - x_j _2^2,
    such that 1/dσ_t(O)_1 = α


where m^(k) are independent realizations drawn from ∏_iℬ(σ_t(O_i)).

Similar to the re-parameterization trick used in variational techniques, such as the VAE <cit.>, we can re-write Equation (<ref>):


    min_O, θ∑_j=1^n 1/K∑_k=1^K  A_θ (F^Hdiag(U^(k)≤σ_t(O)) F x_j) - x_j _2^2, 
    
    such that 1/dσ_t(O)_1 = α

where U^(k) are independent realizations drawn from ∏_i=1^d 𝒰(0,1), a spatially independent set of uniform random variables on [0,1]. The result of the inequality operation is 1 if the condition is satisfied and 0 otherwise.
In Equation (<ref>) the random draws are thus from a constant distribution (independent uniform) and the probabilistic mask P (through its parameterization O), only effects the thresholding operation. 





 §.§ Implementation



In our first version of LOUPE, we used a 2D U-Net architecture <cit.> for the anti-aliasing function, as depicted in Figure <ref> and similar to other prior work <cit.>. 
To enable a differentiable loss function, we relaxed the non-differentiable threshold operation of Equation (<ref>) with a sigmoid. This relaxation is similar to the one used in recent Gumbel-softmax <cit.> and concrete distributions <cit.> in training neural networks with discrete representations.
Our final objective is:


    min_O, θ∑_j=1^n 1/K∑_k=1^K
     A_θ (F^Hdiag(σ_s(σ_t(O)-U^(k))) F x_j) - x_j _2^2,
    
    such that 1/dσ_t(O)_1 = α


where σ_s(a) = 1/1+e^-sa is a sigmoid with slope s.


A critical issue for LOUPE is the enforcement of the hard constraint 1/dσ_t(O)_1=1/dP_1 = α.


To achieve this we introduce a normalization layer, which rescales P  to satisfy the constraint. 
We first define p̅ = P_1/d, which is the average value of the pre-normalization probabilistic mask. 
Note that 1-p̅ is the average value of 1-P_1.
We define the normalization layer as:



    N_α(P)= 
        α/p̅P,  if p̅≥α
    
       1 - 1-α/1-p̅ (1 - P),   otherwise.


It can be shown that Equation (<ref>) yields N_α(P)∈ [0,1]^d and N_α(P)_1/d = α.

This normalization trick allows us to convert the constrained problem of Equation (<ref>) to the following unconstrained problem:

    min_O, θ∑_j=1^n 1/K∑_k=1^K 
    
     A_θ (F^Hdiag(σ_s(N_α(σ_t(O))-U^(k))) F x_j) - x_j _2^2.




Figure <ref> illustrates an instantiation of the LOUPE model. 
In the depicted version of LOUPE, we used a U-Net architecture to implement A_θ, the reconstruction network. 
We have also experimented with an alternative reconstruction model, namely Cascade-Net proposed in <cit.>.
In the supplementary material, we present LOUPE results obtained with this alternative implementation. 



The LOUPE network accepts complex 2D images represented as two channels, corresponding to the imaginary and real components. 

The network applies a forward discrete Fourier transform, F, converting the data into k-space. 
The k-space measurements are under-sampled by m̃ = σ_s(P-u^(k)), which approximates a Monte Carlo realization of the probabilistic mask P.
The values of the probabilistic mask are computed using the sigmoid operation of pixel-wise parameters O that are learned during the training process. 
The mask m̃ is multiplied element-wise with the k-space data (approximating retrospective under-sampling by inserting zeros at missing points), which is then passed to an inverse discrete Fourier transform F^H.


This image, which will typically contain aliasing artifacts, is then provided as input to the anti-aliasing network A_θ(·),  also  called the reconstruction network.
This network takes a complex-valued under-sampled image, represented as two-channels, and aims to minimize the reconstruction error, such as the squared difference between the magnitude images of ground truth and the reconstruction. 
We can view the entire pipeline to be made up of two building blocks: the first optimizing the under-sampling pattern, and the second solving the reconstruction problem.



The model was implemented in Keras <cit.>, with Tensorflow <cit.> in the back-end. 
Custom functions were adapted from the Neuron library <cit.>. ADAM <cit.> with an initial learning rate of 0.001 was used for optimization, and learning was terminated when validation loss plateaued.
Consistent with prior work <cit.>, a single Monte Carlo realization was drawn for each training datapoint. 
This is common practice in variational neural networks as it is a computationally efficient approach that yields an unbiased estimate of the gradient that is used in stochastic gradient descent. 
We used a batch size of 16.

 We used a slope value of s=200 for σ_s that approximates the thresholding operation and a slope value of t=5 for σ_t that squashes the values of O to the range [0, 1]. 
 These values for the slope hyper-parameters were chosen based on a grid search strategy, where we identified the pair of values that yields the smallest validation loss (see Supplementary Material).
We note, however, that we did not observe a significant amount of sensitivity to these values.

Our code is freely available at: <https://github.com/cagladbahadir/LOUPE/>






§ EMPIRICAL ANALYSIS
 



 §.§ Data


We conducted our experiments using the NYU fastMRI dataset <cit.> (fastmri.med.nyu.edu). 


This is a freely available, large-scale, public data set of knee MRI scans.
The dataset originally comprises of 2D coronal knee scans acquired with two pulse sequences that result with Proton Density (PD) and Proton Density Fat Supressed (PDFS) weighted images. 

In our experiments, we used the provided emulated single-coil (ESC) k-space data from the Biograph mMR (3T) scanner, which were derived from raw 15-channel multi-coil data. 
We used 100 volumes from the provided official training dataset as our training dataset, and split the provided validation data into two halves, where we used 10 volumes for validation and 10 volumes for test.
We could not rely on the provided test data for evaluation as they were not fully sampled. 
For the scans we analyzed, the sequence parameters were: echo train length of 4, matrix size of 320x320, in plane resolution of 0.5 mm × 0.5 mm, slice thickness of 3 mm and no gap between slices. The time of repetition (TR) varied between 2200 and 3000 ms and the Echo Time (TE) ranged between 27 and 34 ms. Training volumes had 38 ± 4 slices, where the validation volumes had 37± 3 and test volumes had 38± 4. 

Each set (training, validation and test) had differing slice sizes across volumes. 
After taking the Inverse Fourier Transform of the ESC k-space data, and rotating and flipping the images to match the orientation in the fastMRI paper <cit.>, we cropped the central 320x320 and normalized by dividing to the maximum magnitude within each volume.

The U-Net architecture we employed was also used in the original NYU fastMRI study <cit.> for reconstruction, with the only difference that their implementation accepted a single-channel, real-valued input image. 






 §.§ Evaluation


We used three metrics in the quantitative evaluations of the reconstructions with respect to the ground truth images created from fully sampled measurements: (1) peak signal to noise ratio (PSNR), (2) structural similarity index (SSIM), and (3) high frequency error norm (HFEN).

PSNR is a widely used metric for evaluating the quality of reconstructions in compressed sensing applications <cit.>, 
defined as:


    PSNR(x, x̂) = 10 log_10max(x)^2 d/x - x̂_2^2,


where d is, as above, the total number of pixels in the full-resolution grid, and x̂ denotes the reconstruction of the ground truth full-resolution image x.


SSIM aims to quantify the perceived image quality <cit.>:


    SSIM(x,x̂)=(2μ_xμ_x̂+c_1)+(2σ_xx̂+c_2)/(μ_x^2+μ_x̂^2+c_1) (σ_x^2+σ_x̂^2+c_2),


where μ_x and μ_x̂ are defined as the “local” average values for the original and reconstruction images, respectively,
computed within an N × N neighborhood.
Similarly, σ_x^2 and σ_x̂^2 are the local variances; and σ_xx̂ is local covariance between the reconstruction and ground truth images.
c_1=(k_1L)^2 and c_2=(k_2L)^2 are  constants that numerically stabilize the division, where L is the dynamic range of the pixel values and k_1 and k_2 are user defined.
We use the parameters provided in <cit.> for computing SSIM values: k_1=0.01, k_2=0.03, and a window size of 7×7.
The dynamic range of pixel values was computed over the entire volume.







Finally, we also report the high-frequency error norm (HFEN), which is used to quantify the quality of reconstruction of edges and fine features.
Following  <cit.>, we use a 15 × 15 Laplacian of Gaussian (LoG) filter with a standard deviation of 1.5 pixels.
The HFEN is then computed as the  L2 difference between LoG filtered ground truth and reconstruction images.



 §.§ Benchmarks
 
We compared the LOUPE-optimized mask and other mask configurations, together with one machine learning (ML) based and three non-ML based  reconstruction techniques. 



  §.§.§ Reconstruction Methods


Each reconstruction method was optimized on a representative slice from the data-set to find the values of hyper-parameters that yielded the best quality reconstruction for each of the individual masks. In other words, hyper-parameter tuning was done separately for each mask.


The first benchmark reconstruction method is BM3D <cit.>, which is an iterative algorithm that alternates between de-noising and reconstructions steps, and was shown to yield faithful reconstructions in the under-sampled cases.[We used the code at: <http://web.itu.edu.tr/eksioglue/pubs/BM3D_MRI.htm>.] 
A grid search on the parameters: σ, noise and λ were conducted, over the range 3 - 3× 10^2 and 0-10^2, respectively. 

The second reconstruction benchmark method, called P-LORAKS, is based on low-rank modeling of local k-space neighborhoods <cit.>.[ 
We used the code available at: <https://mr.usc.edu/download/loraks2/>.]

A grid search on the parameters: λ, VCC (Virtual Conjugate Coils) and rank was conducted. The regularization parameter λ was searched between 10^-2 and 1. The rank value which is related to the non-convex regularization penalty was searched between 1 to 45 and the VCC, which is a Boolean parameter was searched for both cases: 1 and 0. 

The third benchmark method is Total Generalized Variation (TGV) based reconstruction with Shearlet Transform <cit.>. The method regularizes image regions while keeping the edge information intact.[We used the code at <http://www.math.ucla.edu/ wotaoyin/papers/tgv_shearlet.html>.]
We conducted a grid search for the parameters β, λ, α_0, α_1 and  μ_1, adopting a range for each parameter based on the suggestions from the original paper <cit.>.

λ was searched between 10^-3 and 5 × 10^-1, β was searched between 10^2 and 10^3. The parameters μ_1, α_0 and α_1 were respectively searched between values: 3×10^2 and 10 ×3, 8 × 10^-4 to 10^-2 and 10^-3 to 10^-1.

The final benchmark reconstruction method is the residual U-Net <cit.>, which is also a building block in our model. 
In this framework, the U-Net is used as an anti-aliasing neural network widely used for biomedical image applications. 
Unlike in LOUPE, the benchmark U-Net implementation is trained for a fixed under-sampling mask, that is provided by the user. 
In our experiments the U-Net models were individually trained and tested for each mask configuration. 








  §.§.§ Under-sampling Masks


We implemented several widely used benchmark under-sampling masks, shown in Figure <ref>. 
The first set of masks we considered are what we refer to as 2D Cartesian under-sampling masks.
These are physically feasible for 2D imaging with (i) two phase-encoding dimensions and no frequency-encoding gradient; or (ii) for 3D acquisition with one frequency encoding dimension and two phase-encoding dimensions.
In this category, we have “Random Uniform” <cit.>, which exhibits the benefits of stochastic under-sampling in creating incoherent, noise-like artifacts <cit.>, thus facilitating the discrimination of artifacts from the signal present in the image. 
We also employed a parametric Variable Density <cit.> (VD) mask that samples lower frequencies more densely, while still enjoying the benefits of creating noise-like artifacts.
We used the publicly available code of <cit.> to identify the optimal values of the VD mask that maximizes incoherence (of aliasing artifacts) in the wavelet domain.
Another 2D under-sampling mask we implemented was a data-driven strategy that was based on <cit.>. 
In our implementation, we averaged the magnitude spectrum of all the fully-sampled training data.
The average spectrum was then thresholded to keep the k-space points that have the largest average magnitude. 
We refer to this mask as “spectrum-based.”



Next, we considered 1D under-sampling masks, which are physically feasible for 2D imaging with one frequency-encoding dimension and one phase-encoding dimension.
Cartesian <cit.> with skipped lines orthogonal to the provided read-out direction is a popular mask in this category.
In addition, we implemented a version of LOUPE where the mask was constrained to be made up of lines along the read-out direction, which is similar to a strategy proposed in <cit.>.
We achieved this with a simple modification in our code by sharing the weights corresponding to the entries of O along each read-out line.
This “line-constrained version” of LOUPE, we believe, is one step closer to being physically realistic.
When we need to be explicit, we refer to the first LOUPE version as unconstrained.


We underscore that it might not be fair to compare the 1D and 2D under-sampling masks, as the two categories would rely on different acquisition protocols.

Thus, any comparisons we present below need to account for this difference.











 §.§ Results


Figure <ref> shows the LOUPE optimized and benchmark masks, for two different acceleration rates R=4 and R=8.



The LOUPE,  spectrum-based, and variable density masks share the behavior of a drop in density from lower to higher frequencies, however differ significantly in their symmetry and shape. 
Although the unconstrained LOUPE-optimized masks are similar to the VD masks in terms of emphasizing lower frequencies, importantly lateral frequencies are favored significantly more than ventral/dorsal frequencies. 
The data-driven spectrum-based masks, on the other hand, exhibit a similar tendency to pick out more lateral frequencies, yet with a dramatically different overall shape that concentrates around the two axes.
The line-constrained LOUPE masks also show a strong preference for lower frequencies.



In an earlier conference paper <cit.>, we had reported an unconstrained LOUPE-optimized mask for T1-weighted brain MRI scans, obtained from a public dataset <cit.>.
Figure <ref> shows a side by side comparison of the two LOUPE-optimized masks for the knee and brain anatomies, respectively, emphasizing the asymmetry present in the knee dataset.  The knee mask favors lateral frequencies more than ventral/dorsal frequencies due to the unique features of knee anatomy, where there is significantly more tissue contrast in the lateral direction. Brain scans, on the other hand, exhibit a more radially symmetric behavior. This comparison highlights the importance of a data-driven approach in identifying the under-sampling mask.


















Figure <ref> shows subject-level quantitative reconstruction quality values computed for the different masks and reconstruction methods, under two different acceleration conditions.
We present the results for the 6 test subjects with PD weighted images and 4 test subjects with PDFS weighted images, separately. 
The fat suppression operation inherently lowers the signal level, as fat has the highest levels of signal in an MRI scan, thus yielding a noisier image, where small details are more apparent. 
In contrast, regular PD weighted scans that contain fat tissue have inherently higher SNR. 
Therefore, the PD scans yield better quality reconstructions compared to the PDFS scans. 

Figure <ref> reveals that overall LOUPE yields significantly better results than competing masks, when used with any of the tested reconstruction methods and under the two acceleration rates. 
Furthermore, the U-Net reconstruction model, coupled with LOUPE-optimized masks yields the best reconstructions, for both acceleration rates and for all test subjects.
Unsurprisingly, the unconstrained version of LOUPE is generally better than the line-constrained version, with the exception of the TGV reconstruction method, where the line-constrained LOUPE outperformed unconstrained LOUPE. 


These results suggest that LOUPE-optimized masks can offer a performance gain even when not used with the U-Net based reconstruction method.
We further observe that the Variable Density  and spectrum-based masks outperforms the uniform mask, which is consistent with prior literature. 




Table <ref> lists the summary quantitative values for each of the reconstruction method and mask configurations, where averages within the PD and PDFS subjects are presented separately.
Supplementary Tables 1 and 2 includes standard deviations for these results.
The lowest HFEN and highest PSNR and SSIM values are bolded for each reconstruction method, metric and subject group. 
We observe that the LOUPE-optimized masks yield the highest reconstruction quality under all considered conditions, and the line-constrained and unconstrained versions are often quite similar. 
For reference, we also added the best values reported in <cit.> that was achieved with a U-Net reconstruction model and a Variable Density Cartesian mask. 
We underscore that their model was trained on the full training set and evaluated on the full validation set. 
Nevertheless, we observe that their results are consistent with ours, often falling between our U-Net models trained for Cartesian and VD masks, but always under-performing compared to the U-Net model trained with LOUPE masks.



































































Figures <ref> and <ref> show example PD weighted slices and U-Net based reconstructions, obtained after under-sampling with R=4 and R=8, using different masks.

As can be appreciated from these figures, much of the structural details that correspond to the femur, tibia, tendons, ligaments, muscle tissue, and meniscus are more faithfully captured in the reconstructions obtained with the LOUPE-optimized masks, whereas reconstructions computed with other masks were more blurry.
 
 For example, 
 edges of the posterior cruciate ligament and lateral meniscus suffer from blurring artifacts in the benchmark mask configurations compared to the LOUPE-optimized masks. 
 
 Overall, in agreement with the quantitative results presented above, the LOUPE-optimized mask yields the best-looking reconstruction results and these reconstructions appear less prone to artifacts and blurring, while preserving more of the anatomical detail. 


[t]
Quantitative results of reconstruction quality, averaged over test cases. 1D masks are separated from 2D masks and comparisons should be done accordingly. MSE: Mean squared error. MAE: Mean absolute error. HFEN: High frequency error norm. PSRN: Peak signal to noise ratio. SSIM: Structural similarity metric.  

Acceleration        Method                      Masks                          MSE PD               MSE PDFS             MAE PD              MAE PDFS            HFEN PD             HFEN PDFS           PSNR PD            PSNR PDFS          SSIM PD           SSIM PDFS     
 25*4-fold     6*BM3D       Cartesian  (1D) Skipped Lines                    0.00796              0.00629              0.0646              0.0587              0.9484              1.0596              21.1               22.1               0.59              0.63          

                                                         LOUPE Line Constrained         0.00023              0.00089     0.0105              0.0234     0.548               0.7096              36.69              30.5      0.95              0.82 

                         3-13 
                                                         Uniform Random                 0.01074              0.00902              0.0707              0.0695              0.9467              1.1416              19.95              20.47              0.63              0.58          

                                                         Variable Density               0.00032              0.00137              0.0128              0.029               0.5772              0.8288              35.2               28.66              0.93              0.75          

                                                         Spectrum                       0.00021              0.00135              0.0107              0.0287              0.537               0.8464              37.05              28.72              0.94              0.75          

                                                         LOUPE Unconstrained            0.00014     0.00094              0.0088     0.0241              0.2972     0.5298     38.81     30.28              0.96     0.81          
 2-13 
                             6*LORAKS     Cartesian   (1D) Skipped Lines                   0.22717              0.26113              0.1785              0.1896              2.8759              3.5001              8.45               6.93               0.46              0.28          

                                                         LOUPE Line Constrained         0.00025              0.00097     0.0112              0.0246     0.5801              0.7668              36.32              30.14     0.94     0.8  

                          3-13 
                                                         Uniform Random                 0.01692              0.01861              0.0913              0.1052              1.1886              1.2677              17.99              17.3               0.42              0.26          

                                                         Variable Density               0.00041              0.00231              0.0152              0.0372              0.6364              1.1334              34.1               26.4               0.89              0.64          

                                                         Spectrum                       0.00036              0.00208              0.0142              0.0355              0.6691              1.184               34.69              26.85              0.91              0.67          

                                                         LOUPE Unconstrained            0.00019     0.00113              0.0105     0.0266              0.3754     0.6854     37.42     29.48              0.94     0.78          
 2-13 
                             6*TV         Cartesian   (1D) Skipped Lines                   0.00789              0.00582              0.0643              0.0564              0.8938              0.9126              21.14              22.46              0.55              0.55          

                                                         LOUPE Line Constrained         0.00089     0.0012      0.0217     0.027      0.9029              0.7823              30.57     29.22     0.8      0.73 

                          3-13 
                                                         Uniform Random                 0.01286              0.01178              0.0755              0.0818              0.9303              0.9184              19.33              19.33              0.62              0.55          

                                                         Variable Density               0.003                0.00333              0.0417              0.0454              1.2104              1.0221              25.26              24.78              0.6               0.56          

                                                         Spectrum                       0.00144              0.00188              0.0289              0.0344              0.8695              0.8492              28.48              27.27              0.71              0.62          

                                                         LOUPE Unconstrained            0.00196              0.00218              0.0332              0.0366              0.8991     0.7559     27.13              26.63              0.71              0.65          
 2-13 
                             6*U-Net      Cartesian   (1D) Skipped Lines                   0.00406              0.00345              0.0381              0.0401              0.8265              0.8648              24.11              24.65              0.83              0.79          

                                                         LOUPE Line Constrained         0.00015              0.00058              0.0084              0.0186              0.3563              0.4567              38.53              32.36              0.97              0.89          

                          3-13 
                                                         Uniform Random                 0.0022               0.00204              0.0303              0.0323              0.8342              0.7689              26.89              26.93              0.87              0.85          

                                                         Variable Density               0.00014              0.00053              0.0086              0.018               0.4088              0.5186              38.79              32.74              0.97              0.9           

                                                         Spectrum                       0.00013              0.00053              0.0084              0.0179              0.4266              0.5727              39.05              32.77              0.97              0.9           

                                                         LOUPE Unconstrained            0.00009     0.00048     0.0071     0.0173     0.1681     0.246      40.71     33.16     0.98     0.91 
 2-13 
                             U-Net NYU                   Variable Density Cartesian                                                                                                                               34.01              29.95              0.82              0.64          
 25*8-fold     6*BM3D       Cartesian   (1D) Skipped Lines                   0.00871              0.00671              0.0665              0.0603              0.9855              1.0604              20.71              21.84              0.53              0.6           

                                                         LOUPE Line Constrained         0.00059              0.00124              0.0152              0.0273              0.7843              0.8844              32.66              29.08              0.9               0.75          

                          3-13 
                                                         Uniform Random                 0.01248              0.00854              0.0759              0.0683              1.0331              1.1762              19.56              20.81              0.59              0.5           

                                                         Variable Density               0.00097              0.00153              0.0217              0.0306              0.706               0.9147              30.45              28.16              0.87              0.73          

                                                         Spectrum                       0.00032              0.00148              0.0128              0.0301              0.6758              0.9278              35.26              28.32              0.92              0.7           

                                                         LOUPE Unconstrained            0.0002      0.00104     0.0103     0.0256     0.4684     0.7022     37.16     29.85     0.93     0.76 
 2-13 
                             6*LORAKS     Cartesian    (1D) Skipped Lines                  0.4987               0.34232              0.2875              0.2739              4.2092              4.1481              4.05               4.77               0.44              0.31          

                                                         LOUPE Line Constrained         0.00046              0.00138              0.0147              0.029               0.8339              0.9979              33.56              28.63              0.9               0.73          

                          3-13 
                                                         Uniform Random                 1.15118              0.86483              0.3866              0.3221              9.7904              6.5299              0.26               0.77               0.35              0.17          

                                                         Variable Density               0.00041              0.00188              0.0152              0.0343              0.7671              1.2851              34.12              27.27              0.89              0.65          

                                                         Spectrum                       0.00038              0.00175              0.0145              0.033               0.7555              1.1352              34.45              27.59              0.9               0.66          

                                                         LOUPE Unconstrained            0.00022     0.00111     0.0111     0.0266     0.5038     0.8037     36.79     29.56     0.92     0.75 
 2-13 
                             6*TV         Cartesian     (1D) Skipped Lines                 0.0086               0.00631              0.0664              0.0585              0.9508     0.9606     20.76              22.11              0.46              0.45          

                                                         LOUPE Line Constrained         0.00114     0.00153     0.0245     0.0302     1.1597              1.0106              29.48     28.15     0.79     0.71 

                          3-13 
                                                         Uniform Random                 0.01392              0.01097              0.0792              0.0785              0.975               0.9729              19.15              19.7               0.6               0.55          

                                                         Variable Density               0.00286              0.0032               0.0414              0.0445              1.2487              1.076               25.48              24.95              0.58              0.56          

                                                         Spectrum                       0.00199              0.0024               0.034               0.0388              1.0703              0.9868              27.07              26.21              0.66              0.59          

                                                         LOUPE Unconstrained            0.00229              0.0027               0.0364              0.0407              1.2016              1.0351              26.44              25.69              0.66              0.57          
 2-13 
                             6*U-Net      Cartesian   (1D) Skipped Lines                   0.007                0.00509              0.0517              0.0489              1.0244              0.988               21.72              22.97              0.73              0.68          

                                                         LOUPE Line Constrained         0.00027              0.00079              0.0108              0.0214              0.6327              0.7301              36.0               31.01              0.95              0.85          

                          3-13 
                                                         Uniform Random                 0.00316              0.00279              0.0363              0.0376              0.8189              0.8267              25.17              25.56              0.84              0.8           

                                                         Variable Density               0.00027              0.00074              0.0118              0.0211              0.4746              0.5759              35.82              31.32              0.96              0.87          

                                                         Spectrum                       0.00023              0.00075              0.0105              0.021               0.5828              0.7233              36.69              31.28              0.96              0.86          

                                                         LOUPE Unconstrained            0.00014     0.00066     0.0085     0.0199     0.3044     0.4303     38.85     31.83     0.97     0.87 
 2-13 
                             U-Net NYU                                                                                                                                                                            31.5               28.71              0.76              0.56      







In the Supplementary Material, we present our results obtained with a version of LOUPE where the U-Net reconstruction model was replaced with Cascade-Net <cit.> with the same architecture details of the reference paper.
Everything else in LOUPE, including all optimization settings were kept the same.
We observe that the optimized masks are very similar between U-Net-LOUPE and Cascade-Net-LOUPE. 
We also quantified the quality of reconstructions with these masks and found that there was no substantial difference. 
These results, we believe, demonstrate that LOUPE can be relatively robust to the choice of the reconstruction architecture. 
We note, however, that we have not experimented with more advanced reconstruction networks, such as the recently proposed <cit.>, which might yield different results.
We consider the exploration of different architectural designs for LOUPE as an important direction for future research.




§ DISCUSSION



Compressed sensing is a promising approach to accelerate MRI and make the technology more accessible and affordable. 
Since its introduction over a decade ago, it has gained in popularity and been used in a range of application domains.
Compressed sensing MRI has two fundamental challenges. 
The first is identifying where in k-space to sample (acquire) from. 

One approach to solve this problem is to take a data-driven strategy, where the optimal under-sampling pattern is identified using, for example, a collection of full-resolution scans that are retrospectively under-sampled.


The second core problem of compressed sensing MRI is solving the ill-posed inverse problem of reconstructing a high quality image from under-sampled measurements.
This is conventionally solved via a regularized optimization approach, and more recently deep learning techniques have been proposed to achieve superior performance.

In this paper, we tackled both of these problems simultaneously.
We leveraged recently proposed deep learning based reconstruction techniques to design an end-to-end learning technique that enabled us to optimize the under-sampling pattern on given training data.
Our approach relies on full-resolution data that is under-sampled retrospectively, which is then fed to a reconstruction network.
The final output is then compared to the original input to assess overall quality.
We formulate an objective that optimizes reconstruction quality with a constrained number of collected measurements.
The stronger the constraint, the more aggressive acceleration rates we can achieve.

Our results indicate that LOUPE, the proposed method, can compute an optimized under-sampling pattern that is data-dependent and performs well with different reconstruction methods.
In this paper, we present empirical results obtained with a large-scale, publicly available knee MRI dataset.
Across all conditions we experimented with,  LOUPE-optimized masks coupled with the U-Net based reconstruction model offered the most superior reconstruction quality.
Even with an aggressive 8-fold acceleration rate, LOUPE's reconstructions contained much of the anatomical detail that was missed by alternative masks as seen in Figure <ref>.

In our primary experiments, we used a U-Net based reconstruction method in the implementation of LOUPE. 
Additionally, we investigated the use of an alternative reconstruction network, namely Cascade-Net. 
Our results, which we present in the Supplementary Material, indicate that both versions of LOUPE achieved similar optimized masks and yielded similar reconstruction quality.
However, we consider the exploration of different architectural designs (such as ADMM-Net <cit.> or MoDL <cit.>) for LOUPE as an important direction for future research.

The LOUPE-optimized under-sampling mask for the analyzed knee MRI scans exhibited an interesting asymmetric pattern, which emphasized higher frequencies along the lateral direction much more than those along the vertical direction. 

This is likely due the anisotropy in the field of view of sagittal knee images.
This is in stark contrast to the LOUPE-optimized mask that we recently presented at a conference for T1-weighted brain MRI scans <cit.>.
The brain mask exhibited a radially symmetric nature, very similar to a variable density mask, which is widely used in compressed sensing MRI.
This comparison highlights the importance of adapting the under-sampling pattern to the data and application at hand.

An important caveat of the presented LOUPE framework is that we largely ignored the costs of physically implementing the under-sampling pattern in an MR pulse sequence.
In  the unconstrained version of LOUPE, the cost of an under-sampling pattern is merely captured by the number of collected measurements.
We underscore that in 3D, the LOUPE-optimized mask can be implemented as a Cartesian trajectory, since the isolated measurement points in the coronal slices line up along the z-direction. 
We also presented results for a readout-line-constrained version of LOUPE, which is consistent with widely used 2D acquisition protocols.
These results, we believe, demonstrate the utility and flexibility of LOUPE.
Nonetheless, integrating the actual physical cost of a specific under-sampling pattern in an MR pulse sequence is a direction that we leave for future research.
One way to achieve this would be to constrain the sub-sampling mask to a class of viable trajectories, for example, by directly parameterizing those trajectories.
Alternatively, one can incorporate sampling constraints that obey a set of predefined hardware requirements (as defined in terms of, e.g., peak currents and maximum slew rates of magnetic gradients), as described recently in <cit.>.

Our paper restricted its treatment to Cartesian sampling and largely ignored non-Cartesian (off-the-grid) acquisition schemes, which can offer important benefits <cit.>.
Extension of LOUPE to non-Cartesian settings will involve the use of the non-uniform (inverse) Fourier transform operator or an approximation of it.
Furthermore, one would need to interpolate off-grid measurements in retrospective sampling.


Another weakness of the presented LOUPE implementation is due to the relaxation of the thresholding operation, which allowed us to use
a back-propagation-based learning strategy.
However, due to this relaxation, the reconstruction model sees input data that have been rescaled with a continuous value between 0 and 1. 
By adopting a relatively large slope parameter (s=200), we ensure that these continuous values are almost always very close to zero or 1. 
However, this relaxation means that once we have obtained the optimized mask, we need to retrain a reconstruction model with a binary mask. Furthermore, we are likely paying a performance penalty due to the relaxation, since it will undoubtedly influence the optimization. 
We see two ways to address this issue. 
The first approach is to implement the so-called “straight-through” gradient estimation strategy, where the continuous relaxation is merely used in the backward pass for computing the gradient and not in the forward pass <cit.>. 
Another approach would involve using non back-propagation based techniques, such as the REINFORCE estimator <cit.>. 
This technique is known to suffer from high variance and there are several recent papers that have attempted to address this issue <cit.>. 
We plan to explore this direction in future work.




In this paper, we used an L2 difference between the magnitude of the ground truth and reconstruction images to train LOUPE.
Our experiments with an L2-loss computed on two-channel (complex-valued) reconstructions produced inferior results.
We also noticed that our LOUPE-optimized masks did not emphasize high frequency content, which can also explain some of the blurriness in reconstructions. 
This is particularly noticeable in the line-constrained LOUPE results. 
We believe that our choice of the L2 loss is the main cause of this effect, which can be remedied with alternate loss functions that care more about high-frequency content. 
We note that we experimented with L1 reconstruction loss, as we reported in our conference paper <cit.>.
 Our experiments (results not shown) revealed that L1 loss can yield better quality results, particularly as measured by the high frequency error norm (HFEN).





We intend to explore such differences in future work.
While these loss functions (L1 or L2) are a widely used global metrics of reconstruction quality, they might miss clinically important anatomical details.
An alternative approach could be to devise loss functions that emphasize structural features that are relevant to the clinical application.
This is another future direction we will explore.

In the current version of LOUPE, we did not consider parallel imaging via multiple coils.
Such hardware-based techniques, in combination with compressed sensing, promise to yield even higher degrees of acceleration.
Investigating how to combine LOUPE with multi-coil imaging will be an important area of research, as was recently demonstrated <cit.>.



Finally, we would like to emphasize that LOUPE, as in any other data-driven optimization strategy, will perform sub-optimally if test time data differ significantly from training data.
For example, if we train only on healthy subjects and then, at test time, are presented with pathological cases, the reconstructed scans might not be clinically useful.
Thus, carefully monitoring reconstruction quality and ensuring the model is adapted to the new data will be of utmost importance in deploying this approach in the real world.
For example, this could be achieved via fine-tuning the pre-trained LOUPE model on new data, as recently described in <cit.>.
We also consider machine learning paradigms that don't rely on high quality full-resolution data as an important direction. 
There have been some recent pre-prints in this domain <cit.>.
Integrating such techniques into LOUPE will be critical for certain applications.
For instance, in some anatomical regions such as the abdomen, a fully sampled scenario may not be possible due to organ motion, breathing, the cardiac cycle or other factors. 
This is also true for populations such as young children or patient groups, where subject motion can make the acquisition of high quality fully sampled data impossible. 











§ ACKNOWLEDGMENT


This work was, in part, supported by NIH R01 grants (R01LM012719 and R01AG053949, to MRS), the NSF NeuroNex grant (1707312, to MRS), and an NSF CAREER grant (1748377, to MRS). This work was also in part supported by a Fulbright Scholarship (to CDB).

Data used in the preparation of this article were obtained from the NYU fastMRI Initiative database <cit.>. As such, NYU fastMRI investigators provided data but did not participate in analysis or writing of this report. A listing of NYU fastMRI investigators, subject to updates, can be found at:<fastmri.med.nyu.edu>. The primary goal of fastMRI is to test whether machine learning can aid in the reconstruction of medical images.

The concepts and information presented in this paper are based on research results that are not commercially available.

Finally, we would like to express our gratitude for the extremely thorough, thoughtful, and constructive feedback we received from the Associate Editor Justin Haldar, and three anonymous reviewers. 




























IEEEtran







§ SUPPLEMENTARY MATERIAL








 §.§ Fine-tuning the slope hyperparameters


Figure <ref> shows the validation set MAE and MSE for LOUPE with masks optimized using a variety of hyper-parameter values for the two sigmoidal slopes s and t.
These errors are relatively stable compared to the variation of results among the baselines in Table I. 
For example, for MAE, changing s by ten-fold and t by 2-fold results in a loss variation of less than ±0.02, which is well within the variability of the results produced by the benchmark masks. 




 §.§ Details of the Cascade-Net Implementation

We implemented the identical architecture described in the original Cascade-Net paper. 
Specifically, for the regularization layers of Cascade-Net, a 5-layer model was used with channel size 64 at each layer. Each layer consists of convolution followed by a ReLU activation function. We used a residual learning strategy which adds the zero-filled input to the output of the CNN. The overall architecture is unrolled such that K=5. In the original paper and in our case, lambda is not learnable but is a hyperparameter that we set to 0.001, which was found by grid search over validation loss and further optimized using a Bayesian Optimization hyper-parameter tuning package: <https://iopscience.iop.org/article/10.1088/1749-4699/8/1/014008/meta>.














[t]
Supplementary Table, Mean Values ± Standard Deviations  

Acceleration                 Method     Masks                      MSE PD                MSE PDFS              MAE PD              MAE PDFS            HFEN PD             HFEN PDFS           PSNR PD          PSNR PDFS        SSIM PD         SSIM PDFS   
 24*4-fold     BM3D       Cartesian                  0.00796 ± 0.00176     0.00629 ± 0.00127     0.0646 ± 0.0092     0.0587 ± 0.0051     0.9484 ± 0.0091     1.0596 ± 0.0199     21.1 ± 0.99      22.1 ± 0.9       0.59 ± 0.03     0.63 ± 0.01 

                             BM3D       LOUPE Line Constrained     0.00023 ± 8e-05       0.00089 ± 7e-05       0.0105 ± 0.0018     0.0234 ± 0.001      0.548 ± 0.0421      0.7096 ± 0.0715     36.69 ± 1.62     30.5 ± 0.36      0.95 ± 0.02     0.82 ± 0.02 

                             BM3D       Uniform Random             0.01074 ± 0.0035      0.00902 ± 0.00081     0.0707 ± 0.0176     0.0695 ± 0.0063     0.9467 ± 0.0738     1.1416 ± 0.1027     19.95 ± 1.56     20.47 ± 0.39     0.63 ± 0.06     0.58 ± 0.02 

                             BM3D       Variable Density           0.00032 ± 0.00011     0.00137 ± 0.00013     0.0128 ± 0.0023     0.029 ± 0.0013      0.5772 ± 0.0449     0.8288 ± 0.0846     35.2 ± 1.54      28.66 ± 0.41     0.93 ± 0.02     0.75 ± 0.02 

                             BM3D       Spectrum                   0.00021 ± 9e-05       0.00135 ± 0.00012     0.0107 ± 0.0021     0.0287 ± 0.0013     0.537 ± 0.0464      0.8464 ± 0.0831     37.05 ± 1.78     28.72 ± 0.4      0.94 ± 0.02     0.75 ± 0.02 

                             BM3D       LOUPE Unconstrained        0.00014 ± 5e-05       0.00094 ± 8e-05       0.0088 ± 0.0015     0.0241 ± 0.0011     0.2972 ± 0.049      0.5298 ± 0.0845     38.81 ± 1.49     30.28 ± 0.39     0.96 ± 0.01     0.81 ± 0.02 
 2-13 
                             LORAKS     Cartesian                  0.22717 ± 0.25121     0.26113 ± 0.13396     0.1785 ± 0.1043     0.1896 ± 0.0521     2.8759 ± 1.0952     3.5001 ± 0.9479     8.45 ± 3.86      6.93 ± 3.64      0.46 ± 0.07     0.28 ± 0.06 

                             LORAKS     LOUPE Line Constrained     0.00025 ± 7e-05       0.00097 ± 0.0001      0.0112 ± 0.0018     0.0246 ± 0.0013     0.5801 ± 0.0295     0.7668 ± 0.0788     36.32 ± 1.41     30.14 ± 0.43     0.94 ± 0.02     0.8 ± 0.02  

                             LORAKS     Uniform Random             0.01692 ± 0.00581     0.01861 ± 0.0002      0.0913 ± 0.0231     0.1052 ± 0.0062     1.1886 ± 0.0497     1.2677 ± 0.0434     17.99 ± 1.57     17.3 ± 0.05      0.42 ± 0.08     0.26 ± 0.02 

                             LORAKS     Variable Density           0.00041 ± 0.00014     0.00231 ± 0.00028     0.0152 ± 0.0025     0.0372 ± 0.0024     0.6364 ± 0.0881     1.1334 ± 0.1639     34.1 ± 1.45      26.4 ± 0.57      0.89 ± 0.04     0.64 ± 0.04 

                             LORAKS     Spectrum                   0.00036 ± 0.00013     0.00208 ± 0.00021     0.0142 ± 0.0024     0.0355 ± 0.0017     0.6691 ± 0.0743     1.184 ± 0.1094      34.69 ± 1.45     26.85 ± 0.43     0.91 ± 0.03     0.67 ± 0.03 

                             LORAKS     LOUPE Unconstrained        0.00019 ± 6e-05       0.00113 ± 9e-05       0.0105 ± 0.0017     0.0266 ± 0.0011     0.3754 ± 0.0745     0.6854 ± 0.0853     37.42 ± 1.37     29.48 ± 0.35     0.94 ± 0.02     0.78 ± 0.02 
 2-13 
                             TV         Cartesian                  0.00789 ± 0.00176     0.00582 ± 0.00129     0.0643 ± 0.0092     0.0564 ± 0.0055     0.8938 ± 0.0028     0.9126 ± 0.0141     21.14 ± 1.01     22.46 ± 0.99     0.55 ± 0.03     0.55 ± 0.03 

                             TV         LOUPE Line Constrained     0.00089 ± 0.00012     0.0012 ± 2e-05        0.0217 ± 0.0019     0.027 ± 0.0002      0.9029 ± 0.1649     0.7823 ± 0.0636     30.57 ± 0.6      29.22 ± 0.07     0.8 ± 0.03      0.73 ± 0.02 

                             TV         Uniform Random             0.01286 ± 0.00505     0.01178 ± 0.0017      0.0755 ± 0.0229     0.0818 ± 0.01       0.9303 ± 0.0111     0.9184 ± 0.0134     19.33 ± 2.01     19.33 ± 0.61     0.62 ± 0.08     0.55 ± 0.01 

                             TV         Variable Density           0.003 ± 0.00037       0.00333 ± 0.00013     0.0417 ± 0.0041     0.0454 ± 0.0012     1.2104 ± 0.2833     1.0221 ± 0.0692     25.26 ± 0.54     24.78 ± 0.18     0.6 ± 0.04      0.56 ± 0.06 

                             TV         Spectrum                   0.00144 ± 0.00025     0.00188 ± 6e-05       0.0289 ± 0.0032     0.0344 ± 0.0007     0.8695 ± 0.1871     0.8492 ± 0.0663     28.48 ± 0.72     27.27 ± 0.14     0.71 ± 0.06     0.62 ± 0.02 

                             TV         LOUPE Unconstrained        0.00196 ± 0.00032     0.00218 ± 0.00011     0.0332 ± 0.0037     0.0366 ± 0.0011     0.8991 ± 0.2634     0.7559 ± 0.0998     27.13 ± 0.75     26.63 ± 0.22     0.71 ± 0.03     0.65 ± 0.02 
 2-13 
                             U-Net      Cartesian                  0.00406 ± 0.00112     0.00345 ± 0.00038     0.0381 ± 0.0081     0.0401 ± 0.0018     0.8265 ± 0.0395     0.8648 ± 0.0182     24.11 ± 1.37     24.65 ± 0.47     0.83 ± 0.04     0.79 ± 0.01 

                             U-Net      LOUPE Line Constrained     0.00015 ± 5e-05       0.00058 ± 3e-05       0.0084 ± 0.0014     0.0186 ± 0.0004     0.3563 ± 0.0148     0.4567 ± 0.018      38.53 ± 1.52     32.36 ± 0.19     0.97 ± 0.01     0.89 ± 0.0  

                             U-Net      Uniform Random             0.0022 ± 0.00078      0.00204 ± 0.00017     0.0303 ± 0.0078     0.0323 ± 0.0009     0.8342 ± 0.0974     0.7689 ± 0.0149     26.89 ± 1.7      26.93 ± 0.35     0.87 ± 0.03     0.85 ± 0.01 

                             U-Net      Variable Density           0.00014 ± 4e-05       0.00053 ± 3e-05       0.0086 ± 0.0015     0.018 ± 0.0005      0.4088 ± 0.0231     0.5186 ± 0.0158     38.79 ± 1.5      32.74 ± 0.21     0.97 ± 0.01     0.9 ± 0.0   

                             U-Net      Spectrum                   0.00013 ± 4e-05       0.00053 ± 3e-05       0.0084 ± 0.0014     0.0179 ± 0.0005     0.4266 ± 0.02       0.5727 ± 0.019      39.05 ± 1.52     32.77 ± 0.22     0.97 ± 0.01     0.9 ± 0.0   

                             U-Net      LOUPE Unconstrained        9e-05 ± 3e-05         0.00048 ± 1e-05       0.0071 ± 0.0011     0.0173 ± 0.0002     0.1681 ± 0.01       0.246 ± 0.0162      40.71 ± 1.44     33.16 ± 0.09     0.98 ± 0.01     0.91 ± 0.0  
 24*8-fold     BM3D       Cartesian                  0.00871 ± 0.00189     0.00671 ± 0.00146     0.0665 ± 0.0094     0.0603 ± 0.0054     0.9855 ± 0.0069     1.0604 ± 0.012      20.71 ± 0.99     21.84 ± 0.95     0.53 ± 0.02     0.6 ± 0.02  

                             BM3D       LOUPE Line Constrained     0.00059 ± 0.00024     0.00124 ± 0.0001      0.0152 ± 0.0029     0.0273 ± 0.0012     0.7843 ± 0.0395     0.8844 ± 0.0559     32.66 ± 1.7      29.08 ± 0.36     0.9 ± 0.02      0.75 ± 0.02 

                             BM3D       Uniform Random             0.01248 ± 0.00567     0.00854 ± 0.00211     0.0759 ± 0.024      0.0683 ± 0.0112     1.0331 ± 0.0744     1.1762 ± 0.0649     19.56 ± 2.22     20.81 ± 0.98     0.59 ± 0.09     0.5 ± 0.04  

                             BM3D       Variable Density           0.00097 ± 0.00038     0.00153 ± 0.00012     0.0217 ± 0.0045     0.0306 ± 0.0011     0.706 ± 0.0639      0.9147 ± 0.0872     30.45 ± 1.6      28.16 ± 0.34     0.87 ± 0.03     0.73 ± 0.03 

                             BM3D       Spectrum                   0.00032 ± 0.00012     0.00148 ± 0.00012     0.0128 ± 0.0023     0.0301 ± 0.0013     0.6758 ± 0.0463     0.9278 ± 0.071      35.26 ± 1.65     28.32 ± 0.36     0.92 ± 0.02     0.7 ± 0.03  

                             BM3D       LOUPE Unconstrained        0.0002 ± 7e-05        0.00104 ± 9e-05       0.0103 ± 0.0017     0.0256 ± 0.0012     0.4684 ± 0.0444     0.7022 ± 0.0835     37.16 ± 1.52     29.85 ± 0.39     0.93 ± 0.02     0.76 ± 0.03 
 2-13 
                             LORAKS     Cartesian                  0.4987 ± 0.32268      0.34232 ± 0.07458     0.2875 ± 0.0871     0.2739 ± 0.0318     4.2092 ± 1.071      4.1481 ± 0.9071     4.05 ± 3.22      4.77 ± 0.99      0.44 ± 0.06     0.31 ± 0.03 

                             LORAKS     LOUPE Line Constrained     0.00046 ± 0.00013     0.00138 ± 0.00011     0.0147 ± 0.002      0.029 ± 0.0013      0.8339 ± 0.0421     0.9979 ± 0.0701     33.56 ± 1.27     28.63 ± 0.36     0.9 ± 0.02      0.73 ± 0.02 

                             LORAKS     Uniform Random             1.15118 ± 0.80465     0.86483 ± 0.21083     0.3866 ± 0.1535     0.3221 ± 0.0275     9.7904 ± 3.022      6.5299 ± 0.6076     0.26 ± 2.62      0.77 ± 1.13      0.35 ± 0.07     0.17 ± 0.03 

                             LORAKS     Variable Density           0.00041 ± 0.00013     0.00188 ± 0.00014     0.0152 ± 0.0025     0.0343 ± 0.0014     0.7671 ± 0.0731     1.2851 ± 0.1822     34.12 ± 1.5      27.27 ± 0.34     0.89 ± 0.03     0.65 ± 0.04 

                             LORAKS     Spectrum                   0.00038 ± 0.00011     0.00175 ± 0.00015     0.0145 ± 0.0023     0.033 ± 0.0015      0.7555 ± 0.0334     1.1352 ± 0.0907     34.45 ± 1.36     27.59 ± 0.37     0.9 ± 0.03      0.66 ± 0.03 

                             LORAKS     LOUPE Unconstrained        0.00022 ± 7e-05       0.00111 ± 0.00011     0.0111 ± 0.0018     0.0266 ± 0.0014     0.5038 ± 0.054      0.8037 ± 0.0823     36.79 ± 1.41     29.56 ± 0.45     0.92 ± 0.02     0.75 ± 0.02 
 2-13 
                             TV         Cartesian                  0.0086 ± 0.00189      0.00631 ± 0.00143     0.0664 ± 0.0094     0.0585 ± 0.0055     0.9508 ± 0.0016     0.9606 ± 0.008      20.76 ± 1.0      22.11 ± 0.99     0.46 ± 0.04     0.45 ± 0.06 

                             TV         LOUPE Line Constrained     0.00114 ± 0.00018     0.00153 ± 5e-05       0.0245 ± 0.0021     0.0302 ± 0.0004     1.1597 ± 0.244      1.0106 ± 0.058      29.48 ± 0.67     28.15 ± 0.15     0.79 ± 0.03     0.71 ± 0.02 

                             TV         Uniform Random             0.01392 ± 0.00634     0.01097 ± 0.00251     0.0792 ± 0.0274     0.0785 ± 0.0128     0.975 ± 0.0079      0.9729 ± 0.0058     19.15 ± 2.41     19.7 ± 0.91      0.6 ± 0.08      0.55 ± 0.03 

                             TV         Variable Density           0.00286 ± 0.00043     0.0032 ± 0.00012      0.0414 ± 0.0043     0.0445 ± 0.0011     1.2487 ± 0.3041     1.076 ± 0.0697      25.48 ± 0.69     24.95 ± 0.17     0.58 ± 0.05     0.56 ± 0.07 

                             TV         Spectrum                   0.00199 ± 0.00032     0.0024 ± 9e-05        0.034 ± 0.0036      0.0388 ± 0.0008     1.0703 ± 0.1989     0.9868 ± 0.0546     27.07 ± 0.73     26.21 ± 0.17     0.66 ± 0.04     0.59 ± 0.07 

                             TV         LOUPE Unconstrained        0.00229 ± 0.00031     0.0027 ± 0.00012      0.0364 ± 0.0036     0.0407 ± 0.0011     1.2016 ± 0.3756     1.0351 ± 0.1059     26.44 ± 0.61     25.69 ± 0.19     0.66 ± 0.04     0.57 ± 0.01 
 2-13 
                             U-Net      Cartesian                  0.007 ± 0.00191       0.00509 ± 0.00075     0.0517 ± 0.0101     0.0489 ± 0.0023     1.0244 ± 0.0283     0.988 ± 0.0041      21.72 ± 1.26     22.97 ± 0.6      0.73 ± 0.02     0.68 ± 0.02 

                             U-Net      LOUPE Line Constrained     0.00027 ± 9e-05       0.00079 ± 6e-05       0.0108 ± 0.0019     0.0214 ± 0.0007     0.6327 ± 0.0499     0.7301 ± 0.0225     36.0 ± 1.53      31.01 ± 0.3      0.95 ± 0.01     0.85 ± 0.01 

                             U-Net      Uniform Random             0.00316 ± 0.00089     0.00279 ± 0.00026     0.0363 ± 0.0081     0.0376 ± 0.0022     0.8189 ± 0.0472     0.8267 ± 0.0133     25.17 ± 1.24     25.56 ± 0.39     0.84 ± 0.03     0.8 ± 0.02  

                             U-Net      Variable Density           0.00027 ± 8e-05       0.00074 ± 4e-05       0.0118 ± 0.0019     0.0211 ± 0.0006     0.4746 ± 0.0268     0.5759 ± 0.0193     35.82 ± 1.31     31.32 ± 0.24     0.96 ± 0.01     0.87 ± 0.0  

                             U-Net      Spectrum                   0.00023 ± 8e-05       0.00075 ± 5e-05       0.0105 ± 0.0019     0.021 ± 0.0007      0.5828 ± 0.0394     0.7233 ± 0.0199     36.69 ± 1.59     31.28 ± 0.29     0.96 ± 0.01     0.86 ± 0.01 

                             U-Net      LOUPE Unconstrained        0.00014 ± 5e-05       0.00066 ± 2e-05       0.0085 ± 0.0014     0.0199 ± 0.0003     0.3044 ± 0.0175     0.4303 ± 0.0229     38.85 ± 1.59     31.83 ± 0.15     0.97 ± 0.01     0.87 ± 0.0 

















































