PS-Net: Learned Partially Separable Model for Dynamic MR Imaging

Based on PS model, we proposed an adaptive module to characterize the low rankness of dynamic MRI. According to the classical Prony’s results[25], there exists a $L$-tap filter $h_{ps}(r,t)$ annihilates the PS model $\gamma(r,t)$:

where $h_{ps}(r,t)$ has the following form:

By the equivalence relation of Hankel matrix multiplication and convolution, (4) is equivalent to:

where $H(\gamma)$ is a linear convolution matrix with Hankel-structure in the $N_{t}$ direction and $h_{ps}(r,t)$ is also called the null space filter. Then the low rankness can be implicitly characterized via this null space filter.

Since the signal intensities of adjacent pixels in dynamic MR images are close to each other, the sparse property is satisfied after performing the gradient on the image[29]:

where $\mu$ is a set of bandlimited functions, $\mu\in\mathbb{C}^{N_{x}\times N_{y}}$. By the convolution theorem, (7) can be equivalently expressed as:

where $\widehat{\nabla\gamma}$ represents the Fourier coefficients of $\nabla\gamma(r)$ (along the $N_{x}N_{y}$ direction), $w$ is the coordinate in frequency domain, $h_{s}(w)$ is the frequency domain signal obtained by the Fourier transform of $\mu(r)$. By the equivalence of convolution and Hankel matrix multiplication, (8) can be rewritten as:

where $H(\widehat{\nabla\gamma}(w))$ is the Hankel-structured matrix in spatial dimension.

With the above PS annihilation module and sparse module, we formulate the dynamic MRI reconstruction problem as the following unconstrained two-norm minimization problem:

where $||H_{1}(\widehat{\nabla\gamma_{i}})h_{si}||^{2}_{F}$ represents the sparse regularizer, $T$ is the number of time frames, and $||H_{2}(\gamma)h_{ps}||^{2}_{F}$ is the null space annihilation regularizer, which is derived by the partially separable property. $\widehat{\nabla\gamma}=[\widehat{\nabla\gamma_{1}},\widehat{\nabla\gamma_{2}},...,\widehat{\nabla\gamma_{T}}],h_{s}=[h_{s1},h_{s2},...,h_{sT}]$, $\lambda_{1}$ and $\lambda_{2}$ are the corresponding regularization parameters. To solve the optimization problem in (10), we introduce auxiliary variables $x_{i}$, $\widehat{U_{i}}$, $Z$. The reformed optimization problem is as follows:

By substituting $\widehat{\nabla\gamma_{i}}=\left(\begin{array}[]{l}iw_{x}F\left(\gamma_{i}\right)\\ iw_{y}F\left(\gamma_{i}\right)\end{array}\right)$, where $w_{x}$ and $w_{y}$ denote the k-space coordinates in $x$ and $y$ direction, respectively. The penalty function of (11) can be written as:

Based on the semi-quadratic splitting method (HQS), we can construct an iterative solution algorithm for the above penalty function:

By the first-order optimization condition, it reads:

$\gamma$ can be obtained by iteratively solving the four sub-problems of (14), but iteratively solving (14) directly requires adjusting many parameters ($\lambda_{1},\lambda_{2},\rho_{1},\rho_{2}$).

We solve the above problem by unrolling (14) into PS-Net. The above four sub-problems correspond to three modules in the network: sparse module $S$ (corresponding to the $\widehat{U}_{i}^{k}$ sub-problem), PS Annihilation module $PS$ (corresponding to the $Z^{k}$ sub-problem), and reconstruction module $R$ (corresponding to the $x_{i}$ and $\gamma^{k}$ sub-problems). The overall structure of PS-Net is shown in Fig. 1. There are 10 iterations in the network. The regularization parameters are learned by the network adaptively. The detailed structures of the three models in the network is shown below.

• •

  The sparse module $S$: Since the Hankel matrix multiplication and the convolution are equivalent, the $\widehat{U}_{i}^{k}$ sub-problem is solved via the convolution layer in the network. The convolution layer is set to 2D convolution along the $N_{x}$ and $N_{y}$ directions because the gradient image of dynamic MRI satisfies the sparse property along these two directions. The number of convolution layers is 5 and the number of convolution channels is 64. We split the real and imaginary parts of the input data into two channels, which means that the dimension size in the time direction changes from $T$ to $2T$.

• •

  The PS annihilation module $PS$: The PS annihilation layer characterizes the low-rank prior by a learnable null space filter mapping. Similarly, $Z^{k}$ is implemented equivalently through the convolution layers. Considering that the image domain has redundancy [29], we consider not only the convolution in the time direction but also along the $N_{x}$ and $N_{y}$ directions. The numbers of convolution layers and convolution channels are the same as the sparse module. The input data is also divided into real and imaginary channels.

• •

  The reconstruction module $R$: According to the first-order optimization condition, the minimal value of the energy function with respect to $\gamma$ can be solved by (14), thus ensuring data consistency.

The pseudo-code for the network unrolling solution is shown in Alg. 1.

PS-Net reconstructs the image via supervised learning. The label of PS-Net is the fully sampled dynamic MR image. Input data is generated by artificially adding the undersampling mask to the label. The network’s output is the reconstructed dynamic MR image, using the minimum mean square error as the loss function, which can be expressed as:

where $\gamma^{*}$ is the output, $\gamma_{ref}$ is the label, $N_{samples}$ is the total number of samples. The activation function in the network is handled by Relu[30], and Adam has been selected as the optimizer[31]. All hyper-parameters in the network are initialized to 1 ($\lambda_{1},\lambda_{2},\rho_{1},\rho_{2}$). We set the epoch of network training to 50. The initial value of the learning rate is 0.001; the learning rate delay is 0.95; the experimental environment is tensorflow2.7[32], cuda11.3, ubuntu20.04, and the graphics card is NVIDIA A100 Tensor Core GPU.

The cardiac cine dataset was acquired from 30 individuals using a 3T MRI scanner (Trio, Siemens, 20-channel receiver coil). The Ethics Review Committee of the Shenzhen Institute of Advanced Technology, Chinese Academy of Sciences, authorized all scanning studies with the subjects’ informed consent. The imaging parameters were: acquisition matrix = 256$\times$256, slice thickness = 6mm, TR/TE = 3.0ms/1.5ms, and FOV = 330$\times$330mm. Total 356 slices were acquired. PS-Net was performed in both single-channel and multi-channel data, where the single-channel data was obtained by fusing the sensitivity of 20 coils and the coil sensitivity in (10) is set to the identity matrix for the single-channel scenario. Data from 25 individuals were randomly selected as the training dataset, and the remaining 5 individuals were used as the test dataset. We cropped the dynamic images along the $x$, $y$, and $t$ direction and used rigid transform cropping for data enhancement. Eventually, the training dataset had 800 images, and the test dataset had 118 images, with a size of 192$\times$192$\times$18.

The data was undersampled by applying undersampling masks on the fully sampled k-space data. Random cartesian undersampling masks[33, 34] were employed for the fully sampled data. The phase coding direction is undersampled and the frequency coding direction is fully sampled. The input and fully sampled images are coupled one by one as the training pairs.

The quantitative evaluation metrics include minimum mean square error (MSE), peak signal-to-noise ratio (PSNR), and structural similarity (SSIM):

where $\gamma^{*}$ is the output and $\gamma_{ref}$ is the label, $l,c,s$ are three comparison measurements between the samples of $\gamma^{*}$ and $\gamma$[35]. The smaller MSE, the larger PSNR and SSIM, indicate a better reconstruction.

In the single-channel scenario, we compared PS-Net with the traditional compressed sensing method L+S[36], the deep learning end-to-end method CRNN[37], and the model-driven based method SLR-Net[14].

The reconstruction results of the single-channel data undersampled by 8-fold acceleration are shown in Fig. 2. From top to bottom (except the y-t image), Fig. 2 shows the reconstruction results of the 1st frame, the 10th frame, and the 18th frame to represent the dynamic images both in diastole and systole. The first row of each frame shows the reconstruction results of L+S, CRNN, SLR-Net, and PS-Net. The second row shows the zoomed-in view of the ROI region (indicated by the yellow box in the first row), and the third row shows the error map of the ROI region. The y-t image in the bottom shows the 92nd slice along the $y$ and temporal direction. The error range of the error map is [0–0.09].

The enlarged ROI region in the second row of the 1st and the 18th frame (the diastolic phase) suggest that PS-Net clearly reconstructs the tiny structures of papillary muscle without artifacts and outperforms other methods in detail retention and artifact removal. The error map of the reconstruction in the fifth row along the temporal dimension shows that PS-Net also outperforms the other methods in the dynamic map reconstruction. The reconstructed images of the 10th frame (the systolic phase) are not as good as the 1st and the 18th frame. Although all images are blurred due to high undersampling rate in the 10th frame, the PS-Net achieves the lowest errors with no significant deformation appears in the reconstructed results.

Compared with SLR-Net, which learns low-rank by the singular value soft threshold (SVT), we found that PS-Net is more accurate for low-rank learning, indicating that our adaptive low-rank learning strategy outperforms the representing low-rank prior by the hand-crafted nuclear norm.

We also performed the experiments at higher acceleration rates with 10-fold and 12-fold. The results are shown in Fig. 3. PS-Net clearly reconstructs the tiny structures of papillary muscle even at 12-fold acceleration. The quantitative evaluation metrics are shown in Table I, and it can be seen that the performance of PS-Net outperforms the state-of-the-art methods.

In the multi-channel scenario, comparisons were made with the L+S method[36], and the model-driven based methods MoDL[15] and SLR-Net[14]. Because multi-channel data contains redundant information such as coil sensitivity, multi-channel data can be accelerated to a higher rate. The reconstruction results of multi-channel data undersampling with 12-fold acceleration are shown in Fig. 4, and the quantitative evaluation metrics are shown in Table II. In the multi-channel case, the reconstruction result of PS-Net outperforms the L+S, MoDL, and SLR-Net in terms of detail description, artifact removal, etc.

We designed ablation experiments to verify that adaptive low-rank learning strategy in PS-Net is effective. Three network were performed. The first one uses the SVT strategy to depict the low rank, implemented as the SVT module, which replaces the PS annihilation module, and the sparse module remains unchanged. This network is referred to as PS-Net (SVT). The second one removes the PS annihilation module, and the other modules and initialized hyper-parameters were kept unchanged, called Sparse-Net. The third one removes the sparse module and still uses the PS annihilation module for low-rank adaptive learning, named LowRank-Net. The initialization hyper-parameters of the above three networks are the same as PS-Net. The ablation experiments were performed in the multi-channel scenario with 12-fold acceleration. The reconstruction results are shown in Fig. 5. PS-Net reconstructs the papillary muscle more clearly and more detailed than PS-Net (SVT), without artifacts. And the error map of PS-Net in third row and fifth row are clearer than PS-Net (SVT), which indicates that our adaptive low-rank learning strategy effectively portrays the low-rank of dynamic MR images more accurately. The reconstruction effect of LowRank-Net is better than that of PS-Net(SVT), Sparse-Net, which also shows that our adaptive low-rank learning strategy is effective. After adding the sparse prior, the reconstruction result of PS-Net achieves the best result among the four networks. Table III shows the quantitative metrics of the reconstruction.