3d-SMRnet: Achieving a new quality of MPI system matrix recovery by deep learning

TL;DR

This paper introduces a deep learning framework called 3d-SMRnet that significantly accelerates the recovery of 3D MPI system matrices, surpassing compressed sensing methods in quality and speed, and adapts to different particle types.

Contribution

The paper presents a novel 3D system matrix recovery network that achieves higher subsampling factors and faster recovery times than existing compressed sensing approaches.

Findings

Recovered 3D system matrix with a subsampling factor of 64 in under a minute.

Outperformed compressed sensing in matrix quality and image reconstruction.

Capable of inferring system matrices for various particle types.

Abstract

Magnetic particle imaging (MPI) data is commonly reconstructed using a system matrix acquired in a time-consuming calibration measurement. The calibration approach has the important advantage over model-based reconstruction that it takes the complex particle physics as well as system imperfections into account. This benefit comes for the cost that the system matrix needs to be re-calibrated whenever the scan parameters, particle types or even the particle environment (e.g. viscosity or temperature) changes. One route for reducing the calibration time is the sampling of the system matrix at a subset of the spatial positions of the intended field-of-view and employing system matrix recovery. Recent approaches used compressed sensing (CS) and achieved subsampling factors up to 28 that still allowed reconstructing MPI images of sufficient quality. In this work, we propose a novel framework with a 3d-System Matrix Recovery Network and demonstrate it to recover a 3d system matrix with a subsampling factor of 64 in less than one minute and to outperform CS in terms of system matrix quality, reconstructed image quality, and processing time. The advantage of our method is demonstrated by reconstructing open access MPI datasets. The model is further shown to be capable of inferring system matrices for different particle types.

Figures

Tables

Table 1. Table 1: Numerical results for all reconstructed phantoms demonstrates that our 3d-SMRnets clearly outperform CS when comparing the same subsampling factors.

	Shape Phantom			Resolution Phantom			Concentration Pht.			Avg.
	NRMSE	SSIM	PSNR	NRMSE	SSIM	PSNR	NRMSE	SSIM	PSNR	$\overline{\text{NRMSE}}$	$\overline{\text{SSIM}}$	$\overline{\text{PSNR}}$
CS8×	0.0275	0.9908	51.64	0.0120	0.9996	67.00	0.0098	0.9995	65.93	0.0164	0.9966	61.52
CS27×	0.0628	0.9394	44.47	0.0267	0.9975	60.06	0.0206	0.9972	59.42	0.0367	0.9780	54.60
CS64×	0.0915	0.8804	41.19	0.0452	0.9922	55.50	0.0349	0.9916	54.86	0.0572	0.9547	50.52
3d-SMRnet8×	0.0186	0.9959	55.03	0.0087	0.9998	69.79	0.0066	0.9998	69.39	0.0113	0.9985	64.74
3d-SMRnet27×	0.0320	0.9866	50.31	0.0208	0.9985	62.25	0.0135	0.9988	63.09	0.0221	0.9946	58.55
3d-SMRnet64×	0.0284	0.9874	51.36	0.0249	0.9978	60.68	0.0150	0.9986	62.22	0.0228	0.9946	58.09

Equations

$$\mathbf{S}=({\hat{s}}_k({\mathbf{r}}_n){)}_{k=0,\dots ,K-1;n=0,\dots ,N-1}\in {\mathbb{C}}^{K\times N}\text{vspace}-0.1cm$$

$$\underset{{\mathbf{s}}_k}{\min }{\Vert \mathbf{\Phi }{\mathbf{s}}_k\Vert }_1\text{subject to }\mathbf{P}{\mathbf{s}}_k={\mathbf{y}}_k.\text{vspace}-0.15cm$$

Full text

1 Introduction

Magnetic particle imaging (MPI) is a young tomographic imaging technique that quantitatively images magnetic nanoparticles with a high spatio-temporal resolution and is ideally suited for vascular and targeted imaging [6]. One common way to reconstruct MPI data is the system matrix (SM)-based reconstruction [7]. It requires a complex-valued SM, which is currently determined in a time-consuming calibration measurement. A delta sample is moved through the field-of-view (FOV) using a robot and the system response is measured in a calibration process. The number of voxels encoded in the SM directly determines the image size but also the scan time. The acquisition of a $37\times 37\times 37$ voxel SM takes about 32 hours, compared to an $9\times 9\times 9$ SM, which takes about 37 minutes. Therefore, a compromise between image size and scan time is usually made. While in principle the calibration of the SM needs to be done only once, the resulting SM is only valid for a very specific set of scan parameters. When changing scan parameters such as the size or the position of the FOV, the SM calibration needs to be redone. Furthermore, the SM highly depends on the type of the particles and their binding state, viscosity and even temperature. This makes it almost impossible to record high-resolution 32-hour system matrices for each combination of scan parameters and particle settings.

The first work that investigated calibration time reduction for MPI applied a compressed sensing (CS) to recover a subsampled 2d SM with 10-fold subsampling [8]. The basic idea is to exploit the fact that the MPI system matrix rows (i.e. frequency components) are consisting of wave-like patterns with a certain oscillation degree. By applying a discrete cosine transform (DCT), each matrix row can be sparsified such that CS with an $L_{1}$ prior can be applied. Motivated by the success, recent work proposed a CS approach by combining CS with a method that exploits symmetries in the SM. For up to 28-fold subsampling, sufficient image quality after reconstruction was archived [14].

While the CS approach for SM recovery is promising, it still leaves room for improvement since the sparsification using the DCT is not perfect and row specific. Furthermore, the CS approach currently cannot take prior knowledge from existing high resolution (HR) system matrix measurements into account. In the present work, we will, for the first time, investigate if deep learning (DL) techniques can be used to improve the SM recovery problem in MPI. DL-based super-resolution techniques have been demonstrated to be superior in the up-scaling of images in computer computer vision [4, 11] and recently also for 3d medical image up-scaling [12, 3, 2]. While Super-Resolution Convolutions Networks (SRCNNs) and Super-Resolution Generative Adversarial Networks (SRGANs) are mostly used directly in the image domain, it is inherently difficult to restore texture and structural details. SRGANs have proven to be successful in CV and medical image processing in modeling visually more appealing images than SRCNNs; SRCNNs typically tend to blur the SR image as a whole. Yet, this property may potentially be beneficial if SR is applied prior to image reconstruction – like in the current case for MPI SM recovery.

To evaluate the potential of SRCNN-based MPI SM recovery, we present a novel framework that comprises three central steps (see proposed method branch in Fig. 1). First, we acquire a low resolution (LR) SM on a specific sampling grid. Secondly, we encode each complex number of the 3d-SM to RGB vectors, allowing us to leverage SRCNNs from CV or medical image processing. Thirdly, we employ a SRCNN, which we call 3d System Matrix Recovery network (3d-SMRnet), to recover a high resolution (HR) SM by adapting the model to work on 3d RGB input data and employing it to each frequency component of the SM. Finally, we decode each RGB vector of the high resolution SM back to a complex number and use this newly recovered SM to reconstruct a high resolution image.

We evaluate our method in Sec. 3 on the Open MPI Data and will show that our framework reaches superior performance in image quality, SM quality, and processing time compared to the current state-of-the-art. All afore-mentioned aspects (introduction of 3d-SMRnet; conversion of MPI raw data to RGB format; comparison of 3d-SMRnet to compressed sensing) are novel contributions.

2 Methods

In MPI, the relation between the particle concentration $c(\cdot)$ and the Fourier coefficients $\hat{u}_{k}$ of the induced voltage $u(t)$ at frequency $f_{k}$ can be described by the linear integral equation $\hat{u}_{k}=\int_{\Omega}\hat{s}_{k}(\mathbf{r})c(\mathbf{r})\text{d}^{3}r$ where $\mathbf{r}$ is the spatial position, $k$ is the frequency index, and $\hat{s}_{k}(\cdot)$ denotes the system function. By sampling the FOV at $N$ positions $\mathbf{r}_{n},n=0,\dots,N-1$, one obtains a linear system of equations written in matrix-vector form $\mathbf{Sc}=\hat{\mathbf{u}}$ where $\hat{\mathbf{u}}=(\hat{u}_{k})_{k=0}^{K-1}\in\mathbb{C}^{K}$ and $\mathbf{c}=(c(\mathbf{r}_{n}))_{n=0}^{N-1}\in\mathbb{R}^{N}$ are the measurement vector and the particle concentration vector, respectively. $K$ is the total number of frequency components. The goal of this work is to recover a high resolution system matrix

[TABLE]

using a subset of the total number of sampling positions $N$. In our work, we treat the SM rows as independent images and train a network based on the entire set of rows of a measured HR SM. Then for SM recovery one measures a LR SM and infers the HR SM from the trained network. The method consists of three main steps, which are outlined in Fig. 1 and explained in detail below.

System Matrix Sampling: Usually, the system matrix is acquired using a 3-axis linear robot to scan the FOV equidistantly on a predefined grid. In this way, both HR and LR system matrix can be acquired. This is comparable to the generation of a LR SM by sampling every $\text{n}^{\text{th}}$ voxel from a HR SM. We employ the latter to get a LR and HR SM pair as training data for our 3d-SMRnet.

RGB Encoding and Decoding of System Matrix: Before we feed the frequency components of the SM into our 3d-SMRnet, we transform each complex number $\hat{s}_{k}(\mathbf{r}_{n})$ of the SM to a RGB color vector $s^{\text{RGB}}_{k}(\mathbf{r}_{n})\in\mathbb{R}^{3}$ in two steps. First, we use the hue-saturation-value (HSV) color model to represent the phase $\arg\hat{s}_{k}(\mathbf{r}_{n})$ with hue following the color wheel. Therefore, we employ the transformation $T_{\text{HSV}}:\mathbb{C}\to\mathbb{R}^{3}$ with $T_{\text{HSV}}(\hat{s}_{k}(\mathbf{r}_{n}))=(H,S,V)=(\arg\hat{s}_{k}(\mathbf{r}_{n}),1,1)$, where we omit S and V by setting them to 1. Secondly, we convert the HSV color vector $s^{\text{HSV}}_{k}(\mathbf{r}_{n})=T_{\text{HSV}}(\hat{s}_{k}(\mathbf{r}_{n}))$ to $s^{\text{RGB}}_{k}(\mathbf{r}_{n})$ with a standard HSV to RGB transformation [5]. Finally, the amplitude $|\hat{s}_{k}(\mathbf{r}_{n})|$ is coded by the intensity. Hence, we linearly scale the RGB color vector by the amplitude.

For decoding, we recover the complex numbers by extracting the scaling factor of the RGB color vector for the amplitude. Afterwards, we normalize the RGB color vector by the amplitude and convert it to a HSV color vector. The phase is now the hue value.

**3d-System Matrix Recovery Network: ** Following [13] and [9], we extend the SRCNN with Residual-In-Residual-Dense-Blocks (RRDBs) from 2d-RGB to 3d-RGB image processing. Our model contains two branches: image reconstruction and feature extraction. The feature extraction branch consists of $R$ stacked RRDBs (here $R=9$). Each RRDB combines three dense connected blocks and four residual connections as illustrated in Fig. 1. The dense connected blocks are built upon five convolutional layers. The image reconstruction branch generates the final up-scaled image by $U$ up-convolution blocks, followed by two convolutional layers. The up-convolution block contains a nearest-neighbor interpolation upsampling and a convolutional layer as proposed by [10] to reduce checkerboard artifacts from deconvolution. In our model, all 2d convolutions are replaced by 3d convolutions. Hence, the 3d-SMRnet recovers a HR SM by employing it to each frequency component $K$ of a LR SM.

3 Materials and Experiments

We apply our framework to the Open MPI Data111https://magneticparticleimaging.github.io/OpenMPIData.jl/latest/ dataset; It contains two HR system matrices, one for the particles Perimag $\mathbf{S}^{\text{Peri}}_{\text{HR}}$ and another for the particles Synomag-D $\mathbf{S}^{\text{Syno}}_{\text{HR}}$. Both are acquired using a $4\leavevmode\nobreak\ \text{}\mu$L delta sample with a concentration of 100 mmol/L and a grid size of $37\times 37\times 37$. Hence, $\mathbf{S}^{\text{Peri}}_{\text{HR}}$ and $\mathbf{S}^{\text{Syno}}_{\text{HR}}$ have the dimensions $37\times 37\times 37\times K$. Furthermore, three different phantom measurements with Perimag are provided in the Open MPI Data: Shape Phantom, Resolution Phantom, and Concentration Phantom. In our experiment, we train our 3d-SMRnet on Synomag-D with frequency components of the subsampled $\mathbf{S}^{\text{Syno}}_{\text{LR}}$ as input against $\mathbf{S}^{\text{Syno}}_{\text{HR}}$ and test it on the subsampled $\mathbf{S}^{\text{Peri}}_{\text{LR}}$. This represents the interesting case where the SM for new particles is inferred from a network trained on an established particle system. In addition, this approach prevents overfitting of the data.

We evaluate the recovered SM results in two steps. First, we compare all recovered system matrices with the ground truth $\mathbf{S}^{\text{Peri}}_{\text{HR}}$ by calculating the normalized root mean squared error (NRMSE) for each frequency component. Secondly, we reconstruct the measurements of the Open MPI shape, resolution and concentration phantoms with all recovered system matrices using the same standard regularization parameter ($\lambda=0.01,\text{iter}=3$). **Implementation Details and Training: ** We implement three versions 3d-SMRnet8×, 3d-SMRnet27×, and 3d-SMRnet64×. For 3d-SMRnet8× and 3d-SMRnet27×, we set $U=1$ and 2-times and 3-times upsampling, respectively. Furthermore, 3d-SMRnet64× uses $U=2$ and 2-times upsampling to finally upsample 4-times. To generate a LR and HR SM pair for training, we zero-pad $\mathbf{S}^{\text{Syno}}_{\text{HR}}$ to $40\times 40\times 40$ with two rows and one row at the beginning and the end, respectively. Afterwards, we apply 8-fold, 27-fold, and 64-fold subsampling, resulting in $20\times 20\times 20$, $13\times 13\times 13$ and $10\times 10\times 10$ spatial dimensions for the input volume, respectively. After applying a threshold with a signal-to-noise ratio (SNR) of 3, $\mathbf{S}^{\text{Peri}}_{\text{HR}}$ and $\mathbf{S}^{\text{Syno}}_{\text{HR}}$ have $K=3175$ and $K=3929$ frequency components, respectively. We split $K$ of $\mathbf{S}^{\text{Syno}}_{\text{HR}}$ into 90% training and 10% validation data.

We use random $90^{\circ}$ rotations and random flipping as data augmentation. In total, we train for $2\cdot 10^{5}$ iterations. Each iteration has a minibatch size of 20 for 3d-SMRnet8× and 64 for 3d-SMRnet27× and 3d-SMRnet64×. For optimization, we use ADAM with $\beta_{1}=0.9$, $\beta_{2}=0.999$, and without regularization. We trained with a learning rate of $10^{-5}$ for 3d-SMRnet8× and $10^{-4}$ for 3d-SMRnet27× and 3d-SMRnet64×. The learning rate is reduced by two every $4\cdot 10^{3}$ iteration. As loss function, we employ the mean squared error (MSE). Our models are implemented in PyTorch and trained on two Nvidia GTX 1080Ti GPUs. To support the reproduction of our results and further research, our framework and code are publicly available at https://github.com/Ivo-B/3dSMRnet.

**Comparison to State-of-the-Art: ** Compressed sensing exploits the fact that the SM becomes sparse when a discrete cosine transform (DCT) is applied to its rows. As shown in [8, 1], such signals can be recovered from an subsampled measurement by solving regularized least squares problems of the form

[TABLE]

Here, $\mathbf{y}_{k}\in\mathbb{C}^{M}$ contains the values of the $k^{\text{th}}$ SM row at the measured points and $\mathbf{P}\in\mathbb{C}^{M\times N}$ is the corresponding sampling operator. Moreover, $\mathbf{s}_{k}\in\mathbb{C}^{N}$ is the SM row to be recovered and $\bm{\Phi}\in\mathbb{C}^{N\times N}$ denotes the DCT-II.

Since CS requires an incoherent sampling, it cannot be applied directly to the regular sampled LR SMs used for our 3d-SMRnet. Instead, we use 3d Poisson disc patterns to subsample $\mathbf{S}^{\text{Peri}}_{HR}$ and obtain incoherent measurements with the same number of samples as used by the 3d-SMRnet. For every frequency component, we then normalize the measurement $\mathbf{y}_{k}$ and solve (2) using the Split Bregman method. The solver parameters are chosen manually such that the average NRMSE for all frequency components is minimized.

4 Results and Discussion

Fig. 2 (left) shows the NRMSE plot and Fig. 2 (right) shows visualizations of recovery for three frequencies. The results show that all our 3d-SMRnets can correct the noisy characteristics in $\mathbf{S}^{\text{Peri}}_{\text{HR}}$ due to smoothing characteristics of training with a MSE loss function, whereas CS cannot. For all three reduction factors, the 3d-SMRnet has a lower mean NRMSE than CS: 0.040 vs 0.044, 0.048 vs 0.051, and 0.048 vs 0.077 for 8-fold, 27-fold, and 64-fold subsampling, respectively. While CS cannot sufficiently recover the SM for 64-fold subsampling, our 3d-SMRnet64× still recovers the SM with a $37.66\%$ lower NRMSE. This low NRMSE is comparable to the results of CS with a 8-fold subsampling. Furthermore, our model 3d-SMRnet64× is over $42$ times faster and takes $\approx 23.3\text{ sec}$ compared to CS64× with $\approx 17\text{ min}$ for SM recovery. Still, some frequency components tend to have a high NRMSE (NRMSE $>0.11$) for our 3d-SMRnet64×, whereas for the 3d-SMRnet8× they do not. This problem can occur because of the equidistant subsampling and the symmetric patterns in the SM.

In Fig. 3, we show one representative slice ($Z=19$) of the reconstructed shape and resolution phantoms. For CS, all reconstructed phantoms show a ”checkerboard“ noise and an overestimation of the voxel intensity, which increases for 27- and 64-fold subsampling. Our proposed 3d-SMRnet results produce smoother reconstructed images, with voxel intensities better resembling the ground truth data. Yet, the shape phantom reconstructed with the S${}^{27\times}_{\text{3d-SMRnet}}$ shows some artifacts, while for S${}^{64\times}_{\text{3d-SMRnet}}$ those are not present. Still, the results for the resolution phantom with S${}^{27\times}_{\text{3d-SMRnet}}$ are visually better than S${}^{64\times}_{\text{3d-SMRnet}}$ (see Fig. 3 second row). Table 1. lists the subject-wise average structural similarity index (SSIM), peak signal to noise ratio (PSNR), and NRMSE. For all three phantoms, the overall best results archived our 3d-SMRnet8× with 0.0113, 0.9985, and 64.74 for $\overline{\text{NRMSE}}$, $\overline{\text{SSIM}}$, and $\overline{\text{PSNR}}$. Compared to the second best CS8×, this is an improvement by $31.1\%$ and $5.2\%$ for $\overline{\text{NRMSE}}$ and $\overline{\text{PSNR}}$. Furthermore, our 3d-SMRnet64× is on par with CS27× for the resolution phantom and considerably better for the shape and the concentration phantom.

5 Conclusion

We presented a novel method based on a 3d-System Matrix Recovery Network to significantly shorten calibration time in MPI. Our method can recover a highly subsampled system matrix: Using 64-times less samples compared to the original SM still allowed sufficient recovery of the SM. We have further shown that our method not only outperforms the current state-of-the-art in SM recovery quality (CS), but also in reconstructed image quality and processing time. Furthermore, our 3d-SMRnet can be applied to different types of particles after training. In the future, it is of interest to evaluate different kinds of sampling methods and multi-color MPI where two particles are simultaneously imaged.