Joint MR sequence optimization beats pure neural network approaches for
  spin-echo MRI super-resolution
Current MRI super-resolution (SR) methods only use existing contrasts
acquired from typical clinical sequences as input for the neural network (NN).
In turbo spin echo sequences (TSE) the sequence parameters can have a strong
influence on the actual resolution of the acquired image and have consequently
a considera-ble impact on the performance of the NN. We propose a
known-operator learning approach to perform an end-to-end optimization of MR
sequence and neural net-work parameters for SR-TSE. This MR-physics-informed
training procedure jointly optimizes the radiofrequency pulse train of a proton
density- (PD-) and T2-weighted TSE and a subsequently applied convolutional
neural network to predict the corresponding PDw and T2w super-resolution TSE
images. The found radiofrequency pulse train designs generate an optimal signal
for the NN to perform the SR task. Our method generalizes from the
simulation-based optimi-zation to in vivo measurements and the acquired
physics-informed SR images show higher correlation with a time-consuming
segmented high-resolution TSE sequence compared to a pure network training
approach.
Hoai Nam Dang, Vladimir Golkov, Thomas Wimmer, Daniel Cremers, Andreas Maier, Moritz Zaiss
Joint MR sequence optimization beats pure neural network approaches for spin-echo MRI super-resolution

###### Abstract

Current MRI super-resolution (SR) methods only use existing contrasts acquired from typical clinical sequences as input for the neural network (NN). In turbo spin echo sequences (TSE) the sequence parameters can have a strong influence on the actual resolution of the acquired image and have consequently a considerable impact on the performance of the NN. We propose a known-operator learning approach to perform an end-to-end optimization of MR sequence and neural network parameters for SR-TSE. This MR-physics-informed training procedure jointly optimizes the radiofrequency pulse train of a proton density- (PD-) and T2-weighted TSE and a subsequently applied convolutional neural network to predict the corresponding PDw and T2w super-resolution TSE images. The found radiofrequency pulse train designs generate an optimal signal for the NN to perform the SR task. Our method generalizes from the simulation-based optimization to in vivo measurements and the acquired physics-informed SR images show higher correlation with a time-consuming segmented high-resolution TSE sequence compared to a pure network training approach.

Keywords:super-resolution, turbo spin echo, joint optimization

## 1 Introduction

Magnetic resonance imaging plays an essential role in clinical diagnosis by acquiring the structural information of biological tissue. Spatial resolution is a crucial aspect in MRI for the precise evaluation of the acquired images. However, there is an inherent trade-off between the spatial resolution of the images and the time required for acquiring them [1]. In order to obtain high-resolution (HR) MR images, patients are required to remain stable in the MR scanner for long time, which leads to patients' discomfort and inevitably introduces motion artifacts that again compromise image quality and actual resolution [2]. Since super-resolution (SR) can improve the image quality without changing the MRI hardware, this post-processing tool has been widely used to overcome the challenge of obtaining HR MRI scans [3]. Using model-based methods like interpolation algorithms [4] and iterative deblurring algorithms [5] or learning-based
methods such as dictionary learning [6], SR achieved the restauration of fine structures and contours. In recent years, deep learning has become a main-stream approach for super-resolution imaging, and a number of neural network-based SR models were proposed [7]. Among the proposed model- or learning-based methods, convolutional neural networks (CNN) produce superior SR results with better clarity and less artifacts [8]. Super-resolution for MRI data has been only recently applied [9-12]. In [10] a CNN is proposed for cardiac MRI to estimate an end-to-end non-linear mapping between the upscaled low-resolution (LR) images and corresponding HR images to rebuild a HR 3D volume. In other work, motion compensation for the fetal brain was achieved by a CNN architecture [11] to solve the 3D re-construction problems. SR MRI has been also applied to low-field MR brain imaging [12]. However, these existing methods only used single contrast MRI images and did not make full use of multi-contrast information. In the clinical routine, T1, T2 and PD weighted images are often acquired together for diagnosis with complementary information. Although each weighted image highlights only certain types of tissues, they reflect the same anatomy, and can provide synergy when used in combination [8].

Fast imaging techniques like Turbo-Spin-Echo (TSE) [13] sequences can also be utilized to sample more data in given timeframe, thus allowing a higher resolution. However, due to the long echo-train duration the T2-decay is significant during the signal acquisition. This process acts as a voxel-T2-dependent k-space filter that lowers the actual resolution w.r.t the nominal resolution due to a broadening of the point-spread-function (PSF) [14]. However, by adjusting the refocusing radiofrequency (RF) pulses, the signal decay can be reduced during the TSE echo-train [15]. The RF pulse train strongly influence the signal dynamic in a highly complex fashion, as each RF pulse affects all future signal contributions. Current MRI super-resolution methods use contrasts acquired from typical clinical protocols as input for the neural network and disregard the influence of the MR sequence parameters for optimization. Using so-called known operator learning [16], we propose an approach that utilizes a MR physics model during the optimization to not only train a neural network for super-resolution, but also adapt the refocusing RF pulses to directly influence the PSF. This approach also allows the use of the uncorrupted theoretical contrast as ground truth, which is only available during the simulation. By using two different encoding schemes in our sequences, we gain additional information from the two different contrasts PD and T2 that are used as input for the CNN and both will have different PSFs, thus provide valuable information for the SR task. Both sequences are optimized jointly to allow generation of optimal contrasts for the SR task of the neural network.

The main contribution and the novelty of our work is the end-to-end optimization of MR sequence and neural network parameters for super-resolution TSE. For this purpose, we use a fully differentiable Bloch simulation embedded in the forward propagation to jointly optimize the RF pulse train of proton density (PD) and T2 weighted TSE sequences and a subsequent applied convolutional neural network to predict the corresponding PDw and T2w super-resolution TSE images. The ground truth targets are directly generated by the simulation and represent the uncorrupted MR contrast. Our
jointly optimized approach is compared to a network trained on a TSE with a 180\({}^{\circ}\) RF pulse train. The optimized sequences and networks are verified at the real scanner system by performing in vivo measurements of a healthy subject and compared to a highly segmented, high-resolution vendor-provided sequence.

## 2 Theoretical Background

### Image Degradation in TSE-Sequences.

When there is no relaxation decay during the echo-train, the k-space signal obtained from a TSE pulse sequence \(S(k_{x},k_{y})\) yields the true spatial distribution of the theoretical transverse magnetization \(M_{\perp}(x,y)\) via the Fourier transform (FT):

\[M_{\perp}(x,y)=\int\int_{k_{x},k_{y}}S\big{(}k_{x},k_{y}\big{)}\cdot e^{i\big{(} k_{x}x+k_{y}y\big{)}}dk_{x}dk_{y} \tag{1}\]

When considering the signal relaxation behavior during acquisition an additional filtering function in k-space, the Modulation Transfer Function (MTF) for each tissue type (tt), has to be applied:

\[\widetilde{M}_{\perp}(x,y) =\sum\nolimits_{tt}M(x,y)\] \[=\sum\nolimits_{tt}\int_{k_{x},k_{y}}S_{tt}\big{(}k_{x},k_{y} \big{)}\cdot MTF_{tt}\big{(}k_{x},k_{y}\big{)}\cdot e^{i\big{(}k_{x}x+k_{y}y \big{)}}dk_{x}dk_{y}\] \[=\sum\nolimits_{tt}\int_{k_{x},k_{y}}S_{tt}\big{(}k_{x},k_{y} \big{)}\cdot e^{i\big{(}k_{x}x+k_{y}y\big{)}}dk_{x}dk_{y}\] \[\qquad\qquad*\int_{k_{x},k_{y}}MTF_{tt}\big{(}k_{x},k_{y}\big{)} \cdot e^{i\big{(}k_{x}x+k_{y}y\big{)}}dk_{x}dk_{y}\] \[=\sum\nolimits_{tt}M_{\perp,tt}(x,y)\cdot B_{tt}(x,y),\]

where \(*\) denotes a convolution and \(B(x,y)\) is a blur kernel equal to the Fourier-transformed MTF. For a single-shot 180\({}^{\circ}\) TSE sequence with constant RF pulses the filtering function in k-space can be described as:

\[MTF_{tt}\big{(}k_{x},k_{y}\big{)}=\sum\nolimits_{tt}e^{-\frac{t\big{(}k_{x}x \cdot k_{y}\big{)}}{T2tt}}\,,\qquad\quad(2)\]

which is unique for each tissue type with different T2 value. Using variable RF pulses, the MTF can become more homogeneous across the k-space and therefore reduce the width of the PSF.
Methods

### Sequences

A single-shot 2D TSE sequence is being used as default sequence for our optimization. The acquisition time for the single-shot 2D TSE is 0.76 s at 1.56 mm in-plane resolution. Single slice acquisition was used for all sequences. The refocusing RF pulses of the PDw TSE with TE=12 ms and T2w TSE with TE=96 ms were optimized jointly. The PDw TSE sequence uses a centric phase-encoding reordering. For T2w imaging the centric phase-encoding reordering is shifted to have the central k-space line encoded at the given echo time TE, which for TE=96 ms is at the 8th echo. Other parameters were as follows: acquisition matrix of 128\(\times\)128, undersampling factor in phase: 2x, reconstructed with GRAPPA[17], FOV=200 mm\(\times\)200 mm, slice thickness of 8 mm and bandwidth of 133 Hz/pixel. For all sequences, the 90\({}^{\circ}\) excitation pulse was kept fixed.

### Simulation & Optimization

All simulations and optimizations were performed in a fully differentiable Bloch simulation framework [18]. The framework generates MR sequences and corresponding reconstruction automatically based on the target contrast of interest. The optimization is carried out in an MR scanner simulation environment mirroring the acquisition of a real MR scanner. The forward simulation consists of a chain of tensor-tensor multiplication operations, representing the Bloch equations, that are differentiable in all parameters and supports an analytic derivative-driven nonlinear optimization. The entire process - MRI sequence, reconstruction, and evaluation - is modelled as one computational chain and is part of the forward and backwards propagation during the optimization, as depicted in Figure 1. The optimization problem is described by:

\[\Psi^{*},\Theta^{*}=\operatorname*{argmin}_{\Psi,\Theta}\left(\sum_{i}\left\| M_{\lambda,i}-\operatorname*{NN}_{\Theta}\left(RECO\left(\operatorname{SCAN}_{\Psi}(P_{i}) \right)\right)\right\|_{p}\right), \tag{3}\]

where \(\Psi\) are the optimized sequence parameters and \(\Theta\) the neural network parameters. For given tissue maps \(P_{i}\) for each Voxel \(i\) the Bloch simulation \(SCAN\) outputs the MR signal and is reconstructed by the algorithm \(RECO\). Signal simulations were performed by a fully differentiable extension of the Extended Phase Graph (EPG) [19] formalism. The simulation is done with PyTorch [20] complex-valued datatype and outputs a complex-valued signal. The forward simulation outputs the TSE signal which is conventionally reconstructed to magnitude images, and in addition the corresponding contrast as ground truth target as given in equation [1]. For the SR network a CNN, DenseNet [21] was adapted, which receives the magnitude TSE images of PDw and T2w TSE as input. To prevent scaling discrepancy, the TSE images are normalized to have maximum value of 1 before applying the CNN for both cases in simulation and in vivo. The DenseNet consist of 4 Dense blocks (Convolution-\(>\)BatchNorm\(>\)PReLu\(>\)Concat) followed by an UpsampleBlock (bicubic upsampling-\(>\)Convolution-\(>\)BatchNorm\(>\)PReLu) and a final CNN Layer. Each convolution had a 3\(\times\)3 kernel size, except for the first layer with a kernel size of 7\(\times\)7. In total, the model had 174,706 trainable
parameters. Of the sequence parameters \(\Psi\), the amplitude of the refocusing RF pulses of the TSE sequences and the NN parameters \(\Theta\) were optimized jointly in an end-to-end training procedure using the Adam optimizer [22] in PyTorch. We follow a known-operator approach [16], where the conventional reconstruction, including parallel imaging by means of GRAPPA, is fixed, but fully differentiable. The simulation is fully differentiable and all parameters except the refocusing RF pulses are fixed. The gradient update propagates back through the whole chain of differentiable operators. The complete RF pulse train and CNN are updated at each iteration step. The RF pulses are initialized with random values around 50\({}^{\circ}\) and standard deviation of 0.5\({}^{\circ}\).

The training data consisted of synthetic brain samples based on the BrainWeb [23] database. The fuzzy model segments were filled with in vivo-like tissue parameters: Proton density PD values were taken from [24], T1 and T2 from [25], T2' was calculated from T2 and T2* values [26] and diffusion coefficient D was taken from [27]. B0 and B1 were assumed to be without inhomogeneities. In total, 19 subject volumes each consisting of 70 slices were used as training data and one separate subject volume as test dataset. The simulation uses coil sensitivity maps acquired at the MR system and calculated using ESPIRiT [28]. The optimizations were performed on an Intel Xeon E5-2650L with 256GB RAM. A full optimization on CPU took 4 days with memory consumption of 230GB RAM. The learning rate for the model parameters and sequence parameters were lr_model=0.001, lr_rf=0.01, respectively. Other hyperparameters of the optimization were: batch size = 1, n_epoch=10, damping factors of Adam (0.9, 0.999).

### Data acquisition at a real MR system

After the optimization process, all sequences were exported using the Pulseq standard [29] and the populseq tool [30]. Pulseq files could then be interpreted on a real MRI scanner including all necessary safety checks, and were executed on a PRISMA 3T scanner (Siemens Healthineers, Erlangen, Germany) using a 20-channel head coil. Raw

Figure 1: Overview of the proposed processing pipeline: The MR signal of PDw and T2w TSE is simulated for a given RF pulse train; GRAPPA reconstruction and SR CNN are applied subsequently. The output is compared to the actual theoretical uncorrupted HR contrasts at TE\({}_{\text{eff}}\) and a gradient descent is performed to update refocusing FA and NN parameters, simultaneously.
data were automatically sent back to the terminal and the same reconstruction pipeline used for the simulated data was used for measured images. As high-resolution reference a vendor-provided TSE sequence was acquired with following parameters: 32-shot segmented, GRAPPA2, TE=12/96 ms, TR=12 s, FOV=200 mm\(\times\)200 mm, matrix of 256\(\times\)256, FA=180\({}^{\circ}\). All MRI scans were under approval of the local ethics board and were performed after written informed consent was obtained. Measurements were performed on a healthy volunteer.

### Reconstruction and Evaluation

Signals of the TSE sequences were reordered, and reconstructed with GRAPPA. The optimization was based solely on magnitude images. Structural similarity index measure (SSIM) [31] and the Peak Signal-to-Noise Ratio (PSNR) were calculated for evaluation of simulation and in vivo measurements w.r.t. the simulated ground truth and the HR segmented in vivo measurement, respectively. The evaluation was performed in Matlab [32] with the build-in functions for SSIM and PSNR.

## 4 Results

### Qualitative Visual Results

The original LR TSE image with the zero-filled image and the reconstructed SR image are compared to our optimized RF pulse train design and a conventional 180\({}^{\circ}\) RF pulse train TSE sequence for each contrast in Figure 2. The optimization process can be seen in Supporting Figure S1 and Supporting Animation S3. Starting from the initialized values, the RF pulses converge to the optimal RF pulse train, while the NN parameters are optimized simultaneously. The converged RF pulse state has been found to be independent from the initialization. The final optimized RF pulse design for the PDw and T2w TSE sequences is shown in Figure 2a. It can be observed, that in all cases the SR image leads to an improvement over the LR TSE image by showing clearer resolved borders between white and gray matter. The optimized RF pulse train further improves the nominal resolution, which can be observed by a clear increase of sharpness of the sulcus between Gyrus cinguli and Gyrus frontalis superior as indicated by the red arrows. The optimized sequence and CNN translate well to in vivo measurements, where similar improvements as seen in the simulated images can be observed (Figure 2d,e).

### Quantitative Metrics Results

Table 1 and Table 2 report the quantitative metrics scores of PSNR and SSIM for the images shown in Figure 2. The quantitative metrics agree with our visual observations and show that our end-to-end optimization approach performs better than the SR based on existing conventional 180\({}^{\circ}\) TSE sequence data only. Compared to the acquisition time of the segmented reference sequence with 192.85s, our optimized single-shot sequence only requires an acquisition time of 0.76s. Thus, the SR performance can
potentially be further increased by a multi-shot scan sacrificing a little more time. To find the best procedure multiple ablation studies were performed (Supporting Information).

\begin{table}
\begin{tabular}{|l|l|l|l|l|} \hline  & PSNR - PDw & SSIM - PDw & PSNR - T2w & SSIM – T2w \\ \hline LR – 180\({}^{\circ}\) & 10.4 & 0.39 & 11.6 & 0.40 \\ \hline ZF – 180\({}^{\circ}\) & 14.7 & 0.73 & 22.2 & 0.79 \\ \hline HR NN – 180\({}^{\circ}\) & 28.7 & 0.93 & 29.4 & 0.93 \\ \hline LR – opt. FA & 10.8 & 0.38 & 11.2 & 0.39 \\ \hline ZF – opt. FA & 17.5 & 0.75 & 23.6 & 0.80 \\ \hline \end{tabular}
\end{table}
Table 1: PSNR and SSIM for simulation in Figure 2 using the high-resolution GT as reference.

Figure 2: (a) Optimized RF pulse train and phase encoding for both contrast (due to centric reordering, the central k-space line k_y=0 is acquired at repetition 0 and 7 for the PDw and T2w TSE, respectively). (b) Simulation of a static 180\({}^{\circ}\) RF pulse train and (c) optimized RF pulse train compared to the uncorrupted ground truth. (d) In vivo measurements of a static 180\({}^{\circ}\) RF pulse train and (e) optimized RF pulse train compared to the vendor’s TSE sequence is shown as a high-resolution reference. In both cases the improvement of the optimized RF pulses over the constant 180\({}^{\circ}\) RF pulse train can be observed by a better resolved border between white and gray matter as indicated by the red arrows. SSIM and PSNR values are shown in Table 1 and 2.
## 5 Discussion

We demonstrated a new end-to-end learning process for TSE super-resolution by jointly optimizing refocusing RF pulse trains and neural network parameters. This approach utilizes a differentiable MR physics simulation embedded in the forward and backward propagation. The joint-optimization outperforms a pure neural network training. Although our approach is solely based on simulated data, the optimized sequence and trained CNN translate well to in vivo data. By using simulation-based training data, we are able to use the theoretical uncorrupted contrast as ground truth target. Apart from the expensive acquisition of HR in vivo data, real measured target data also have inherent drawbacks compared to their LR counterpart. Due to the longer scan time motion artifacts become more significant and to acquire the same contrast the bandwidth has to be increased, leading to a decrease of SNR [33]. However, we also admit the limitation of a simulation-based optimization, as the performance is bound to the accuracy of the model behind the simulation. We can observe in our results, that our approach is not able to resolve small vessel structures, as these are not existing in our synthetic brain database. Fortunately, the NN does not hallucinate details, when encountering these structures. Using real measured data to finetune the trained network could be a possible solution for this problem. Another way could be including uncertainty quantification layers [34] in the CNN to handle unknown structures. Our approach is compatible with any network architecture e.g. [35-37] to further improve the SR task. Furthermore, the training objective can be also extended to requirements on the MR sequence by including constraints in the loss function e.g. reduced RF pulse amplitudes for decrease of energy deposition SAR or increase of SNR.

To conclude, we propose an end-to-end optimization of MR sequence and neural network parameters for TSE super-resolution. This flexible and general end-to-end approach benefits from a MR physics informed training procedure, allowing a simple target-based problem formulation, and outperforms pure neural network training.

\begin{table}
\begin{tabular}{|l|l|l|l|l|} \hline  & PSNR - PDw & SSIM - PDw & PSNR - T2w & SSIM – T2w \\ \hline LR \(-180^{\circ}\) & 12.6 & 0.35 & 16.2 & 0.38 \\ \hline ZF \(-180^{\circ}\) & 14.1 & 0.59 & 22.5 & 0.70 \\ \hline HR NN – \(180^{\circ}\) & 23.2 & 0.84 & 26.2 & 0.89 \\ \hline LR – opt. FA & 13.2 & 0.37 & 15.9 & 0.38 \\ \hline ZF – opt. FA & 15.5 & 0.65 & 22.8 & 0.70 \\ \hline HR NN – opt. FA & **26.0** & **0.87** & **27.0** & **0.90** \\ \hline \end{tabular}
\end{table}
Table 2: PSNR and SSIM for in vivo measurements in Figure 2 using the segmented high-resolution TSE as reference.