Complex-valued neural networks to speed-up MR Thermometry during
  Hyperthermia using Fourier PD and PDUNet
Hyperthermia (HT) in combination with radio- and/or chemotherapy has become
an accepted cancer treatment for distinct solid tumour entities. In HT, tumour
tissue is exogenously heated to temperatures between 39 and 43 $^\circ$C for 60
minutes. Temperature monitoring can be performed non-invasively using dynamic
magnetic resonance imaging (MRI). However, the slow nature of MRI leads to
motion artefacts in the images due to the movements of patients during image
acquisition. By discarding parts of the data, the speed of the acquisition can
be increased - known as undersampling. However, due to the invalidation of the
Nyquist criterion, the acquired images might be blurry and can also produce
aliasing artefacts. The aim of this work was, therefore, to reconstruct highly
undersampled MR thermometry acquisitions with better resolution and with fewer
artefacts compared to conventional methods. The use of deep learning in the
medical field has emerged in recent times, and various studies have shown that
deep learning has the potential to solve inverse problems such as MR image
reconstruction. However, most of the published work only focuses on the
magnitude images, while the phase images are ignored, which are fundamental
requirements for MR thermometry. This work, for the first time, presents deep
learning-based solutions for reconstructing undersampled MR thermometry data.
Two different deep learning models have been employed here, the Fourier
Primal-Dual network and the Fourier Primal-Dual UNet, to reconstruct highly
undersampled complex images of MR thermometry. The method reduced the
temperature difference between the undersampled MRIs and the fully sampled MRIs
from 1.3 $^\circ$C to 0.6 $^\circ$C in full volume and 0.49 $^\circ$C to 0.06
$^\circ$C in the tumour region for an acceleration factor of 10.
Rupali Khatun, Soumick Chatterjee, Christoph Bert, Martin Wadepohl, Oliver J. Ott, Rainer Fietkau, Andreas Nürnberger, Udo S. Gaipl, Benjamin Frey
# Fourier PD and PDUNet: Complex-valued networks to speed-up MR Thermometry during Hyperthermia

###### Abstract

Hyperthermia (HT) in combination with radio- and/or chemotherapy has become an accepted cancer treatment for distinct solid tumour entities. In HT, tumour tissue is exogenously heated to temperatures of 39 to 43 \({}^{\circ}\)C for 60 minutes. Temperature monitoring can be performed noninvasively using dynamic magnetic resonance imaging (MRI). However, the slow nature of MRI leads to motion artefacts in the images due to the movements of patients during image acquisition time. By discarding parts of the data, the speed of the acquisition can be increased - known as Undersampling. However, due to the invalidation of the Nyquist criterion, the acquired images have lower resolution and can also produce artefacts. The aim of this work was, therefore, to reconstruct highly undersampled MR thermometry acquisitions with better resolution and with less artefacts compared to conventional techniques like compressed sensing. The use of deep learning in the medical field has emerged in recent times, and various studies have shown that deep learning has the potential to solve inverse problems such as MR image reconstruction. However, most of the published work only focusses on the magnitude images, while the phase images are ignored, which are fundamental requirements for MR thermometry. This work, for the first time ever, presents deep learning based solutions for reconstructing undersampled MR thermometry data. Two different deep learning models have been employed here, the Fourier Primal-Dual network and Fourier Primal-Dual UNet, to reconstruct highly undersampled complex images of MR thermometry. It was observed that the method was able to reduce the temperature difference between the undersampled MRIs and the fully sampled MRIs from 1.5 \({}^{\circ}\)C to 0.5\({}^{\circ}\)C.

keywords: Hyperthermia, MRI, MR Image Reconstruction, Deep Learning, Undersampled MRI, Undersampled, MR Reconstruction, Complex Image, MR Thermometry Msc: [2022] 00-01, 99-00 +
Footnote †: journal:

## 1 Introduction

Hyperthermia (HT) has become one of the well-accepted cancer treatments, in combination with radio / chemotherapy. In HT, tumour tissue is exogenously heated to temperatures of 39 to 44 \({}^{\circ}\)C for 60 minutes to sensitise tumour cells for chemo- and/or radiotherapy (Chorric et al., 2015; Datta et al., 2019; Van der Zee, 2002; Kok et al., 2020). Temperature monitoring is an important part of quality-controlled HT and can be performed noninvasively by Magnetic Resonance Imaging (MRI). However, a major challenge is that the MRI is an inherently slow process. Consequently, the scan time for high-resolution imaging is long, compromising with the temporal resolution. Longer scan times can also lead to an increase in motion artefacts due to the movements of the patient during image acquisition time. The speed of image acquisition can be increased by discarding parts of the data, known as undersampling (Chatterjee et al., 2022). But that leads to loss of resolution and can also produce artefacts due to invalidation of the Nyquist criterion (Nyquist, 1928; Shannon, 1949) Hence, MRI image reconstruction and reduction of motion artefacts are in high demand. This work aims to reconstruct highly undersampled MR thermometry acquisitions of patients with sarcoma with better resolution and with less artefacts compared to conventional techniques such as compressed sensing. The use of deep learning in the medical field is spreading, including for undersampled MR reconstruction. Using the ReconResNet model as the network backbone, the NCC1701 pipeline has been revealed to be capable of removing artefacts from highly undersampled images (Chatterjee et al., 2022) with acceleration factors as high as 20. However, this work does not only focus on the magnitude images, while the phase images are ignored, which are fundamental requirements for MR thermometry. (Chatterjee et al., 2022).

### Thermal Therapy and Thermometry

One of the oldest and easiest thermal analysis techniques is thermometry. It involves the measurement of temperatures and the measurement of time. Magnetic Resonance Imaging (MRI)
has the ability to map temperatures (Cline et al., 1994), and it has been more than 30 years that several extensive studies have been performed to understand the quality of temperature monitoring in thermal treatment. The hotspot must be located correctly during ablation therapy with the use of MR guidance. It is necessary to locate the ablation site extremely precisely in order to burn just the unhealthy cells and spare the normal ones. Temperatures are achieved using microwave (MW) or radio frequency (RF) or ultrasound (US) or infrared (IR) techniques. Thermal therapy can be divided into two techniques. Low temperature or hyperthermia (HT), where tumour tissue is heated to a minimum temperature of 40 to 44\({}^{\circ}C\) for 60 minutes with the aim of directly killing cancer cells, increasing oxygenation, and thus also increasing radiosensitisation of the cancer cell (Kim and Hahn, 1979). Local, regional, and whole body hyperthermia can be classified on the basis of the size of the heated area. External heat sources, as well as intraluminal or interstitial insertion of microwave guided or wires, can be used to apply heat to the tumour. High-temperature thermal ablation, where the tumour tissues are heated to temperatures of 50-80\({}^{\circ}C\) or higher for a shorter period of time with the aim of killing cancer cells(Thomsen, 1991).

### Deep learning in medical imaging

The use of deep learning in the medical field, especially in the field of medical imaging, is increasing rapidly. Deep learning achieved outstanding performance in the task of undersampled MR image reconstruction or the elimination of artefacts present in those MRIs (Qin et al., 2018; Lyu et al., 2021). Wang et al. (Wang et al., 2016) applied deep learning to compressed sensing MRI. Deep Residual Network (ResNet) (He et al., 2016) was proposed in 2016 to optimise and improve the accuracy of deep learning architecture. ResNet was able to solve the vanishing gradient problem and open up the door to a deeper network. ResNet was proposed mainly for image classification, but later it has been used for many other applications, such as image classification (Mou et al., 2017; Zhang et al., 2019), image segmentation (Pakhomov et al., 2019) and image denoising (Jifara et al., 2019; Zhu et al., 2017). The residual learning model proved to be very efficient in MRI reconstruction as well (He et al., 2016).

One of the most commonly used network architectures for MRI reconstruction is UNet (Hyun et al., 2018), which was first employed for the task of MR reconstruction in 2018. UNet is capable of reconstructing highly undersampled images. In 2021 (Chatterjee et al., 2022a) Chatterjee et al. came up with the NCC1701 pipeline with the ReconResNet model as the backbone. This has shown an improvement in reconstruction of undersampled Cartesian and radial MRIs over UNet, and demonstrated that it is capable of reconstructing up to acceleration factors of 20 and 17, for Cartesian and radial MRIs, respectively.

UNet and ReconResNet work only in the image space (magnitude images), completely discarding the phase images. Even though these methods work with both image and k-space, they apply real-valued convolution operations on the complex-valued image space and k-space data, which leads to a destruction of the rich geometric relationship present within the complex data. In 2021, the first time complex-valued convolutions were applied for the task of undersampled MRI reconstruction directly in the k-space (Chatterjee et al., 2021).

On the other hand, Adler and Oktem 2018 (Adler and Oktem, 2018) proposed the primal-dual network or PDNetwork or PDNet for the reconstruction of sparse computed tomography (CT) data. Given that CT and radial MRI reconstructions have mathematical similarities due to the Fourier slice theorem, Ernst et al. (Ernst et al., 2023) in 2023 applied PDNet successfully for the task of undersampled radial MRI reconstruction, and also extended PDNet to PDUNet - which outperformed PDNet with statistical significance. Both networks employ two types of network blocks, filtering in the image space and singogram space. In 2022, these models were further extended using complex-valued convolutions into Fourier-PDNet and Fourier-PDUNet by Chatterjee et al. (Chatterjee et al., 2022c). These two models work in the k-space (i.e. the raw data space of MRI) instead of the sinogram space, along with working in the image space.

Although several models have been proposed for reconstructing undersampled MRIs, including models that work in both image and k-space, and models that work directly with complex data, the main focus of the evaluations carried out on those manuscripts was on the magnitude images. Current research aims to bridge this gap by focussing on both magnitude, as well as phase images, and then further evaluating the quality of MR thermometry from the reconstructed MRIs.

### MR-guided thermometry

With the benefit of obtaining 3D temperature maps, MR-guided hyperthermia provides a non-invasive approach for temperature monitoring(Wust et al., 2006).On the basis of proton density, T1 or T2 relaxation time, the water molecule's molecular diffusion coefficient, magnetisation transfer, temperature-sensitive contrast agents, proton resonance frequency (PRF) shift imaging, or spectroscopy, various methods of measuring temperature with an MR system have been reported(Kuroda, 2005; Ludemann et al., 2010; Quesson et al., 2000; Rieke and Pauly, 2008; Wlodarczyk et al., 1999). Techniques such as measuring longitudinal and transverse relaxation times (Parker et al., 1983), the diffusion coefficient, or the proton density rely heavily on the characteristics of the tissue. Because of this, the PRF shift technique is now the preferred technique for MRI-based temperature measurements due to its potential for online imaging and tumour control during treatments, ie, good linearity and high temperature sensitivity(Cernicanu et al., 2008; Gellermann et al., 2005). The PRF shift method's pre-clinical calibrations and uses are outlined in(McDannold, 2005). The current standard for noninvasive temperature assessments in daily clinical practice is the PRF shift measurement. The goal of the guided system is to understand the real-time temperature distribution and deliver quality controlled treatment and also be able to co-relate treatment temperature with treatment outcome in terms of actual thermal tissue damage. Fig.1 1 shows an example of MR-based temperature monitoring at different time points.
## 2 Methodology

Most of the previous work took either of these two directions: working only with the magnitude images (ignoring the phase images completely) or working with the complex image by splitting the data into real and imaginary parts before supplying it to the network as two separate channels. The first approach is not suitable for the current task at hand, while the second approach destroys the rich geometric structure present in the complex data. Both these approaches apply real-valued convolution operations. As the data are complex-valued, applying complex-valued convolution should be better suited, which is capable of working directly with the complex-valued data without splitting them into channels, effectively preserving the geometric structure.

### Experiment Design

MR images of 44 patients with different sarcoma cancers who have received the HT treatment in a combination of radiation/chemotherapy were used in this study (refer to Fig. 2). One key goal of this work is the reconstruction of the temperature, so phase images are as important as magnitude images. As the next step, the magnitude and phase images were combined together, and complex images were created. These complex images are artificially undersampled. Afterwards, the undersampled complex images were randomly divided into three different sets - Training, Validation and Test sets and the number of subjects was 26, 7 and 11, respectively. In the following steps, we have used a training and validation set to train our PD Net (Adler and Oktem, 2018) / PDUNet (Ernst et al., 2023) models, and a test set has been used to test the models. After testing, to validate the results produced by the models, the Structural Similarity Index (SSIM) (Renieblas et al., 2017)) has been used.

### Network Architectures

Primal-Dual network or PD Net: The Primal-Dual network is a deep learning-based technique for computed tomography data with sparse sampling (Adler and Oktem, 2018). The algorithm unrolls a proximal-dual method with convolutional neural networks in place of the proximal operators to accommodate for (potentially non-linear) forward operators in deep neural networks. The algorithm is trained end-to-end, using only the raw measured data, and is not dependent on any initial reconstruction, such as filtered backprojection. This not only raises the standard of the final reconstruction, but also ensures data consistency. The quality of PD Net depends on the number of iterations, just like many iterative algorithms, for coverage of all the parameters of the network, an optimal number of iterations is needed. The fewer the parameters of the convolutional block, the more iterations are needed for convergence.

Primal-Dual UNet or PD UNet: Primal-Dual UNet(Ernst et al., 2023) is the improved version of the Primal-Dual network in terms of accuracy and reconstruction speed. A UNet has been used in place of a convolutional block of PD net for image space to obtain a higher number of parameters with low processing time.

### Data consistency step

In the data consistency step, the actual acquired undersampled data replaces the network's output. The network only helps to fill the data which were ignored before during the undersampled data, this is how the final output is not totally dependent on the network.

Following (Hyun et al., 2018), a data consistency step was performed for the after reconstructing the undersampled Cartesian data. To obtain the corresponding k-space, FFT was performed on the output image first. Then, to identify the k-space values that were not acquired, an inverted mask was applied to this. The missing estimated k-space data from the network were combined with the measured data. To obtain the final output, iFFT was applied to this combined k-space.

### Dataset

In this work, MRI images of 44 patients treated at the Department of Radiation Oncology of the University Hospital Er

Figure 1: Example of non-invasive PRF shift MR thermometry to monitor and control temperature during clinical hyperthermia (HT). The initial temperature map at time step 0 requires two MR images, which are calculated by voxel-wise subtracting the second phase image from the first reference phase image. This reference image is subtracted from the phase images of further acquired MR images, which are taken every 10 minutes during HT therapy. In the initial high-precision magnitude MR image, the temperature map is shown as a colour overlay (blue: relative temperature decrease; green: constant temperature; red: relative temperature rise).
Figure 3: The network design known as the Primal-Dual Network (PDNet) was utilised to address the missing data problem in computed tomography. Primal iterates are displayed in blue boxes, whereas dual iterates are displayed in green boxes. The architecture of all the blue boxes is the same and is shown in the matching large boxes. When several arrows lead to the same box, concatenation occurs. As the data are transmitted to the dual iterates, the initial estimates enter from the left. Primal-Dual Network (PDNet): Network architecture used to solve the tomography problem. The dual iterates are shown in green boxes, whereas the primal iterates are shown in blue boxes. All blue boxes have the same architecture, which is illustrated in the corresponding large boxes. Several arrows pointing to one box indicate concatenation.

Figure 2: Experiment Design
langen, acquired from 2015 to 2020, who underwent HT treatment with MR thermometry, were used. Image sets of 44 different patients with different types of sarcoma cancer, originating mainly in the leg of the patients (details are in Table 3), have been acquired at 1.5T, the scanning sequence is GR or Gradient Recalled, Sequence Name fl2D or Fast low angle shot (FLASH 2D). The age range of the patients is from 23 to 80 years. 26 subjects out of 44 were used for the training set selected randomly; seven subjects were used for the validation set, and 11 subjects were used for testing the model. MR Thermometry images have two types of acquisition, static acquisition and dynamic acquisition, and the parameters of these two types of acquisition have been shown in Table 1 and in Table 2 respectively.
datasets do not contain any raw MR data, using the MR-Under (Chatterjee, 2020) pipeline, the single-channel fully sampled raw data and various undersampled datasets were generated artificially.

Cartesian raw data have been artificially undersampled using the k-space sampling pattern (sampling mask)(Lustig et al., 2007) which was created by randomly choosing completely sampled readout lines in the phase encoding direction, with the centre of the distribution following a one-dimensional normal distribution the distribution 5 (Fig. 5 (a)) that matches the k-centre space (referred to as 1D Varden). Furthermore, another mask is designed to include a densely sampled centre consisting of eight lines while gradually decreasing the sampling density toward the edges of k-space(Lustig et al., 2007)(referred as 2D Varden)5(Fig 5(b)) densely sampled centre covers 2. 5% of the k space randomly and distributed according to a two-dimensional normal distribution pattern. Three distinct Cartesian undersampling patterns were used in the first round of trials. 1D and 2D Varden masks were generated by randomly sampling 25% or 10% of the k-space - achieving acceleration factors of 4 and 10, respectively.

### Implementation

As input, the complex undersampled images are used to train the network to obtain the reconstructed complex images. However, most deep learning networks are implemented only on real-valued data, not complex-valued data. Therefore, the use of complex-valued convolution or CV-CNN(Chatterjee et al., 2022c) was a necessity. The convolution operation is the main component of CNNs, and it is computed by the sum of the product of two functions - the input (x) and the kernel (w) and the outcome is referred to as the feature map or activation map (s), and it is given by:

\[s(t)=(w\star x)(t)=\sum_{a}x(t+a)w(a) \tag{1}\]

In this case, w and x are both real-valued. Complex-valued convolutional networks, commonly referred to as CV-CNNs, improve on this by using the complex-valued convolution operation, which is defined as:

\[\begin{split} C_{r}(t)&=(w_{r}\star x_{r})\left(t \right)-(w_{i}\star x_{i})\left(t\right)\\ C_{i}(t)&=(w_{i}\star x_{r})\left(t\right)+(w_{r} \star x_{i})\left(t\right)\end{split} \tag{2}\]

where \(x_{r}\) and \(x_{i}\) are the real and imaginary components of the complex-valued input \(x\), respectively. Similarly, \(w_{r}\) and \(w_{i}\) are components of the complex-valued kernel \(w\), and \(C_{r}\) and \(C_{i}\) are components of the generated complex-valued feature map s. This can also be expressed in matrix notation:

\[\left[\begin{array}{c}\Re(\mathbf{w}\star\mathbf{x})\\ \Im(\mathbf{w}\star\mathbf{x})\end{array}\right]=\left[\begin{array}{cc} \mathbf{w}_{r}&-\mathbf{w}_{i}\\ \mathbf{w}_{i}&\mathbf{w}_{r}\end{array}\right]\star\left[\begin{array}{c} \mathbf{x}_{r}\\ \mathbf{x}_{i}\end{array}\right] \tag{3}\]

CV-CNNs can learn more sophisticated representations while preserving the algebraic structure of complex-valued data. Fig.6 shows the working mechanism of the proposed framework.

### Training and Inference

Fig.6 demonstrates the operational principles of the proposed framework, which include a network backbone called ReconResNet and for reconstructions using Cartesian undersampling. During the training process, only the network backbone was used. However, the entire framework is used during inference. L1 or mean absolute error (MAE) has been calculated to train ReconResNet.

The absolute difference between a prediction and the actual value is calculated for each sample in a dataset known as the L1 loss or Mean Absolute Error (MAE).

\[f(y,\hat{y})=\frac{1}{N}\sum_{i=1}^{N}|y_{i}-\hat{y}_{i}| \tag{4}\]

Where \(y\) is the actual value or ground truth, \(\hat{y}\) is the predicted value, and \(N\) is the whole dataset. The model was trained for 100 epochs, with a batch size of one, and the loss value was then minimised using Adam Optimiser (Initial learning rate 0.0001, decayed by 10 after every 50 epochs; \(\beta_{1}=0.9,\beta_{2}=0.999,\epsilon=1e-09\)). This network was implemented using PyTorch (Paszke et al., 2019), Python version 3.10.9 was used and was trained using NVIDIA GeForce RTX 2080 Ti.

### Evaluation Criteria

Structural Similarity Index 5(SSIM) (Renieblas et al., 2017), Normalised root-mean-squared error (NRMSE, 7) and Universal Image Quality Index (UIQI) (Wang and Bovik, 2002) have been used to evaluate the results.

The range of SSIM values is between zero and one, where the higher the SSIM value, the higher the similarity between two images.

\[SSIM(x,y)=\frac{(2\mu_{x}\mu_{y}+C_{1})(2\sigma_{x}y+C_{2})}{(\mu_{x}^{2}+\mu_ {y}^{2}+C_{1})(\sigma_{x}^{2}+\sigma_{y}^{2}+C_{2})} \tag{5}\]

where \(x\) and \(y\) are the two images between which the structural similarity is to be calculated, \(\mu_{x},\mu_{y},\sigma_{x},\sigma_{y}\) and \(\sigma_{xy}\), are the pixel means of \(x\), pixel means of \(y\), standard deviations, and cross-covariance for images x and y, respectively. \(c_{1}=(k_{1}L)^{2}\) and \(c_{2}=(k_{2}L)^{2}\) where \(L\) is the dynamic range of the pixel-values, \(k_{1}=0.01\) and \(k_{2}=0.03\).

Figure 5: (a) 1D Varden Mask and (b) 2D Varden Mask. All of them are for image size 256x256, taking 30thc k-space.
To statistically compare the two images (output and ground truth), NRMSE was used, calculated as:

\[MSE=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-\hat{Y}_{i})^{2} \tag{6}\]

\[NRMSE=\frac{\sqrt{MS\,E}}{\sqrt{\frac{1}{n}\sum_{i=1}^{n}Y^{2}}} \tag{7}\]

where the pixels of the fully sampled ground truth image have been denoted as \(Y_{i}\), the pixels of the undersampled image or the reconstruction (depending on the comparison performed) have been denoted as \(\hat{Y}_{i}\) and \(n\) denotes the number of pixels in the image.

Universal Image Quality Index (UIQI)(Wang and Bovik, 2002) Loss of correlation, luminance distortion, and contrast distortion are the three components that make up any image distortion when it is modelled. The proposed index is simple to calculate and adaptable to numerous image processing applications, as opposed to using conventional error summation techniques.

\[Q=\frac{\sigma_{x_{y}}}{\sigma_{x}\sigma_{y}}\cdot\frac{2\bar{x}\bar{y}}{(\bar {x})^{2}+(\bar{y})^{2}}\cdot\frac{2\sigma_{x}\sigma_{y}}{\sigma_{x}^{2}+ \sigma_{y}^{2}} \tag{8}\]

\[\bar{f}=\frac{1}{MN}\sum_{i=0}^{M-1}\sum_{j=0}^{N-1}x[i,j]\quad\bar{y}=\frac{1 }{MN}\sum_{i=0}^{M-1}\sum_{j=0}^{N-1}y[i,j]\]

\[\sigma_{xy}=\frac{1}{M+N-1}\sum_{i=0}^{M-1}\sum_{j=0}^{N-1}(x[i,j]-\bar{x})(y[ i,j]-\bar{y})\]

\[\sigma_{x}^{2}=\frac{1}{M+N-1}\sum_{i=0}^{M-1}\sum_{j=0}^{N-1}(x[i,j]-\bar{x})^ {2}\]

\[\sigma_{y}^{2}=\frac{1}{M+N-1}\sum_{i=0}^{M-1}\sum_{j=0}^{N-1}(y[i,j]-\bar{y}) ^{2}\]

where x and y are two images, considered as matrices having M and N number of columns and rows with x[i,j], y[i,j] pixels where (0\(\geq\)i \(>\) M, 0 \(\geq\)j \(>\) N ) and \(Q\) is the Universal image quality index.

\(Q\) can be obtained by multiplying three components together. The correlation coefficient is the initial component, which quantifies the level of linear correlation between the images x and y; the range varies [-1,1]. The second component assesses the similarity of mean luminance between images and has a value range of [0, 1]. With a range of [0, 1], the third component quantifies how closely the contrasts of the images match.

### Temperature Map

Temperature maps were generated to evaluate the retrieved temperatures by the Fourier-PDUNet and Fourier-PDNet models, as well as from the undersampled inputs and ground-truth images for comparison. To create these temperature maps, the proton Resonance Frequency Shift (PRFS) (McDannold, 2005) method has been used. For MRI-based temperature measurements, the PRF shift approach is currently the clinically recommended practice.

Linear proton Resonance Frequency Shift (PRFS) of water molecules with temperature is one of the widely used methods for MR-based temperature monitoring (McDannold, 2005). In Gradient Recalled Echo (GRE) images, the change of resonance frequency is expressed as phase change. The temperature difference can be derived by calculating the phase difference between a GRE image at a certain temperature and a reference temperature (Ishihara et al., 1995). The linear relationship between the temperature difference and the phase change can be expressed as the following equation.

\[\Delta T=\frac{\phi(T)-\phi(Tref)}{\Upsilon\alpha\beta TE} \tag{9}\]

where \(\Delta T\) = Temperature difference, \(\alpha\) = Temperature sensitivity of PRFS, \(\Upsilon\) = Gyro-magnetic constant, \(\delta\) = Main Magnetic

Figure 6: Working Model
field strength and TE = Echo time. According to Eq.10, a complex calculation has been performed to construct the phase difference, which could avoid the phase wrapping problem during the heating cycle (Peters, 2000).

\[\Delta\phi=\mathrm{atan}\left(\frac{\mathrm{Re}(I_{\mathrm{ref}})\cdot\mathrm{Im} (I_{H})-\mathrm{Im}(I_{\mathrm{ref}})\cdot\mathrm{Re}(I_{H})}{\mathrm{Re}(I_{ \mathrm{ref}})\cdot\mathrm{Re}(I_{H})+\mathrm{Im}(I_{\mathrm{ref}})\cdot \mathrm{Im}(I_{H})}\right) \tag{10}\]

where Re and \(I_{m}\) are the real and imaginary components of the heated (\(I_{H}\)) and reference (\(I_{ref}\)) images.

## 3 Results

In this research, images from 44 patients with sarcoma cancer were undersampled with an acceleration factor of 4, resulting in average SSIM values of 1D varden 25% is 63% and 31%, for magnitude and phase images where the Fourier-PDNet and Fourier-PDUNet models managed to reconstruct those data with average SSIM values of 91% and 90% for magnitude images, while achieving 44% and 40% for phase images, respectively. The results have been displayed using the violin plots in Fig. 10. Example outputs from two different subjects for 1D varden 25% sampling patterns are shown in Fig. 7 for qualitative evaluation.

The average SSIM values of the 2D varden undersampled MRIs with 25% undersampling were 43% and 29% for the magnitude and phase images. The Fourier-PDNet and Fourier-PDUNet models managed to reconstruct those with average SSIM values of 94% and 93% for the magnitude images while achieving 47% and 46% for the phase images, respectively. The SSIM values using violin plots and example results are shown in Figures 10 and 8, respectively.

Finally, with an acceleration factor of 10, the SSIM values of the undersampled images of the 2D varden 10% were 39% and 28% for the magnitude and phase images, respectively. Fourier-PDNet and Fourier-PDUNet models improved the SSIM values to 87% and 86% for the magnitude images while achieving 43% and 41% for the phase images, respectively. The violin plots of the SSIM values and the qualitative comparison for two subjects are shown in Figures 10 and 9, respectively.

Table 5 provides a complete qualitative overview of the results using NRMSE and UIQI, along with the SSIMs.

The temperature difference between the ground truth and the highest undersampled images (2D varden 10%) was 1.299\(\pm\)0.032, which is 1.3 \({}^{\circ}\)C more than the ground truth. But the models managed to reduce the difference to 0.618\(\pm\)0.016 and 0.643\(\pm\)0.022, using Fourier-PDNet and Fourier-PDUNet models, respectively, which are only around half a \({}^{\circ}\)C more than the ground truth 4. So, the models give 60% better accuracy in reconstructing the temperature maps compared to undersampled MRIs. This means that the model can speed up MR acquisition by a factor of 10 with only half \({}^{\circ}\)C of temperature difference. Examples of the reconstructed temperature map are shown in Fig. 11.

## 4 Discussion

The assessment of the proposed framework for reconstructing MRI images from undersampled data has shown that it can be applied effectively not only to MRI but also to MR thermometry images, as the model was also capable of reconstructing the temperatures. To the best of the knowledge of the authors, this manuscript is the first one to deal with undersampled MR thermometry, and in extension, MR-guided hyperthermia, using deep learning based methods; while this is also the first research discussing the need and possibility of accelerating MR-guided hyperthermia.

From the results, it can be observed that the framework seems to be robust against various undersampling patterns. For example, the SSIM values of 1D varden 25% are 63% and 31%, for the magnitude and phase images where the Fourier-PDNet and Fourier-PDUNet models managed to reconstruct those data with average SSIM values of 91% and 90% for the magnitude images, while achieving 44% and 40% for the phase images, respectively. SSIM values of the 2D varden 25% are 43% and 29% for the magnitude and phase images. The Fourier-PDUNet and Fourier-PDUNet models managed to reconstruct those data with average SSIM values of 94% and 93% for the magnitude images while achieving 47% and 46% for the phase images. 2D varden 10% is 39% and 28% for the magnitude and the phase images. The Fourier-PDNet and Fourier-PDUNet models managed to reconstruct those data with average SSIM values of 87% and 86% for the magnitude images while achieving 43% and 41% for the phase images, respectively, and the result has been displayed in violin plot10 as well as in table 5

The results show that both the Fourier Primal-Dual network (PD Net) and Fourier Primal-Dual UNet (PDUNet) were able to alleviate the undersampling problem and show that the deep learning model has the potential to improve the novel hyperthermia treatment. From the quantitative result Table 5, it has been clear that Fourier PD net outperformed Fourier PDUNet in SSIM. The same phenomenon has been observed in UIQI, but a slightly different phenomenon has been reported for NRMSE. For the 1D varden 25% and the 2D varden 25%, the value of NRMSE for PDUNet outperformed the output PDNet.

Also, from the result, it has been clear that both of the models performed way better for magnitude images than the phase images, which could be the reason for the present 0.5\({}^{\circ}\)C temperature difference in the reconstructed temperature from the ground truth. Improvement of the models for phase images can also decrease the difference in temperature.

Finally, this research demonstrates the possibility of accelerating MR thermometry during MR-guided hyperthermia. By using undersampling patterns such as 2D varden 10%, the acquisition can be ten times faster, and the methods presented here reduce the compromise in terms of temperature accuracy considerably. The faster acquisition would not only reduce the chances of having motion artefacts due to movements of the patient during the scan, but it would also significantly improve the temporal resolution of the imaging, making it possible to capture subtle changes in temperature over time.
\begin{table}
\begin{tabular}{|c|c|c|c|c|} \hline
**Type of Undersampling** & **Undersample** & **Fourier-PDNet** & **Fourier-PDUNet** \\ \hline
1D Varden 25\% & 1.221\(\pm\)0.035 & 0.637\(\pm\)0.016 & 0.641\(\pm\)0.018 \\ \hline
2D Varden 25\% & 1.296\(\pm\)0.031 & 0.596\(\pm\)0.014 & 0.657\(\pm\)0.019 \\ \hline
2D Varden 10\% & 1.299\(\pm\)0.032 & 0.618\(\pm\)0.016 & 0.643\(\pm\)0.022 \\ \hline \end{tabular}
\end{table}
Table 4: The mean value of MSE of reconstructed temperature map from different Models

\begin{table}
\begin{tabular}{|c|c|c|c|c|c|c|c|c|c|} \hline \multirow{3}{*}{**Input**} & & \multicolumn{3}{c|}{**2D Varden 10\%**} & \multicolumn{3}{c|}{**2D Varden 25\%**} & \multicolumn{3}{c|}{**1D Varden 25\%**} \\ \cline{3-10}  & & SSIM & NRMSE & UIQI & SSIM & NRMSE & UIQI & SSIM & NRMSE & UIQI \\ \hline \multirow{2}{*}{**Input**} & Undersample & _Magnitude_ & 0.385\(\pm\)0.009 & 0.541\(\pm\)0.022 & 0.536 \(\pm\)0.003 & 0.343\(\pm\)0.008 & 0.478\(\pm\)0.015 & 0.563\(\pm\)0.002 & 0.653\(\pm\)0.005 & 0.358\(\pm\)0.012 & 0.705\(\pm\)0.003 \\  & _Phase_ & 0.28\(\pm\)0.011 & 1.147\(\pm\)0.006 & 0.356 \(\pm\)0.013 & 0.294\(\pm\)0.011 & 1.128\(\pm\)0.005 & 0.363\(\pm\)0.012 & 0.312\(\pm\)0.009 & 0.172\(\pm\)0.005 & 0.382\(\pm\)0.011 \\ \hline \multirow{3}{*}{**Output**} & Fourier-PDUNet & _Magnitude_ & _0.886\(\pm\)0.004_ & _0.121\(\pm\)0.002_ & _0.807\(\pm\)0.006_ & _0.941\(\pm\)0.002_ & _0.675\(\pm\)0.001_ & _0.551\(\pm\)0.003_ & _0.509\(\pm\)0.001_ & _0.129\(\pm\)0.002_ & _0.843\(\pm\)0.004_ \\  & & _Phase_ & _0.429\(\pm\)0.01_ & _0.935\(\pm\)0.006_ & _0.433\(\pm\)0.01_ & _0.47\(\pm\)0.009_ & 0.983\(\pm\)0.004 & _0.501\(\pm\)0.009_ & _0.443\(\pm\)0.001_ & 0.961\(\pm\)0.005 & _0.477\(\pm\)0.01_ \\ \cline{2-10}  & Fourier-PDUNet & _Magnitude_ & 0.864\(\pm\)0.012 & 0.864\(\pm\)0.012 & 0.799\(\pm\)0.008 & 0.933\(\pm\)0.011 & 0.084\(\pm\)0.005 & 0.526\(\pm\)0.009 & 0.905\(\pm\)0.002 & 0.144\(\pm\)0.002 & 0.834\(\pm\)0.003 \\  & _Phase_ & 0.405\(\pm\)0.013 & 1.095\(\pm\)0.003 & 0.452\(\pm\)0.011 & 0.462\(\pm\)0.014 & **0.963\(\pm\)0.01** & 0.498\(\pm\)0.012 & 0.401\(\pm\)0.012 & _0.934\(\pm\)0.005_ & 0.466\(\pm\)0.011 \\ \hline \end{tabular}
\end{table}
Table 5: Result of Cartesian undersampling patterns, while being trained separately

Figure 8: Qualitative Result: 2D varden 25%

Figure 7: Qualitative Result: 1D varden 25%
Figure 10: (a)SSIM Value of reconstructed magnitude images of different models. (b) SSIM values of reconstructed phase images of different models. Where the blue violin diagram shows the SSIM value of Fourier-PDUNet, the orange violin plot shows the SSIM value of Fourier-PDNet, and the green violin plot shows the SSIM value of Undersampled images.

Figure 9: Qualitative Result: 2D varden 10%
## 5 Conclusion and future works

This paper introduced deep learning based reconstruction of undersampled MR thermometry acquired during hyperthermia using Fourier-PDNet and Fourier-PDUNet models. After all the different experiments with different types of undersampling methods of different percentages, the results show that the methods were able to alleviate the undersampling problem and managed to get SSIM Score of 0.886\(\pm\)0.004 for magnitude images and 0.429\(\pm\)0.01 for phase images for highest undersampling pattern 2D varden 10% - which means that the MR acquisition is now ten times faster and it also manages to bring the temperature difference close to the ground-truth which is 0.618\(\pm\)0.016. Still, half a \({}^{\circ}\)C temperature difference can be seen in the deep learning results. This can be attributed to the performance difference of the models between the magnitude and phase images.

Future work will focus on improving the networks' performance on the phase images, which should also reduce the temperature difference. Furthermore, combining the Fourier-PDNet and Fourier-PDUNet models with dynamic MRI-centric pipelines (Sarasaen et al., 2021; Chatterjee et al., 2022b) could allow these models to better exploit the spatio-temporal nature of MR thermometry data, improving the overall reconstruction quality.

Another future direction for research is to focus on the latent space. Exploring the latent space can be useful for improving the image reconstruction quality in undersampled image reconstruction tasks where the input images contain artefacts; the latent space represents, in theory, a low-dimensional representation of the input images without the artefacts. The input images with artefacts can be considered as augmented versions of the input images. Different types of variational auto-encoder(Makhzani et al., 2015) methods can be used, such as Factorised Variational Auto-encoder (FactorVAE) (Kim and Mnih, 2018), Vector Quantised Variational Auto-encoder (VQ-VAE) (Van Den Oord et al., 2017), Masked autoencoders (MAE)(He et al., 2022) etc. Use of post hoc explainability methods like Saliency (Simonyan et al., 2013), Occlusion (Zeiler and Fergus, 2014), Guided Backpropagation (Mahendran and Vedaldi, 2016) etc. can give a better understanding of what went wrong with phase images, which can help the authors to improve the network accordingly.

## Acknowledgement

This research has received support from the European Union's Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie (MSCA-ITN) grant "Hyperboost" project, no. 955625.
