\section{Introduction}
\IEEEPARstart{M}{agnetic} Resonance Imaging (MRI) is one kind of widely used medical imaging modalities for clinical diagnosis, which requires a long scan time and has the risk of motion-related artifacts in the reconstructed image. Nyquist sampling theory, though having guided the signal sampling process for many years, has strict sampling requirements and puts a heavy burden on hardware equipment. CS-MRI that only requires much lower sampling rate than Nyquist sampling theory has been proposed to reconstruct image from the sparse characteristics of signals. Although it may cause aliasing artifacts in a spatial domain, the quality of reconstructed image is not significantly reduced.
\IEEEpubidadjcol
Classical CS-MRI methods can learn directly more flexible sparse representation from under-sampled data by restricting solution in a specific transformation domain \cite{Lustig2007} or in a generic dictionary-based subspace \cite{Zhan2016}. Discrete wavelet transform \cite{Lai2016} and discrete cosine transform \cite{Lingala2013} have been used for CS-MRI reconstruction. Due to simple and effective, Total Variation (TV) regularization has been widely used in MRI reconstruction although it introduces staircase artifacts in reconstructed image \cite{Block2007}. Some sophisticated non-local sparsity methods that use groups of local similar patches to exploit the non-local self-similarity properties, can capture more texture priors and improve significantly CS reconstruction performance \cite{Qu2012,Qu2014}. The BM3D denoiser \cite{Metzler2016} has been integrated into CS reconstruction as a new approximate message passing (AMP) framework. Many studies formulate CS a reconstruction problem as the sparsity-regularized optimization problem and then solve it by using different iterative algorithms such as Iterative Shrinkage Threshold Algorithm (ISTA) \cite{Beck2009}, Primal-Dual algorithm \cite{Chambolle2010} and Alternating Direction Multiplier Method (ADMM) \cite{Boyd2010} etc. All the above CS-MRI methods that are based on strongly interpretable and well predefined sparsity image prior, have the advantages of theoretical support and strong convergence. But they usually require expensive computations, but have to face the problems of  selection of suitable  regularizers and model parameters. Consequently, some reconstructed results are not unsatisfactory.


Recently, deep learning (DL) has achieved great success in computer vision community. A generalized inverse problem (e.g., image super-resolution (SR) \cite{Dong2016}, etc.) is exploited to solve CS-MRI reconstruction. A deep Convolutional Neural Network (CNN) has been proposed to reconstruct MRI from under-sampled data \cite{Wang2016b}.
Based on learning a direct inversion of image inverse problem, the network \cite{Jin2017} combines multi-resolution decomposition and residual learning to remove artifacts while preserving image structure. RefineGAN is a variant of fully-residual convolutional auto-encoder and generative adversarial networks (GANs) with cyclic data consistency loss for CS-MRI \cite{Quan2018}. DAGAN couples adversarial loss with an innovative content loss to better preserve texture and edges in CS-MRI, and also incorporates frequency-domain information to enforce similarity in both the image and frequency domains \cite{Yang2018}. W-net can work both in $k$-space and image domain for CS-MRI \cite{Souza2019}, which is composed of a complex-valued residual U-net in $k$-space, an iFFT operation, and a real-valued U-net in image domain. DuDoRNet with deep embedded T1 prior simultaneously recovers $k$-space information and image for accelerating the acquisition of MRI \cite{Zhou2020}. These existing DL-based CS-MRI methods are data-driven based on a large amount of training data  and without any model prior.

Deep unfolding learnable framework of inheriting the merits of model-based and DL-based CS-MRI methods, has sufficient theoretical support and also good performance \cite{Yang2020}. It is first proposed to learn optimal sparse codes in the Learned Iterative Shrinkage-Thresholding Algorithm (LISTA) \cite{Gregor2010}. Later, Yang et al. \cite{Yang2016} presented a novel
deep ADMM-Net framework to supervise data flow graph
in image reconstruction network. A cascaded dilated dense network with two-step data consistency operation in $k$-space is designed for CS-MRI reconstruction \cite{Zheng2019}. Also a conjugate gradient image reconstruction with a CNN-based regularization prior is employed to build MoDL architecture \cite{Aggarwal2019}. Then by unfolding the iterative process of variable splitting optimization scheme, Duan et al. proposed a novel end-to-end trainable deep neural network (DNN) denoted as VS-Net \cite{Duan2019}. DC-CNN \cite{Schlemper2018} using a deep cascade of CNNs with data consistency layers is proposed to reconstruct MRI from under-sampled data. ISTA-Net+ is designed by mapping ISTA into deep CNN framework to learn proximal mapping \cite{Zhang2018}. To embedded FISTA algorithm \cite{Lustig2007} into a deep network, FISTA-Net \cite{Xiang2020} that consists of three-step update blocks including gradient descent, proximal mapping, and acceleration is designed.

Recently, the techniques of solving linear/nonlinear system have been widely used to design effective DNNs. ResNet \cite{He2016} is partly motivated by the hierarchical correction of residual in classical iterative algorithm. PolyNet designs a PolyInception module to enhance feature extraction of network \cite{Zhang2017}. RevNet\cite{Gomez2017} and LM-ResNet \cite{Lu2018} can be interpreted as different reversible Euler-type discrete dynamic systems of ordinary differential equations (ODEs). A deep multigrid method is proposed to optimize restriction and prolongation operations in two-grid scheme, and is straightforwardly extended to the geometric multigrid method \cite{Katrutsa2017}. A multigrid extension of CNNs is proposed to improve accuracy and computational efficiency on CIFAR and ImageNet classification tasks \cite{Ke2017}. MgNet is also designed for image classification to explore the connection between multigrid and CNNs \cite{He2019}. These models can be effective for improvements of deep learning models, and in particular for the mathematical understanding and analysis of network architecture.


However, there are still two main shortcomings: i) the derivation from mathematical theory to network design is not always natural enough for these existing deep unfolding methods. ISTA-Net and FISTA-Net directly replace the nonlinear transformation by several convolution layers, but no reasonable explanation is given. The proximal-point sub-problem has not well been analyzed and explained in DC-CNN, MoDL, and VS-Net. It leads that CNN is directly used to reconstruct the low CS-ratio image;
2) most of existing CS-MRI methods are not flexible enough to handle multi-sampling-ratio reconstruction assignments. For each CS ratio, these methods usually have to train an independent model by considering it as a single task, which results in expensive training cost and large storage space. And this is often contradicted with the fact that real scenarios often contain different CS ratios.

To overcome such drawbacks, inspired by efficient correction of multigrid technique in multiple resolutions, we start from classical CS-MRI optimization problem, and embed it into a two-scale correction-distillation architecture which consists of pre-relaxation, correction and geometric prior distillation. The major differences between this work and the previous work \cite{Fan2021} are that we propose a unified framework to train/predict multi-sampling-ratio CS-MRI tasks in a single model, and two-grid-cycle correction and geometric priors are refined to design network architecture for further promoting reconstruction performance. The main contributions of this work can be summarized as follows
\begin{itemize}
\item[1)] A novel deep unfolding unified framework which enjoys much flexibility to conduct multi-sampling-ratio CS-MRI through a single model is proposed.
\item[2)] A multi-grid inspired   unfolding correction-distillation scheme which can not only
incorporate frequency-domain information to compensate low frequency error in $k$-space, but also learn geometric priors of MR image by adding a geometric prior distillation module. The proposed method gives us a novel perspective to design network architectures.
\item[3)]
A condition module is employed to transmit CS sampling ratio to step-length and to characterize noise level in every stage;
\item[4)] A linkage of  model-based and learning-based methods which integrates the theoretical guarantee of model-based methods and the superior reconstruction performances of DL-based methods, to deliver the expected generalizability of deep learning;
\item[5)] A learning framework to make all physical-model parameters learnable to ensure that suitable choices are automated, as shown in extensive experiments in favourably comparing with state-of-the-art methods in terms of visualizations and quantitative evaluations on flexibility and stability.
\end{itemize}

The rest of paper is organized as follows. In Section II,
the proposed correction-distillation CS-MRI framework is introduced in detail. The
experimental results are shown in Section III. The conclusion of
this paper is presented in Section IV.

\begin{figure*}[t]
\centerline{\includegraphics[width=2\columnwidth]{fig/figs/TCGD_TNNLS_Reconstruction_fig01.pdf}}
\caption{The overall architecture of the proposed unifying multi-sampling-ratio CS-MRI framework with two-grid-cycle correction and geometric prior distillation (CGPD-CSNet). It consists of three major ingredients, i.e. pre-relaxation module $\mathcal{D}_\ell$, correction module $\mathcal{C}_\ell$, and geometric prior distillation module $\mathcal{P}_\ell$.}
\label{fig_architecture}
\end{figure*}


\section{Methodology}

The general sampling equation in Fourier space can be formulated as follows
\begin{equation}
\mathcal{T}\bm{x}=\bm{y},
\label{eq0}
\end{equation}
where $\bm{x} \in \mathbb{R}^{d}$ ($d=mn$) is a target MRI, $\mathcal{T}=PF\in \mathbb{R}^{q\times d} $ is an under-sampled Fourier measurement matrix, $P$ is an under-sampling matrix and $F$ is a discrete Fourier transform, $\bm{y} \in \mathbb{R}^{q} (q\ll d)$ is the under-sampled $k$-space data.

In this section, we aim to design an explainable unifying deep learning CS-MRI reconstruction framework, which is denoted as CGPD-CSNet. Our method inherits the main advantages of model-based and DL-based methods. The overall CNN architecture of the proposed network, which is inspired for learning a unified multi-sampling-ratio CS-MRI optimization problem, is shown in Fig.\ref{fig_architecture}, and more details are provided hereafter.



\subsection{Inverse Problems for Image Reconstruction}
The inverse problem of finding reliable perception from fuzzy images can be solved through the process of "unconscious inference" (\textbf{Helmholtz Hypothesis}) \cite{GHatfield2002}.
This implies that human vision is incomplete, and that contextual details are inferred by the unconscious mind to create a super-resolution image. Let $p(\bm{x}|\bm{y};\gamma)$ be a probabilistic distribution of image $\bm{x}$ reconstructed from the under-sampled data $\bm{y}$. Here $\bm{x}$ can be obtained from $\bm{y}$ by solving the log-MAP estimation as
\begin{equation}\label{eq_feature_1}
\underset{\bm{x}}{\arg \max }\; \log p(\bm{x}|\bm{y};\gamma),
\end{equation}
where $\gamma$ is an environmental parameter in classical "unconscious inference" or inverse problem. If the prior probability can be modeled as a distribution $p(\bm{x})\propto\exp(-\mathcal{R}(\bm{x})/\gamma)$ of the geometric characteristics $\mathcal{R}(\bm{x})$, and the conditional probability can be exploited as $p(\bm{y}|\bm{x})\propto\exp(-\mathcal{H}(\bm{x};\bm{y}))$, then the classical CS-MRI reconstruction optimization problem can be formulated as follows \cite{Fan2021}
\begin{equation}
\min_{\bm{x}}\left\{\mathcal{E}(\bm{x}):= \mathcal{H}(\bm{x};\bm{y})+\gamma \mathcal{R}(\bm{x})\right\},
\label{eq1}
\end{equation}
where $\mathcal{H}(\bm{x};\bm{y})=\frac{1}{2}\left\|\mathcal{T} \bm{x}-\bm{y}\right\|_{2}^{2}$ is a data fidelity term, $\mathcal{R}(\bm{x})$ is a regularizer with either a known or unknown geometric prior, $\gamma$ is a regularization parameter that balances data fidelity term and regularizer.

\subsection{Proposed Learnable CS-MRI Unrolled framework}
The problem presented in (\ref{eq1}) can be effectively resolved using a two-step optimization procedure, which not only avoids time-consuming numerical optimization  but also reconstructs more spatial-frequency geometric priors. In details, we split (\ref{eq1}) into two subproblems as follows
\begin{itemize}
\item[(P1)] linear reconstruction subproblem:
\[\bm{m}_{\ell}\in \{\bm{x}: \mathcal{T} \bm{x}=\bm{y}\};\]
\item[(P2)] geometric prior reconstruction subproblem:
\[\bm{x}_{\ell}:=\underset{\bm{x}}{\arg \min } \left\{\frac{1}{2}\|\bm{x}-\bm{m}_{\ell}\|_{2}^{2}+\lambda_\ell\mathcal{R}(\bm{x})\right\}.\]
\end{itemize}
To learn more detailed information about many specific features that are invariant to contextual information and uninformative
intensity-variations, we consider a generative feature correction
module to solve the linear system
\[\mathcal{T} \bm{x}=\bm{y},\]
which involves a two-grid-cycle scheme of classical multi-grid algorithm \cite{McCormick1987,KChen05}.

One basic two-grid-cycle scheme with an initial guess $\bm{\bar{x}}^{h}$ and right-hand side $\bm{y}^{h}$ on fine grid $h$ can be denoted by
$$
\bm{x}^{h} \leftarrow \text{MG}^{h}\left(\bm{\bar{x}}^{h}, \bm{y}^{h},\mathcal{T}^h\right),
$$
which is given by the following steps:\\
\textbf{Step 1. (pre-relaxation)}: Perform $\nu_0$ sweeps to approximate the solution $\bm{x}^h_l$ of $\mathcal{T}^{h} \bm{x}^{h}=\bm{y}^{h}$ with initial guess $\bm{\bar{x}}^{h}$;\\
\textbf{Step 2. (coarse grid solution)}: Compute $\bm{r}^{h}=\bm{y}^{h}-\mathcal{T}^{h} \bm{x}^{h}_l$ and $\bm{r}^{2h}=I_{h}^{2h} \bm{r}^{h}$, and then seek an approximation $e^{2 h}$ to the solution of $\mathcal{T}^{2h}\bm{e}^{2h}=\bm{r}^{2h}$;\\
\textbf{Step 3. (coarse grid correction)}: Let $\bm{x}^{h}_m \leftarrow \bm{x}^{h}_l+I_{2h}^{h}\bm{e}^{2h}$;\\
\textbf{Step 4. (post-relaxation)}: Run $\nu_1$ sweeps to approximate the solution $\bm{x}^h_r$ of $\mathcal{T}^{h}\bm{x}^{h}=\bm{y}^{h}$ with initial guess $\bm{x}^{h}_m$ \\
where $\nu_0$ and $\nu_1$ are the numbers of relaxation sweeps to be done before and after nested process on fine grid $\Omega^{h}$, respectively. $l$ and $m$ denote the pre-relaxation step and coarse grid correction step, respectively. $h$ and $2h$ represent fine and coarse grids, respectively. $I_{h}^{2h}$ is a restriction operator which restricts residual from fine to coarse grid, and $I_{2h}^{h}$ is a prolongation operator which restricts the corrected residual from coarse to fine grid.

In the following we  bridge the correspondences between well established
principles in conventional modeling methods and CNN reconstruction
networks.

\subsubsection{Correction-distillation architecture}
Inspired by the excellent performance of multigrid in multi-scale error correction, we propose a learnable correction-distillation framework which is embedded to quickly extract low frequency increment in $k$-space (\emph{see} Fig.\ref{fig_architecture}(a)). At each stage $\ell$, we thus construct two modules based on \textbf{Steps 1 -- 3} to generate
solution iterations for  subproblem (P1) as follows
\begin{align}
\bm{h}_{\ell}&=\bm{x}_{\ell-1}-\eta_{\ell}{\mathcal{T}^{H}}(\mathcal{T}\bm{x}_{\ell-1}-\bm{y}),
\label{eq4}\\
\bm{e}^{2h}_\ell&={{F}^{H}}\mathcal{S}_{\ell}\mathcal{I}_{\ell}^{R}(\bm{y}-\mathcal{T} \bm{h}_{\ell}),
\label{eq5}\\
\bm{m}_{\ell}&=\text{Concat}(\bm{h}_{\ell},\mathcal{I}_{\ell}^{P}\bm{e}^{2h}_\ell,\mathcal{M}_\ell^{\beta_\ell}),
\label{eq6
\end{align}
where $\bm{x}_0=\mathcal{T}^{H}\bm{y}$, $\eta_{\ell}$ is a learnable step-length, { ${F}^{H}$ is an inverse discrete Fourier transform.} $\text{Concat}(\cdot,\cdot,\cdot)$ is a correction operation, and $\mathcal{M}_\ell^{\beta_\ell}$ is filled with noise level $\beta_\ell$ learned from condition module $\mathcal{CM}$.

More importantly, we observed that correction step (\ref{eq5}) with error-based solution operator $\mathcal{S}_\ell$ (\emph{see} Fig.\ref{fig_architecture}(e)) can effectively learn significant details to yield improved solution of subproblem (P1). It turns out  that the results have almost no significantly differences   with or without using the  post-relaxation \textbf{Step 4}. To this end, we set the unifying CS-MRI framework without post-relaxation as default in this work (\emph{see} Fig.\ref{fig_architecture}(a)).

However, if we incorporate the learnable geometric prior distillation stage (P2) into the solution fidelity of image inverse problem (\ref{eq0}),
  the network can learn the expected MR image priors. For this purpose, we add the geometric prior fidelity module (\emph{see} Fig.\ref{fig_architecture}(d)) in our network, which is designed by solving the equation
\begin{align}
\bm{x}_{\ell}-\mathcal{N}(\bm{x}_{\ell})=\bm{m}_{\ell},
\label{eq7}
\end{align}
where the nonlinear function $\mathcal{N}(\bm{x})=-\lambda\frac{\partial\mathcal{R}(\bm{x})}{\partial\bm{x}}$ denotes the geometric characteristics of $\bm{x}$.

The above correction-distillation architecture (\ref{eq4}) - (\ref{eq7}) can be performed by a pre-relaxation module $\mathcal{D}_\ell$, a correction module $\mathcal{C}_\ell$, and a geometric prior distillation module $\mathcal{P}_\ell$. It's well known that the module $\mathcal{D}_\ell$ often results in heavy artifacts, while the correction module $\mathcal{C}_\ell$ is used to quickly generate the low frequency correction in $k$-space. Then the explainable geometric prior distillation module $\mathcal{P}_\ell$ is learned to produce more texture priors.

\subsubsection{Pre-relaxation module $\mathcal{D}_\ell$} At each stage $\ell$,
to provide the basis-solution guarantee in high-frequency layer, we can obtain an approximate solution $\bm{h}_{\ell}$ of subproblem (P1)
by solving the least-square problem
\begin{equation*}\label{least-square_prob}
\min\limits_{\bm{x}} \left\{\mathcal{H}(\bm{x};\bm{y})=\frac{1}{2}\left\|\mathcal{T} \bm{x}-\bm{y}\right\|_{2}^{2}\right\}.
\end{equation*}
Especially, the gradient descent \eqref{eq4} with the parameter $\eta_\ell$ is implemented to perform
a linear reconstruction, which is also called as the pre-relaxation module $\mathcal{D}_\ell$ defined by
\begin{align}
\begin{split}
\bm{h}_{\ell}=\mathcal{D}_\ell(\bm{x}_{\ell-1},\eta_{\ell},\bm{y},\mathcal{T})=\bm{x}_{\ell-1}-\eta_{\ell}{ \mathcal{T}^{H}}(\mathcal{T} \bm{x}_{\ell-1}-\bm{y}).
\end{split}
\label{eq8}
\end{align}

 The module $\mathcal{D}_\ell$ corresponding to Eq.(\ref{eq8}) is directly used to generate the preliminary approximation from the previous $\bm{x}_{\ell-1}$ (\emph{see} Fig.\ref{fig_architecture}(b)). It is well known that the step-length $\eta_\ell$ should be positive, and decrease smoothly as the increase of iterations in traditional model-based methods. To enhance network flexibility, we set the step-length $\eta_\ell$ to be
learnable during iterations. There are a variety of ways to use training data to learn step-length $\eta_\ell$. Here, we adopt the learned step-length $\eta_{\ell}$ from condition module $\mathcal{CM}$ in \S \ref{section_CM}.

\subsubsection{Correction module $\mathcal{C}_\ell$} The proposed correction module corresponds strictly to a composition of coarse grid solution \eqref{eq5} and coarse grid correction \eqref{eq6} in two-grid-cycle update scheme. The implementation process of solution correction $\mathcal{C}_\ell$ is shown in Fig.\ref{fig_architecture}(c).

Similar to {\bf Step 2} in classical two-grid-cycle scheme, we first design the reconstruction residual block in $k$-space as
\[\bm{r}_\ell^{h}=\bm{y}^h-\mathcal{T}^h\bm{h}^h_{\ell}.
\]
Then we restrict the reconstructed $k$-space residual $\bm{r}_\ell^{h}$ on fine level to $\bm{r}_\ell^{2h}$ on the coarse level, which is denoted by
\[\bm{r}_\ell^{2h}=\mathcal{I}_\ell^{R}\bm{r}_{\ell}^{h}=\mathcal{I}_\ell^{R}(\bm{y}^h-\mathcal{T}^h\bm{h}^h_{\ell}):=\mathcal{I}_\ell^{R}(\bm{y}-\mathcal{T}\bm{h}_{\ell}),\]
where a convolution $\mathcal{I}_\ell^{R}$ with stride 2 pixel { and $2\times2$ kernel} is learnt to exploit the restriction operator $I_{h}^{2 h}$. Note that $\mathcal{I}_\ell^{R}$ and $I_{h}^{2 h}$ have the similar meaning but the former is for learning implementation while the latter is for mathematical description.

To compute a more accurate correction $\bm{e}^{2 h}_\ell$ as the solution of $\mathcal{T}^{2h}_\ell \mathbf{e}^{2h}_\ell=\bm{r}^{2h}_\ell$ and reduce
the computational costs, we propose
to learn an error-based solution operator $\mathcal{S}_{\ell}$ to quickly generate low frequency correction $\mathcal{S}_{\ell}\bm{r}_\ell^{2h}$ in $k$-space, and then use ${ {F}^{H}}$ to convert corrected error from $k$-space to image space. Mathematically, we can formulate it as
\[\bm{e}^{2h}_\ell={ {F}^{H}}\mathcal{S}_{\ell}\bm{r}_\ell^{2h},\]
where the operation $\mathcal{S}_\ell$ is composed of four residual blocks (ResBlocks) \cite{He2016} with $p$ channels, two $3\times 3$ convolution layers with $p$ channels and a skip connection as shown in Fig.\ref{fig_architecture}(e). Obviously, the block ${ {F}^{H}}\mathcal{S}_{\ell}$ can be seen as an approximation to inverse operation $(\mathcal{T}^{2 h}_\ell)^{-1}$.
We also remark that Batch Normalization (BN) is not adopted because some recent papers showed that BN layer is more likely to yield undesirable representations when the network becomes deeper and more complex
\cite{Wang2018,Zhang2018a}.

Next, we interpolate the correction $\bm{e}^{2h}_\ell$ on coarse level to the fine level by
\[\bm{e}_{\ell}^{h}=\mathcal{I}_\ell^{P}\bm{e}^{2h}_\ell,\]
where $\mathcal{I}_\ell^{P}$ is implemented by a learnable transpose convolution with stride 2 pixel to represent the prolongation operator $I_{2 h}^{h}$.

Using channel concatenation
to add contextual information $\bm{e}^{h}_\ell$ has been widely used
in traditional CNN architectures.
According to {\bf Step 3}, we thus apply
a concatenate operation $\text{Concat}(\cdot,\cdot,\cdot)$ to produce the refined
correction module $\mathcal{C}_\ell$ as
\begin{equation*}
\begin{split}
\bm{m}_{\ell}&=\mathcal{C}_\ell(\bm{h}_{\ell},\bm{y},\mathcal{T},\mathcal{I}_{\ell}^{R},\mathcal{S}_{\ell},\mathcal{I}_{\ell}^{P},\beta_{\ell})\\
&=\text{Concat}(\bm{h}_{\ell},\bm{e}_{\ell}^{h},\mathcal{M}_\ell^{\beta_\ell})=\text{Concat}(\bm{h}_{\ell},\mathcal{I}_{\ell}^{P}\bm{e}^{2h}_\ell,\mathcal{M}_\ell^{\beta_\ell}),
\end{split}
\label{eq9}
\end{equation*}
where the additional noise level map $\mathcal{M}_\ell^{\beta_\ell}$, which is filled with the output $\beta_\ell$ of condition module $\mathcal{CM}$ \cite{Zhang2020,You2021}, has the same size as $\bm{h}_{\ell}$ and can make network more flexible to train multi-sampling-ratio CS-MRI task.

\subsubsection{Geometric prior distillation module $\mathcal{P}_\ell$} We employ a geometric prior distillation module $\mathcal{P}_\ell$ \cite{Fan2021} to refine the compromised image structure and geometric features (\emph{see} Fig.\ref{fig_architecture}(d)).

Unfortunately, it is difficult to directly obtain the close-form solution $\bm{x}_{\ell}=(I-\mathcal{N})^{-1}(\bm{m}_{\ell})$ of the optimal condition \eqref{eq7} of the geometric prior
distillation stage (P2). Naturally, if the operator $\mathcal{N}$ satisfies the contraction condition $\|\mathcal{N}\|<1$, $(I-\mathcal{N})^{-1}$ can be simplified by using series expansion as follows
\begin{equation}
\begin{split}
\bm{x}_{\ell}&=(I-\mathcal{N})^{-1}(\bm{m}_{\ell})=\left(\sum_{i=0}^{k}\mathcal{N}^{i}+R(\mathcal{N}^{k})\right)(\bm{m}_{\ell})\\
&\in\text{span}(\bm{m}_{\ell},\mathcal{N}(\bm{m}_{\ell}),...\,\mathcal{N}^{k}(\bm{m}_{\ell}), R(\mathcal{N}^{k})),\\
\end{split}
\label{eq11}
\end{equation}
where $R(\mathcal{N}^{k})$ is the remainder, and can be approximated by $\mathcal{N}^{k+1}(\bm{m}_{\ell})$.

A more flexible representation of the non-linear operation $\mathcal{N}^j$ in (\ref{eq11}) could be approximated by CNN block $\mathcal{N}^j_{\mathcal{K},\ell}$ with many embedded convolution blocks and ReLU layers, where $\mathcal{G}_{\ell,k}$ is designed to replace $\mathcal{N}$ for learning multi-scale geometric priors of $\bm{x}_{\ell}$.
Finally, a $1\times1$ convolution $\mathcal{W}_{\ell}=\left\{\mathcal{W}_{\ell,j}\right\}_{j=0}^{k}$
is used to fuse features $\mathcal{F}_{\ell,k}$, i.e.
\begin{equation}
\begin{split}
\bm{x}_{\ell}&=\mathcal{P}_\ell(\bm{m}_\ell)=\mathcal{K}_{\ell,0}\ast\bm{m}_{\ell}+\sum_{j=1}^{k} \mathcal{N}^j_{\mathcal{K},\ell}\left(\bm{m}_{\ell}\right)\end{split}\label{eq11A}\end{equation}
\[\begin{split}
&=\mathcal{K}_{\ell,0}\ast\bm{m}_{\ell}+\sum_{j=1}^{k}\left(\mathcal{K}_{\ell,j}\ast\mathcal{B}^{j-1}_{\ell}(\mathcal{A}_{\ell}(\bm{m}_{\ell}))\right)\\
&=\mathcal{W}_{\ell}\ast\mathcal{F}_{\ell,k}=\mathcal{W}_{\ell}\ast\text{Concat}\left(\mathcal{F}_{\ell,k-1},\mathcal{G}_{\ell,k}*\mathcal{F}_{\ell,k-1})\right),
\end{split}\]
where $\mathcal{A}_{\ell}$ and $\mathcal{B}_{\ell}$ are two learnable blocks with the $3\times3$ convolution $\bm{\kappa}$ and ReLU such that
$\mathcal{A}_{\ell}(\cdot):=\text{ReLU}(\bm{\kappa}(\cdot)):\mathbb{R}^{m\times n}\rightarrow \mathbb{R}^{m\times n\times p}$, $\mathcal{B}_{\ell}(\cdot):=\text{ReLU}(\bm{\kappa}(\cdot)):\mathbb{R}^{m\times n\times p}\rightarrow \mathbb{R}^{m\times n\times p}$, $\mathcal{B}_{\ell}^{k-1}=\underbrace{\mathcal{B}_{\ell}\circ\dots\circ\mathcal{B}_{\ell}}_{k-1}$, $\mathcal{G}_{\ell,k}=(\underbrace{\mathcal{A}_{\ell},\mathcal{B}_{\ell},\dots,\mathcal{B}_{\ell}}_{k})$, each $1 \times 1$ convolution $\mathcal{K}_{\ell,j}$ in  $\mathcal{K}_{\ell}=\left(\mathcal{K}_{\ell,0}, \mathcal{K}_{\ell,1}, \ldots, \mathcal{K}_{\ell,k}\right) $ is a linear combination of $\left\{\mathcal{W}_{\ell,j}\right\}_{j=0}^{k}$.  We  also refer the interested readers to \cite{Fan2021} for more details.

\subsubsection{Condition module $\mathcal{CM}$}\label{section_CM} The pre-relaxation reconstruction step (\ref{eq4}) involves choosing/tuning hyper-parameter $\eta_\ell$
that significantly affects model performance.
A suitable step-length $\eta_\ell$ can achieve a satisfactory reconstruction as long as the CS $k$-space sampling-ratio is fixed.
However, it is generally difficult to determine a suitable $\eta_\ell$ without which the reconstructed image can be unsatisfactory,
especially when switching to other medical data with multi-sampling-ratio or other sampling masks. In general,  manual  tuning to determine optimal task-specific hyper-parameters
require high computational costs and often lead  to sub-optimal results.

To make multi-sampling-ratio tasks more flexible in our unifying model and tackle the inefficiency of hyper-parameter tuning (which
includes both adaptive model parameters and hyper-parameters), we learn the   step-length $\eta_{\ell}$ and noise level $\beta_\ell$ in our condition module $\mathcal{CM}$. We remark that CS sampling-ratio $\alpha$ influences the degree of ill-posedness, so it is suitable as the input of condition module $\mathcal{CM}$, resulting in the commonly used parameter prediction module \cite{You2021}, i.e.,
\begin{equation}
[\bm{\eta},\bm{\beta}]=\mathcal{CM}(\alpha),
\label{eq12}
\end{equation}
where $\bm{\eta}=[\eta_1,\eta_2,...,\eta_{n_s}]$ and $\bm{\beta}=[\beta_1,\beta_2,\dots,\beta_{ n_s}]$. $\mathcal{CM}$ is composed of three hidden fully connected layers with $p$ neurons,
ReLU is employed as the first two activation functions and softplus \cite{Glorot2011} is used as the last activation function.

\subsubsection{Loss function}
For the training data pairs $\left\{\left(\bm{y}^{i}, \bm{x}^{i}\right)\right\}_{i=1}^{n_t}$ (with $n_t$   the total number of samples), under-sampled $k$-space data $\bm{y}$ and  initialization $\bm{x}_{0}$ are used as inputs, the final output $\bm{x}_{n_s}$ is obtained by the proposed framework. The loss function is commonly employed to seek the real target image $\bm{x}^*$ by minimizing the distance measure between $\bm{x}_{\scriptscriptstyle n_s}$ and $\bm{x}$. Here we adopt $\ell_1$-loss rather than $\ell_2$-loss which is insufficient to capture perceptually relevant components (e.g., high-frequency geometric details) \cite{KJiang2018} to enlarge the original loss.

Inspired by deep supervision technique \cite{Lee2015}, we add an auxiliary constraint on a network branch to supervise the trunk network   to make the optimization process more flexible. The total loss is defined as follows
\begin{equation}
\mathcal{L}=\frac{1}{n_t mn} \sum_{i=1}^{n_t}\left(\left\|\bm{x}^{i}_{\scriptscriptstyle\overline{ n_s}}-\bm{x}^{i}\right\|_{1}+\left\|\bm{x}^i_{\scriptscriptstyle n_s}-\bm{x}^{i}\right\|_{1}\right),
\label{eq13}
\end{equation}
where $\overline{ n_s}=\frac{ n_s+1}{2}$, $mn$ is the size of $\bm{x}^{i}$ from \eqref{eq11A}, and $ n_s$ is the total stage number of the proposed framework. $\bm{x}^{i}_{\scriptscriptstyle\overline{ n_s}}$ and $\bm{x}^i_{\scriptscriptstyle n_s}$ are two outputs of the sample $\bm{x}^{i}$ on stages $\overline{ n_s}$ and $n_s$, respectively.


\subsubsection{Parameters and Initialization}
The three modules in   stage $\ell$ of the proposed framework strictly inplement our correction-distillation updates (\eqref{eq4} to \eqref{eq7}). The learnable parameter set $\Theta=\left\{\eta_{\ell},\beta_{\ell}, \mathcal{A}_{\ell},\mathcal{B}_{\ell},\mathcal{K}_{\ell},\mathcal{I}_{\ell}^{R},\mathcal{S}_{\ell},\mathcal{I}_{\ell}^{P}\right\}_{\ell=1}^{ n_s}$ consists of step-length $\eta_{\ell}$, noise { level} $\beta_{\ell}$, CNN blocks $\mathcal{A}_{\ell}$ and $\mathcal{B}_{\ell}$ with $3\times3$ convolution and ReLU, $1 \times 1$ convolution $\mathcal{K}_{\ell}$, restriction operator $\mathcal{I}_{\ell}^{R}$, solution operator $\mathcal{S}_{\ell}$ and prolongation operator $\mathcal{I}_{\ell}^{P}$. All these parameters are learned as neural network parameters by minimizing the loss (\ref{eq13}).

Similar to traditional model-based method, the proposed framework also requires an initial input $\bm{x}_{0}={ \mathcal{T}^{H}}\bm{y}$.
The convolution network is initialized with Kaiming Initialization \cite{He2015}. The model parameters $\{p, k, n_{s}\}$ are initialized as $\{32, 8, 13\}$ respectively.

\section{Experiments and results}
In this section, we verify the advantages of the proposed CGPD-CSNet in comparing with the state-of-the-art methods through various experiments. Peak Signal to Noise Ratio (PSNR) and Structural Similarity Index Measure (SSIM) are employed to evaluate their performances, the values (M-PSNR and M-SSIM)
are all obtained by averaging on tested data.

\begin{table*}
\centering
\caption{Quantitative assessment with M-PSNR values of different stage number $ n_s$ and different geometric distillation depths $k$ using Cartesian sampling brain dataset with different CS ratios.}
\label{table_stage_number}
\setlength{\tabcolsep}{2.4mm}{
\begin{tabular*}{\hsize}{@{}@{\extracolsep{\fill}}lcccccc ccccccc@{}}     
\toprule[1.5pt]
        \multirow{2}*{\textbf{CS Ratio}}&\multicolumn{6}{c}{\textbf{Stage number $ n_s$}}   &\multirow{2}*{ }&\multicolumn{6}{c}{\textbf{Geometric distillation depths $k
        }}\\
        \cmidrule[0.7pt](lr){2-7} \cmidrule[0.7pt](lr){9-14}
          & 5 & 7 & 9 & 11 & \textbf{13} & 15  &   & 4 & 5 & 6 & 7 & \textbf{8} & 9\\
          \midrule[0.8pt]
        10\% & 32.54 & 33.59 & 33.80 & 34.16 & \textbf{34.23} & 34.25    &  & 33.69 & 33.91 & 34.02 & 34.22 & \textbf{34.23} & 34.37\\
        20\% & 37.43 & 38.31 & 38.61 & 38.88 & \textbf{39.02} & 39.00     &  & 38.46 & 38.65 & 38.75 & 38.93 & \textbf{39.02} & 39.02\\
        30\% & 39.56 & 40.46 & 40.76 & 41.01 & \textbf{41.16} & 41.15     &  & 40.64 & 40.83 & 40.89 & 41.10 & \textbf{41.16} & 41.14\\
        40\% & 43.45 & 44.38 & 44.71 & 44.92 & \textbf{45.07} & 45.08    &   & 44.71 & 44.84 & 44.80 & 44.99 & \textbf{45.07} & 45.01\\
        50\% & 45.99 & 46.72 & 47.05 & 47.19 & \textbf{47.31} & 47.38    &   & 47.05 & 47.13 & 47.10 & 47.28 & \textbf{47.31} & 47.29\\
       \bottomrule[1.5pt]
    \end{tabular*}}
\end{table*}

\subsection{Implementation Details}
We evaluate the performance of the proposed and comparison methods on the widely used { brain MR dataset \cite{NNA2015}} using three classical types of sampling masks \cite{Yang2020}
for training within 800 MR brain medical images and testing within $50$ images as ISTA-Net+ \cite{Zhang2018}. These MR brain images are T1-weighted 2D images from health and Alzheimer's disease patients from different devices.

Besides, we employ a cardiac dataset \cite{Zheng2019}, which is established based on \cite{Andreopoulos2008} and contains 4480 cardiac real-valued MR images from 33 patients, to further evaluate the performance of comparison methods. The first 30 patients' samples of 4180 MR images are set as training set while the last 3 patients' of 300 MR images are set as testing set. The training data is augmented by conducting horizontal and vertical flipping.

We employ Pytorch to implement the proposed framework. We use Adam optimization \cite{Kingma2014} with a learning rate of $0.0001$ and batch size $1$ to train network for 500 epochs to ensure convergence. In our unifying multi-sampling-ratio CS-MRI framework, the mixed CS sampling ratios are set as $\{10\%, 20\%, 30\%, 40\%, 50\%\}$ for Cartesian sampling mask and pseudo radial sampling mask, and are set as $\{5\%, 10\%, 20\%, 30\%, 40\%\}$ for 2D random sampling mask. All experiments are performed on a workstation with Intel Xeon CPU E5-2630 and Nvidia Tesla V100 GPU.

\subsection{Intra-Method Evaluation}
We first conduct four groups of experiments to investigate the role of different network components in the proposed CGPD-CSNet on reconstruction performance, including stage number $ n_s$, geometric prior distillation depths $k$, as well as ablation study about the proposed correction module $\mathcal{C}_\ell$ and condition module $\mathcal{CM}$, different shared settings of CGPD-CSNet.

\subsubsection{Test of stage number $ n_s$}
To evaluate the effectiveness of different stage number on reconstruction performance, we tune the stage number $ n_s$ from 5 to 15 at 2 intervals. Using the CNN architectures with different stage number $ n_s$ to Cartesian sampling brain dataset with different CS sampling ratios, the M-PSNR values of the reconstructed results are summarized in Table \ref{table_stage_number}. We can observe that the reconstruction performance gradually improves with the increase of stage number $ n_s$ while the M-PSNR value becomes stable after $ n_s\geq13$. Based on this observation, the 13-stage configuration is a preferable setting to balance the reconstruction performance and computational costs, and we fix $ n_s=13$ throughout all the experiments.

\subsubsection{Test  of geometric distillation depths $k$}
To explore the relationship between different geometric prior distillation depths and reconstruction performance, we tune the geometric distillation depths $k$ from 4 to 9. The M-PSNR values of the reconstructed images on brain dataset using different geometric distillation depths and  different Cartesian CS sampling ratios are summarized in Table \ref{table_stage_number}. The reconstruction performance improves slowly after $k\geq8$. Considering the tradeoff between network complexity and reconstruction performance, we set geometric distillation depths $k=8$ in all configurations.

\subsubsection{Ablation study}
Next, we conduct a group of ablation studies to better evaluate the effectiveness of correction module $\mathcal{C}_\ell$ and condition module $\mathcal{CM}$ on the CS-MRI reconstruction performance. The comparisons are shown in Table \ref{table_ablation}.

\begin{table}
\centering
\caption{Quantitative assessment with M-PSNR values of different combinations of correction module $\mathcal{C}_\ell$ and condition module $\mathcal{CM}$ using Cartesian sampling mask on brain dataset. The best and second places are highlighted in Bold font and underline ones, respectively.}
\label{table_ablation}
\setlength{\tabcolsep}{2mm}{
\begin{tabular}{ccccccccc}
\toprule[1.5pt]
        \multirow{2}*{\textbf{Variant}} &\multirow{2}*{\textbf{$\mathcal{C}_\ell$}} &\multicolumn{2}{c}{\textbf{$\mathcal{CM}$}} & \multicolumn{5}{c}{\textbf{CS Ratio}}\\ \cmidrule[0.7pt](lr){3-9}
         & & $\eta_\ell$& $\beta_\ell$ &10\%&20\%& 30\% & 40\% & 50\% \\ \midrule[0.8pt]
        (a)  & -&-&-& 32.53&38.32& 40.63&44.60&46.66 \\
        (b)  & -&+&+&32.61&38.40&40.69&44.80&47.16\\
        (c)  & +&-&-&33.54 &38.54&40.76&44.66&46.89\\
        (d)  & +&-&+&33.90 &38.63& 40.81&44.78&47.11\\
        (e)  & +&+&-&\underline{34.00} &\underline{39.00}&\underline{41.13}&\underline{45.04}&\underline{47.30}\\
        (f)  & +&+&+&\textbf{34.23}&\textbf{39.02}&\textbf{41.16}&\textbf{45.07}& \textbf{47.31} \\
        \bottomrule[1.5pt]
    \end{tabular}}
\end{table}

\begin{table}
\centering
\caption{Quantitative assessment with M-PSNR values with different shared settings of proposed CGPD-CSNet using Cartesian sampling mask on brain dataset. Bold font
in the table is the best of the eight variants.}
\label{table_shared}
\setlength{\tabcolsep}{1.8mm}{
\begin{tabular}{ccccccc}
\toprule[1.5pt]
        \multirow{2}*{\textbf{Variant}} &\multirow{2}*{\textbf{Shared setting}} & \multicolumn{5}{c}{\textbf{CS Ratio}}\\ \cmidrule[0.7pt](lr){3-7}
        &&10\%&20\%& 30\% & 40\% & 50\% \\ \midrule[0.8pt]
        (a)  &Shared $\mathcal{C}_\ell$, $\mathcal{P}_\ell$ & 32.84&38.31& 40.55&44.65&47.02 \\
         (b)  &Shared $\mathcal{P}_\ell$& 33.55&38.55&40.76&44.79&47.09\\
        (c)  & Shared $\mathcal{C}_\ell$ &33.56 &38.71&40.93&44.91&47.17\\
        (d)  &Unshared (default) &\textbf{34.23}&\textbf{39.02}&\textbf{41.16}&\textbf{45.07}& \textbf{47.31}\\
        \bottomrule[1.5pt]
    \end{tabular}}
\end{table}


\begin{table*}
\centering
\newsavebox{\mybox}
\begin{lrbox}{\mybox}
\setlength{\tabcolsep}{2.4mm}{
\begin{tabular*}{\hsize}{@{}@{\extracolsep{\fill}}llcccccc@{}}      
\toprule[1.5pt]
\multirow{2}*{\textbf{Dataset}} &\multirow{2}*{\textbf{Method}} & \multicolumn{5}{c}{\textbf{CS Ratio}} &  \textbf{Time (s)} \\ \cmidrule[0.7pt](lr){3-7}

& & 10\% & 20\% & 30\% & 40\% & 50\% & \textbf{CPU/GPU} \\
\midrule[0.8pt]
\multirow{11}*{\textbf{Brain}} &Zero-filling\cite{Bernstein2001} & 23.85/0.5743& 26.40/0.6611&   28.64/0.7249&   31.12/0.7863 &  32.93/0.8230 & 0.0019/---\\
&TV\cite{Lustig2007} & 25.56/0.6685& 29.92/0.7939& 32.88/0.8582&    37.19/0.9285 & 39.91/0.9472 & 1.2330/---     \\
&{SIDWT} & 25.19/0.6782& 28.96/0.7782& 31.87/0.8476& 36.50/0.9234   & 39.70/0.9493 & 18.033/---\\
&{PBDW}\cite{Qu2012}  & 26.82/0.7355& 32.42/0.8661& 35.04/0.9116& 39.76/0.9599  & 42.52/0.9727   & 65.1656/---\\
&{PANO}\cite{Qu2014}  &28.98/0.7897&34.64/0.8934& 36.79/0.9244& 41.54/0.9661 & 44.14/0.9733& 38.2577/---\\
&{DC-CNN}\cite{Schlemper2018}  & 30.00/0.8232&35.73/0.9204& 38.48/0.9493&42.87/0.9765& 45.53/0.9842& ---/0.0161 \\
&{ISTA-Net}\cite{Zhang2018} & 30.28/0.8256 &36.59/0.9334 &39.14/0.9566&43.93/0.9820&  46.77/0.9885&---/0.0189\\
&{ISTA-Net+}\cite{Zhang2018}& 30.86/\underline{0.8455}&37.05/0.9356 &39.62/0.9597&44.24/0.9826& 46.78/0.9882 &{---/0.0200}\\
&{FISTA-Net}\cite{Xiang2020} & \underline{31.06}/0.8444&    \underline{37.27}/\underline{0.9395} &  \underline{39.86}/\underline{0.9607} &   \underline{44.50}/\underline{0.9835} & \underline{46.99}/\underline{0.9889} &{---/0.0297}\\
&{CGPD-CSNet} &\textbf{34.23}/\textbf{0.9126}&\textbf{39.02}/\textbf{0.9553}&\textbf{41.16}/\textbf{0.9685}&\textbf{45.07}/\textbf{0.9852}& \textbf{47.31}/\textbf{0.9897} &{---/0.0721}\\
\midrule[0.8pt]
\multirow{11}*{\textbf{Cardiac}} &Zero-filling\cite{Bernstein2001} & 20.29/0.5286&   22.99/0.6059&   25.50/0.6698&   28.05/0.7371&   29.97/0.7736 & 0.0017/---\\
&TV\cite{Lustig2007} & 22.26/0.6139& 27.94/0.7786& 31.32/0.8535&    36.27/0.9274 & 38.87/0.9486 & 1.2446/---     \\
&{SIDWT} & 21.56/0.6233& 26.23/0.7480& 29.62/0.8279& 35.02/0.9102   & 37.97/0.9416 & 21.1174/---\\
&{PBDW}\cite{Qu2012}  & 23.73/0.7123& 30.46/0.8666& 33.18/0.9122& 38.58/0.9642  & 41.20/0.9761   & 69.8070/---\\
&{PANO}\cite{Qu2014}  &25.94/0.7606&32.45/0.8809& 34.61/0.9155& 40.81/0.9697 & 43.07/0.9730& 39.0812/---\\
&{DC-CNN}\cite{Schlemper2018}  & 26.31/0.7863&33.51/0.9153& 36.29/0.9484&42.56/0.9836& 45.25/0.9888& ---/0.0131 \\
&{ISTA-Net}\cite{Zhang2018} & 27.81/0.8185&34.23/0.9299&36.70/0.9557&43.94/0.9892&  47.51/0.9941&---/0.0134\\
&{ISTA-Net+}\cite{Zhang2018}& \underline{27.84}/0.8236&34.57/0.9339&37.20/0.9589&44.56/0.9903&  48.07/0.9946 &{---/0.0145}\\
&{FISTA-Net}\cite{Xiang2020} & 27.74/\underline{0.8280}&     \underline{34.70}/\underline{0.9347}&  \underline{37.41}/\underline{0.9600} &  \underline{44.83}/\underline{0.9905} & \underline{48.37}/\underline{0.9949} &{---/0.0229}\\
&{CGPD-CSNet} &\textbf{28.05/0.8516}&\textbf{35.32}/\textbf{0.9433}&\textbf{37.99}/\textbf{0.9649}&\textbf{45.47}/\textbf{0.9919}& \textbf{48.66}/\textbf{0.9952} &{---/0.0539}\\
\bottomrule[1.5pt]
\end{tabular*}}
\end{lrbox}
\caption{Quantitative assessment with M-PSNR and M-SSIM values using Cartesian sampling mask with different CS ratios on brain and cardiac datasets. The best and second best results are highlighted in Bold font and underlined ones, respectively.}
{\usebox{\mybox}}\label{table_comparison}
\end{table*}

Using condition module $\mathcal{CM}$ with two learnable parameters, our method enjoys the flexibility of handling CS-MRI problems with different sampling ratios through a single model. Contrast to the single variant (a) trained with the same Cartesian sampling masks without condition module $\mathcal{CM}$, variant (b) with the proposed condition module $\mathcal{CM}$ consistently outperforms variant (a) across all five CS ratios and brings average 0.18 dB improvement. When the proposed correction module $\mathcal{C}_\ell$ is embedded in reconstruction model, variant (f) with condition module $\mathcal{CM}$ can further improve reconstruction performance than variant (c) without condition module $\mathcal{CM}$ and increases average $0.48$ dB. Furthermore, comparison between variants (a) and (c) can also demonstrate the advantage of the proposed correction module $\mathcal{C}_\ell$. Especially in low sampling rate $10\%$, reconstruction performance improves $1.01$ dB.

From variants (c) and (d), it is clear to descript that noise level $\beta_\ell$ learned from condition module $\mathcal{CM}$ can improve the M-PSNR value by $0.17$ dB. We can also observe from variants (c) and (e) that step-length $\eta_\ell$ learned from condition module $\mathcal{CM}$ can greatly boost the reconstruction performance across all five CS ratios and achieve average $0.42$ dB improvement. Compared with variant (a), our variant (f) combining with correction module $\mathcal{C}_\ell$ and condition module $\mathcal{CM}$ can achieve average $0.81$ dB improvement across all five sampling ratios and $1.7$ dB in low CS ratio $10\%$ especially. These comparisons adequately verify the effectiveness of the proposed correction module $\mathcal{C}_\ell$ and condition module $\mathcal{CM}$.



\subsubsection{Module sharing configurations}
To demonstrate the flexibility of the proposed framework that does not have to be the same network parameter configurations in different stages, we conduct several variants of CGPD-CSNet that have different shared settings among stages. Table \ref{table_shared} lists the M-PSNR values for different shared architectures in case of using Cartesian sampling mask on brain dataset. Note that the best M-PSNR scores are achieved when using the default unshared version (d), which is the most flexible with largest number of parameters. The variant (a) that shares both $\mathcal{C}_\ell$ and $\mathcal{P}_\ell$ in all stages is least flexible with smallest number of parameters and achieves the worst performance. Especially if only $\mathcal{P}_\ell$ or $\mathcal{C}_\ell$ is shared, the variants (b) and (c) increase average 0.27 dB and 0.38 dB over the variant (a). So we adopt the default unshared version (d) to perform the following experiments.




Actually, we attribute the superiority of our method to the following three factors. Firstly, our method has a correction module $\mathcal{C}_\ell$ which can compensate the low frequency of reconstruction error in $k$-space. Secondly, the additional noise level $\beta_\ell$ learned from condition module $\mathcal{CM}$ adds the anti-interference, and the geometric prior distillation module $\mathcal{P}_\ell$ distills the lost contextual details. Thirdly, the learnable step-length $\eta_\ell$ can be adapted and optimized according to different CS sampling ratios, such that multi-sampling-ratio CS-MRI can be jointly learnt through a single model.



\subsection{Comparison to the state-of-the-art Methods}
Since our proposed CGPD-CSNet framework  merges the advantages of traditional methods and modern deep learning methods,
it is necessary to compare with related methods from both categories on performances and hence  generalizability.

We compare our new framework with the state-of-the-art methods including traditional methods (Zero-filling \cite{Bernstein2001}, TV\cite{Lustig2007}, { SIDWT}\footnote{https://github.com/ricedsp/rwt}), patch-based methods (patch-based directional wavelets (PBDW) \cite{Qu2012}, patch-based nonlocal operator (PANO) \cite{Qu2014}),
deep unfolding methods (DC-CNN \cite{Schlemper2018}, ISTA-Net \cite{Zhang2018}, ISTA-Net+ \cite{Zhang2018}, FISTA-Net \cite{Xiang2020}. { Following \cite{Zhang2018,Fan2021}}, the stage number of ISTA-Net, ISTA-Net+ and FISTA-Net is configured as 11. { Following \cite{Schlemper2018},} the stage number of DC-CNN is configured as 5. Cartesian sampling mask which is widely used in clinic \cite{Lustig2007,Qu2012} is employed to acquire under-sampling $k$-space data. Apart from other comparison methods that have to train different model according to every CS sampling ratio, we cast the training process for five different CS sampling ratios as a unified joint learning throughout all the experiments.

\begin{figure*}[htbp]
\centering
\includegraphics[width=0.98\textwidth]{fig/figs/TCGD_TNNLS_Reconstruction_fig03_v1.pdf}
\caption{{ Qualitative comparisons of our method with other popular methods using three sampling masks on brain and cardiac datasets.}}
\label{fig_results}
\end{figure*}

\begin{figure*
\centering
\includegraphics[width=0.98\textwidth]{fig/figs/TCGD_TNNLS_Reconstruction_fig04_v1.pdf}
\caption{Violin visualizations of PSNR scores for different methods using Cartesian sampling mask with CS ratio 10\% on brain dataset (top row) and CS ratio 30\% on cardiac dataset (bottom row).}
\label{fig_violinplot}
\end{figure*}

\subsubsection{Quantitative evaluation}
The M-PSNR/M-SSIM scores of the proposed and comparison methods using Cartesian sampling mask with different CS ratios on brain and cardiac datasets are listed in Table \ref{table_comparison}.
 As a result, a significant increase of all M-PSNR/M-SSIM scores compared to the other methods across all CS ratios can be observed. In particular, our method achieves 3.18 dB improvement at low CS ratio 10\% and average 1.42 dB improvement over other comparison methods on brain dataset. And our method achieves average 0.47 dB improvement over other comparison methods on cardiac dataset. Thanks to Nesterov acceleration and learning nonlinear transformations with deeper convolutional layers, FISTA-Net outperforms ISTA-Net+ across all CS ratios. The reconstruction results of all deep unfolding methods are better than traditional CS-MRI methods and patch-based methods. In addition, Table \ref{table_comparison} lists the average testing time (in seconds) of reconstructing a 256 $\times$ 256 image for all methods. Our method achieve consistently better reconstruction performance while maintaining a fast reconstruction real-time speed.

\begin{table*}
\centering
\begin{lrbox}{\mybox}
\setlength{\tabcolsep}{1.5mm}{
\begin{tabular*}{\hsize}{@{}@{\extracolsep{\fill}}lcccccccccccc@{}}       
\toprule[1.5pt]
\multirow{2}*{\textbf{Method}}&\multicolumn{1}{c}{} & \multicolumn{5}{c}{\textbf{Pseudo radial CS Ratio $\alpha$}}&\multicolumn{1}{c}{}&\multicolumn{5}{c}{\textbf{2D random CS Ratio $\alpha$}} \\ \cmidrule[0.7pt](lr){3-7} \cmidrule[0.7pt](lr){9-13}
&  & 10\% & 20\% & 30\% & 40\% & 50\% & & 5\%& 10\% & 20\% & 30\% & 40\% \\
\midrule[0.8pt]
Zero-filling\cite{Bernstein2001}& & 26.64& 30.28&  32.89&  35.01 & 36.92  &   & 25.98& 27.72&  29.34&  31.12&  33.01\\
TV\cite{Lustig2007}& & 30.83 & 35.16 & 38.03 &  40.13 & 41.94  &  & 31.99& 34.34& 36.94& 39.34 & 41.58\\
{SIDWT}& & 30.88 & 35.66 & 38.67 & 40.79& 42.57   &  & 31.79& 34.43& 37.36& 39.95 & 42.23\\
{PBDW}\cite{Qu2012}&  & 32.48 & 36.54 & 38.90& 40.67& 42.30   &  & 33.34&35.67& 38.48& 40.83 & 42.94\\
{PANO}\cite{Qu2014}&  &33.61 &37.02 & 39.28 & 41.05& 42.72   &  &34.31&36.62& 39.57& 42.04& 42.94\\
{{ADMM-Net}\cite{Yang2016}}&  &{---} &{37.17} & {39.84} & {41.56}& {43.00}   &  &{---}&{---}& {---}&{---}&{---}\\
{DC-CNN}\cite{Schlemper2018}&  & 34.30&38.43& 40.74&42.51 & 44.07 &  & 35.01&37.65& 40.80&43.34 & 45.56\\
{ISTA-Net}\cite{Zhang2018}& & 34.69 &38.75 &41.00&42.26& 44.25  & & 35.59&38.23 &41.44&43.99& 46.21\\
{ISTA-Net+}\cite{Zhang2018}&  & 34.83&38.75 &40.99&42.64& 44.22   &  & \underline{35.92}&\underline{38.41} &41.66&44.18& 46.51\\
{FISTA-Net}\cite{Xiang2020}& & \underline{35.09}&\underline{39.05}& \underline{41.20} & \underline{42.85}& \underline{44.42}    &  & 35.92&38.36& \underline{41.74} & \underline{44.27}& \underline{46.69}\\
{CGPD-CSNet}& &\textbf{36.22}&\textbf{39.54}&\textbf{41.52}&\textbf{43.08}& \textbf{44.58}  &  &\textbf{36.72}&\textbf{39.08}&\textbf{42.10}&\textbf{44.52}& \textbf{46.84}\\
\bottomrule[1.5pt]
\end{tabular*}}
\end{lrbox}
\caption{Quantitative assessment with M-PSNR values using pseudo radial and 2D random sampling mask on brain dataset with different CS ratios for CS-MRI. The best and second places are highlighted in Bold font and underlined ones, respectively.}
{\usebox{\mybox}}\label{table_comparison_mask}
\end{table*}

\subsubsection{Qualitative evaluation}
The visual comparisons of all methods using Cartesian sampling mask with CS ratio 10\% on brain CS-MRI and CS ratio 30\% on cardiac CS-MRI are shown in Fig.\ref{fig_results}. We can observe that our method is able to restore more details (e.g., sharper edges, texture) which yield much better visual quality and achieve much higher PSNR scores than the comparison methods.


To take a deeper perspective of the reconstruction performance, we illustrate the violin plot of PSNR values using Cartesian sampling mask with CS ratio 10\% on brain dataset and CS ratio 30\% on cardiac dataset in Fig.\ref{fig_violinplot}. We can observe that compared with the conventional CS-MRI methods, the deep unfolding methods achieve better and more stable reconstruction performance. We remark that our method strikes a good balance between preserving model-based domain knowledge and powerful learning ability, hence it performs significantly better than other methods and has more potential to achieve better performance.


It is interesting and meaningful to qualitatively analyze the presentations of the pre-relaxation module $\mathcal{D}_\ell$ and the geometric prior distillation module $\mathcal{P}_\ell$ in different stages. Fig.\ref{fig_iteration_image} shows the intermediate visualizations of the proposed method at different stages by using Cartesian sampling mask with CS ratio 10\% on brain dataset and CS ratio 30\% on cardiac dataset. We can observe that the details of intermediate results become gradually more complete and that the reconstruction performance become better as iterations increase. It shows that both modules $\mathcal{D}_\ell$ and $\mathcal{P}_\ell$ can facilitate each other for alternating artifact removal and details recovery. These performances are also due to the deep supervision loss, which ensures that the reconstruction results can be steadily improved even with a large number of iterations.

\begin{figure}[t]
\centerline{\includegraphics[width=\columnwidth]{fig/figs/TCGD_TNNLS_Reconstruction_fig05_v1.pdf}}
\caption{{ Stage representations of the proposed method using three sampling masks on brain and cardiac datasets.}}
\label{fig_iteration_image}
\end{figure}


In summary, from the perspectives of visualization and objective
evaluations, our method can achieve a prominent improvement across all CS ratios on different datasets by using widely used Cartesian sampling mask for clinical $k$-space data acquisition. More importantly, such excellent performance can be achieved by reconstructing multiple CS ratios tasks through a single model. Apparently, it can be attributed to the decoupling of prior term and the bonus of two-grid-cycle architecture that quickly generate the low frequency correction in $k$-space.




\begin{figure}[htbp]
\centering
\includegraphics[width=\columnwidth]{fig/figs/TGC_Reconstruction_fig06_v1.pdf}
\caption{Performance comparisons on Cartesian under-sampled k-space datasets with
untrained CS ratios for brain dataset (a) and cardiac dataset (b). "-Left" denotes the reconstruction of untrained CS ratio 15\%, 25\%, 35\%, 45\% sampled data with trained models of CS ratios 10\%, 20\%, 30\%, 40\%, respectively. "-Right" denotes the reconstruction of untrained CS ratio 15\%, 25\%, 35\%, 45\% sampled data with trained models of CS ratios 20\%, 30\%, 40\%, 50\%, respectively.}
\label{fig_unsampled_compare}
\end{figure}

\subsection{Stability of Proposed Method}
 To achieve comprehensive evaluation, we conduct two groups of experiments including CS-MRI with untrained CS sampling ratios and CS-MRI with different sampling masks (pseudo radial sampling mask and 2D random sampling mask).

\subsubsection{CS-MRI with untrained CS ratios}
Fig.\ref{fig_unsampled_compare} depicts the test comparisons on Cartesian under-sampled $k$-space datasets with untrained CS ratios for brain and cardiac MR images. Here all comparison methods
trained five models according to five CS ratios 10\%, 20\%, 30\%, 40\% and 50\%, respectively, while the proposed unifying framework jointly train a single model on the above five CS ratios. Compared with other methods, our unifying framework improves at least 2.49 dB, 1.56 dB, 0.59 dB, and 0.58 dB for CS sampling ratio 15\% ,25\%, 35\%, 45\% on brain dataset respectively, with an average 1.30 dB improvement. And our unifying framework achieves 0.55 dB, 0.76 dB, 0.43 dB improvement for CS ratio 25\%, 35\%, 45\% on cardiac dataset. It is worth noting that our method can use the single multi-sampling-ratio model to reconstruct more satisfactory result from sampling data with any untrained CS ratio than other comparison methods.


\subsubsection{CS-MRI with different sampling masks}
Table \ref{table_comparison_mask} lists the M-PSNR values of various methods by using pseudo radial sampling mask and 2D random sampling mask with different CS sampling ratios on brain dataset. We can see that the proposed method achieves optimal reconstruction results across all CS ratios. 
In general, it can achieve average 0.47 dB improvement using pseudo radial sampling mask and average 0.45 dB improvement using 2D random sampling mask. Furthermore, we can observe that our method achieves much better visual quality and higher PSNR than other comparison methods in Fig.\ref{fig_results},
and that the details of intermediate results become gradually better as iterations increase in Fig.\ref{fig_iteration_image}
 using pseudo radial sampling mask with CS ratio 10\% and 2D random sampling mask with CS ratio 10\% on brain dataset.


\section{Conclusions}
In this study, we proposed a novel deep unfolding unified framework that can deliver generalizability by flexibly handling multi-sampling-ratio CS-MRI through a single model. Inspired by efficient correction of multigrid technique, we start from classical CS-MRI optimization problem and transform it into a two-grid-cycle correction architecture which consists of pre-relaxation, correction and geometric prior distillation. The proposed method inherits the merits of model-based and DL-based CS-MRI methods, has sufficient theoretical support and also good performance. It can give us a new perspective to design the explainable network. Extensive experiments demonstrate that the proposed framework outperforms other state-of-the-art methods in terms of qualitative and quantitative evaluations.



\bibliographystyle{IEEEtran}
