Magnetic particle imaging visualizes the spatial distribution of superparamagnetic nanoparticles. Because of its key features of excellent sensitivity, high temporal and spatial resolution and biocompatibility of the tracer material it can be used in multiple medical imaging applications. The common reconstruction technique for Lissajous-type trajectories uses a system matrix that has to be previously acquired in a time-consuming calibration scan, leading to long downtimes of the scanning device. In this work, the system matrix is determined by a hybrid approach. Using the hybrid system matrix for reconstruction, the calibration downtime of the scanning device can be neglected. Furthermore, the signal to noise ratio of the hybrid system matrix is much higher, since the size of the required nanoparticle sample can be chosen independently of the desired voxel size. As the signal to noise ratio influences the reconstruction process, the resulting images have better resolution and are less affected by artefacts. Additionally, a new approach is introduced to address the background signal in image reconstruction. The common technique of subtraction of the background signal is replaced by extending the system matrix with an entry that represents the background. It is shown that this approach reduces artefacts in the reconstructed images.
Original content from this work may be used under the terms of the Creative Commons Attribution 3.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.
Magnetic particle imaging (MPI) is an emerging medical imaging technique that visualizes the distribution of superparamagnetic iron oxide nanoparticles (Gleich and Weizenecker 2005). It is a highly sensitive method that is capable of acquiring multiple three-dimensional (3D) frames per second with a good spatial resolution of about 1–4 mm (Panagiotopoulos et al 2015) and the potential for submillimetre resolution using tailored nanoparticles (Ferguson et al 2015). The nanoparticles are excited by oscillating magnetic fields, referred to as drive fields. As the nanoparticles have a non-linear magnetisation behaviour, the receive signal is a distorted version of the excitation signal. The fundamental frequency and its higher harmonics can be detected and the introduction of a magnetic gradient field enables spatial encoding. This gradient field is called the selection field. As different types of nanoparticles have unique magnetization curves, the spectral response of the receive signal is characteristic for the different particles.
The receive signal has to be reconstructed into an image in order to visualize the distribution of the particles. There are two types of reconstruction method, which are equivalent in theory (Grüttner et al 2013) but differ in practice. First, reconstruction can be performed in x-space (Goodwill and Conolly 2010). The temporal receive signal is corrected for the applied trajectory and mapped back into the spatial domain. For 1D excitation (used when realizing a Cartesian sampling scheme) the x-space approach leads to very good reconstruction results. However, the reconstruction of MPI data has not been shown to be feasible using x-space reconstruction techniques when acquisition is along Lissajous-type trajectories (Goodwill and Conolly 2011).
Alternatively, the particle distribution can be reconstructed in the frequency domain using a system matrix (Rahmer et al 2009). A system matrix encodes the particle response and all the imperfections of the system at previously discretized spatial positions of the field of view (FOV). As the spatially dependent receive signal can be described by a linear system of equations (Knopp et al 2010b), it can be reconstructed using the system matrix. The particle behaviour is encoded in the system matrix together with the gradient field distortions and the sensitivity of the receive coils. As the system matrix is ill-conditioned, the reconstruction needs to be regularized. Although automatic determination of the regularization factor has already shown promising results (Gräfe et al 2016), it is still often determined manually. A second challenge with the use of a system matrix is its acquisition.
The system matrix is most often measured using the scanning device itself. After the FOV has been discretized, a nanoparticle sample is shifted by a robot from each spatial position defined by the grid inside the FOV to the next one (see figure 1). At each position, a full measurement is performed and stored as a column vector in the system matrix. This procedure is very time-consuming as it is limited by the robot's movement. Optimistically, assuming an acquisition time of 1.3 s per pixel, a 3D FOV of voxels could be acquired in about 12 h without averaging the receive signal. When a large number of averages are applied the measurement can take up to several days. Furthermore, the signal quality in terms of the signal to noise ratio (SNR) is limited. The SNR decreases with a larger discretization of the FOV, because the sample should not be larger than the voxel in order to guarantee the spatial resolution of the system matrix.
In a clinical environment a scanning device needs to be available for treatment and diagnosis most of the time, with short maintenance times. The acquisition of a system matrix can be considered as maintenance time as it calibrates the reconstruction procedure of the system. Compressed sensing for MPI (Candes et al 2006, Weber and Knopp 2013, von Gladiss et al 2015) has been researched in order to reduce the number of voxels that need to be measured for acquiring a system matrix in a shorter calibration time. Knopp et al have shown that it could be sufficient to acquire about 10% of the voxels and then reconstruct the system matrix based on these undersampled data. Although this reduces the acquisition time of a single system matrix, there would still be a long downtime of a clinical MPI scanning device as system matrices need to be acquired for every type of acquisition parameter and nanoparticle. Furthermore, system matrices for nanoparticles of the same type but different batches may differ.
An alternative to measuring the system matrix is its simulation (Knopp et al 2010b). This features both an arbitrary SNR and a neglectable acquisition time, respectively simulation time. However, a precise model of the particle behaviour in multidimensional magnetic fields has not yet been formulated.
A second alternative to measuring the system matrix was introduced in 2011 by Graeser et al (2011). In their hybrid approach, the particle behaviour is measured in real magnetic fields. The static magnetic gradient field of the scanning device is assumed or measured previously and set as a consecutively changing offset field to the measurements in a separate device. A similar approach has been researched by Halkola et al (2012). By applying focus fields on the scanning device, different spatial positions of the nanoparticle sample have been emulated. A dedicated set of receive coils has been used in order to increase the SNR of the system matrix (Halkola et al 2013).
In this work, a magnetic particle spectrometer (MPS) is used to acquire a hybrid system matrix. As shown in figure 1, the acquisition of a hybrid system matrix with a MPS differs from acquisition inside the scanning device. In a scanning device, the spatial position of the nanoparticle sample is changed by robot movement. For the MPS, the nanoparticle sample stays at the same spatial position inside the sample chamber. Different spatial positions are emulated inside the MPS by changing the magnetic offset field that is applied to the particles. Second, one defined voxel of the FOV in a MPI scanning device covers an area of the gradient field, which is due to the spatial dependence of the magnetic field. Therefore, the particle sample experiences different field strengths in different directions within a single voxel. As a consequence, there are multiple non-uniform field sequences in one single voxel and the receive signal is a sum of particle responses to different field sequences. When acquiring a hybrid system matrix with a MPS the whole sample chamber represents one infinitesimal volume of the FOV in a MPI scanning device. The application of homogeneous offset fields generates the field sequence of this volume for all the particles inside the sample chamber. The SNR of the measurement inside the MPS can be adjusted by the volume of the nanoparticle sample without reducing the spatial resolution of the hybrid system matrix. It has been shown by Gruettner et al (2011) that 1D MPI data can successfully be reconstructed using a hybrid system matrix.
Recently, a novel device has been introduced by Graeser et al (2015) and Graeser et al (2016) that enhances the spectrometry for magnetic particles. The trajectory analysis particle spectrometer (TAPAS) is able to apply multidimensional excitation to the sample chamber. Furthermore, an offset field can be applied in order to emulate the magnetic gradient field. In this work, a hybrid 2D system matrix that has been measured with TAPAS will be compared with one that has been acquired using a pre-clinical MPI Scanner (Bruker BioSpin GmbH, Ettlingen, Germany; subsequently MPI scanner). First, the measurement procedures using TAPAS and the MPI scanner will be discussed. Then, the hybrid system matrix will be compared with the measured one in terms of SNR and measurement time. In order to reconstruct measurement data that have been acquired with the MPI scanner, a transfer function has to be applied to the hybrid system matrix (Knopp et al 2010b). A measured dataset will be reconstructed using both a hybrid matrix and the device's own system matrix. A novel approach for handling the background signal inside the reconstruction procedure will be introduced. Until now, the background has been corrected for by subtracting the background signal from the measurement signal. In this work, the system matrix will be extended by an additional entry representing the background signal. Using this approach, a changing background signal can be handled, as the reconstruction allows for scaling of the background. The reconstruction results will be analysed in terms of spatial resolution, SNR and reconstruction artefacts.
2.1. System matrix acquisition
When reconstructing a MPI dataset with a system matrix, the reconstruction problem is formulated as with the spatial particle distribution , the measurement in dependence on the frequency and the system matrix. In order to acquire the system matrix, the FOV is discretized into voxels. A sample of nanoparticles is placed into the center of one voxel by, for example, a robot (see figure 1). One measurement is performed and the receive signal is stored in the system matrix at the column that corresponds to voxel . Then the sample is moved to the next voxel. Following this procedure, measurements have to be performed in order to acquire a complete system matrix.
Although a single measurement itself may be taken in several milliseconds without averaging, the acquisition time per voxel is about a second due to movement of the robot. The acquisition of a full 3D system matrix could take several hours up to days, depending on the number of averages; this is unacceptable in a clinical scenario.
A second drawback of this type of system matrix acquisition is the dependence of the SNR on the voxel size. The sample size should be the maximum possible in order to achieve a high SNR, but it has to fit into one voxel in order to ensure precise allocation of the sample. Therefore, the sample size should correspond to the size of a voxel, which limits either the spatial resolution or the SNR of the system matrix.
2.2. Hybrid system matrix acquisition
The TAPAS device approaches these issues in an alternative way. TAPAS features a very small distance between the measurement chamber and the receive coils, which guarantees a high SNR. The static magnetic field of a scanning device is emulated by applying variable direct currents (DC) to sending coils (see figure 1). This has two advantages. (1) The whole measurement chamber represents one infinitesimal volume in the scanning device. Therefore, both a high SNR and arbitrary discretization of the emulated FOV can be ensured. The dependence between SNR and discretization of the FOV is resolved. (2) Switching the DC currents can be realized very quickly. As there is no robot movement, the average measurement time per voxel drops from seconds to milliseconds. Up to 400 measurements per second can be performed using TAPAS.
2.3. System matrix dataset
First, a set of system matrices was acquired with the MPI scanner. Both the excitation field amplitudes of the x- and y-direction where set to 14 mT with a gradient of 1.12 T m−1 in the x and y directions. A mm sized FOV was discretized by pixels, resulting in a pixel size of mm. As the FOV is centred around the field-free point (FFP), its boundaries lie within in both directions, x and y. One microlitre of undiluted Resovist was been used for the nanoparticles. The receive signal was averaged times to generate system matrices of different SNRs (see figure 2).
A second set of system matrices was acquired in the TAPAS device. The excitation field amplitudes of TAPAS were also set to 14 mT. The static magnetic field of the MPI scanner was assumed to be linear. Therefore, the offset field was adapted in equidistant steps to emulate the respective spatial positions inside the scanning device. Two different grids of and pixels were measured, each with extensions of in both x and y directions. Ten microlitres of undiluted Resovist (of the same batch as used with the MPI scanner) was been used as the nanoparticle sample. Three different system matrices with averages of were acquired (see figure 2 and 3).
Furthermore, the system matrices with a discretization of pixels were interpolated linearly to match the discretization of pixels (see figure 3). This allows us to differentiate the influence of grid resolution and SNR on the image reconstruction result.
2.4. Phantom data
A resolution phantom was constructed consisting of two cylinders with an inner diameter of 0.7 mm and a variable distance from each other. Beginning with a distance of 15 mm between the centres of the dots, the distance was reduced in steps of 1 mm down to 3 mm and a measurement was performed at each step. The receive signal was averaged 50 000 times. As the phantom is extracted from and reinserted into the scanner at each step, the position of the dots varies in the reconstructed images. The excitation field amplitude was set to 14 mT in the xy-plane to match the system matrix measurements.
2.5. Transfer function
As described in Knopp et al (2010b), a transfer function must be calculated by a least squares approach that represents the different receive path parameters of the systems. The transfer function of the TAPAS system to the magnetic moment has already been determined and its data have been corrected for. The missing transfer function from the magnetic moment to the MPI scanner's receive signal has been calculated by the acquired system matrices with a grid size of pixels and highest SNR. The result is visualized in figure 4.
The data measured with the MPI scanner are corrected with the calculated transfer function. After applying the transfer function on the data from the MPI scanner, the spectra match the measured data from TAPAS for low harmonics and mixing frequencies.
The spatial distribution of the nanoparticles is reconstructed by solving the linear system of equations using Kaczmarz's algorithm (Kaczmarz 1937), where may either be the system matrix acquired with the MPI scanner or the hybrid system matrix acquired with TAPAS. As the system matrix is usually ill-conditioned, the system is regularized (Knopp et al 2010a). The parameter λ indicates the regularization and needs to be found for each measurement dataset. Furthermore, the number of iterations of Kaczmarz's algorithm can be adjusted in order to reach the best reconstruction result. Table 1 shows the reconstruction parameters that were chosen for the different datasets.
Table 1. Parameters that have been used for the reconstruction. One set of parameters has been determined for each type of system matrix and used for reconstructing all the phantom measurements. The same set of parameters has been used for every system matrix that only differs in the number of averages.
2.6.1. Background correction.
In order to cope with the background of the scanning device, an empty measurement may be taken and subtracted from both the system matrix and the measurement data. In this work, it is suggested that the system matrix is extended by an empty measurement . This enables the reconstruction algorithm to reconstruct the background and noise data of the measurement into an additional voxel that is represented by the empty measurement. Furthermore, the empty measurement may be weighted by a parameter . Introducing allows the background to be scaled in order to cope with, for example, a temperature-driven changing background signal. This reconstruction will be referred to as background extended system matrix reconstruction (BESM).
3.1. System matrix measurement
The magnitude of single frequency components of the measured system matrices is visualized in figures 2, 3 and 6. Apparently, the signal quality of the system matrices acquired with TAPAS exceeds that of the system matrices measured with the MPI scanner. Even when averaging 10 000 times, the system functions of the MPI scanner show artefacts. The artefacts get worse with increasing frequency index and less averaging (see figure 2). By visual inspection of the system functions acquired with TAPAS, no difference can be observed between using 100 averages and not averaging at all.
The acquisition times of the system matrices were between 8 and 25 min for a FOV discretized by pixels dependent on the number of averages for TAPAS. As for the MPI scanner, the acquisition times of the pixels system matrices have been as short as 14 min and as long as 85 min depending on the number of averages. The system matrices with higher discretization ( pixels) were acquired in 12 min without averaging and 30 min when averaging 100 times by TAPAS. For measurement of a pixel system matrix with the MPI scanner the acquisition time would be about 122 min as a minimum without averaging. Comparing the measurement time per SNR, TAPAS achieves much quicker measurements, as a low averaged measurement with TAPAS still has a higher SNR than a highly averaged measurement with the MPI scanner. Therefore, the desired result can be obtained much faster using TAPAS.
The acquisition time of TAPAS could even be decreased further. For now, one of the magnetic offset field paths may be swept, which enables a line scan mode when acquiring a hybrid system matrix. Sweeping the second offset field path would enable a quadrant scan mode. The theoretical limit for acquiring a 2D system matrix of pixels is about 1.5 min without averaging.
Figure 5 shows the reconstruction results using system matrix of the MPI scanner and the hybrid system matrix of TAPAS. The two cylinders of the phantom were between 3 mm and 6 mm apart. In all of the images, the cylinders that were 5 mm and 6 mm apart can be separated visually. The spatial resolution of the system matrix is insufficient to separately visualize the cylinders that are 3 mm and 4 mm apart.
The left column has been reconstructed using the standard reconstruction. The background has not been corrected for because subtracting an empty measurement from the measurement signal and the system matrix did not improve the reconstruction result. The reconstruction results of the centre column were achieved using the system matrix that was extended by the empty measurement as a pixel. It can be seen that the reconstructed images have fewer artefacts using the background signal pixel that is cut off before visualizing. As the new approach for reconstructing the particle distribution has achieved better reconstruction results, it will be used further on.
The single cylinders of the phantom can best be distinguished when using the hybrid system matrix (right column) for reconstruction.
3.2.1. Discretization of the FOV.
The system functions of figure 2 have been interpolated linearly to a grid of pixels (see figure 3). Although the interpolated system functions shown in figure 3 are similar to the system functions of the higher-discretized grid, differences can be observed for higher frequencies, i.e. higher mixing orders (see figure 6). With increasing mixing order, the structures inside a system function become more tenuous, and are not sufficiently sampled by a grid of low discretization. The linear interpolation of the undersampled grid cannot represent the structures that are acquired by a system function with higher discretization (see figure 6). When the magnetic gradient field increases, the structures of the system function become even smaller (Rahmer et al 2012). Therefore, the FOV needs to be highly discretized when using a field with a high magnetic gradient, as a high discretization of the FOV corresponds to a high spatial resolution of the system matrix and the reconstructed spatial particle distribution.
Figure 7 shows the reconstruction of the same data used in figure 5 using the system matrix of the higher grid of pixels. As the spatial grid discretization of the system matrix has increased, so has that of the reconstructed images. The two cylinders of the phantom can clearly be separated for all the inter-cylinder distances.
The acquisition of a system matrix on an highly discretized grid is superior to linear interpolation of a system matrix to a higher grid. Although figure 3 promises high resolution of the linearly interpolated system matrices, figure 7 shows the limitation of this approach: for small inter-cylinder distances (3 mm), the reconstructed shape of the cylinders is distorted. Furthermore, strong artefacts appear all around the reconstructed cylinders.
The system matrices of both the MPI scanner and TAPAS were acquired with different numbers of averages. Increasing the number of averages, the SNR of a measurement increases; this has a direct influence on the image quality (Knopp et al 2011). Figure 8 compares the reconstruction results when the number of averages of the system matrix is varied. The number of averages has a strong influence when using the system matrices of the MPI scanner. With increasing number of averages the shape of the reconstructed cylinders becomes sharper. This effect is virtually negligible when using hybrid system matrices having the same spatial resolution, as the SNR is already very high without averaging. However, there is a positive effect regarding the background noise of the images when increasing to 10 averages and regarding the shape of the reconstructed cylinders when increasing to 100 averages using a hybrid system matrix of high spatial resolution.
As the MPI scanner has a low SNR in comparison with TAPAS, increasing the number of averages has a big effect on the measurement data. TAPAS has a high native SNR that exceeds that of the MPI scanner. Increasing the number of averages has a small positive effect on the reconstruction result as long as the SNR of the phantom measurement (i.e. taken with the MPI scanner) is lower than the SNR of the system matrix.
In this work, a multidimensional hybrid system matrix was acquired with the TAPAS device that matched a system matrix of an MPI scanner.
A phantom dataset was measured with the MPI scanner and reconstructed with both the device's own system matrix and the hybrid system matrix. It has been shown that multidimensional datasets measured on an MPI scanning device can be reconstructed successfully with hybrid system matrices. Comparing the results with those using the device's own system matrix, there are no drawbacks. Instead, the reconstructed images using the hybrid system matrix have a higher image quality, due to the higher SNR of the hybrid system matrices.
It is important to acquire system matrices of very high SNR when performing, for example, multicolor MPI (Murase et al 2015, Rahmer et al 2015, von Gladiss et al 2016). As the influence of single particle samples needs to be distinguishable, the noise level of the system matrices should be minimal. Using system matrices of high SNR, the quality of the reconstructed images is mainly influenced by the SNR of the single measurements.
The reconstruction algorithm has been extended in order to cope with the background signal. With addition of a background voxel to the system matrix, the reconstruction algorithm allows for reconstructing background data into the voxel that is discarded later for visualization.
As TAPAS resolves the dependence of SNR on the FOV discretization, another hybrid system matrix has been acquired with a much higher discretization of the FOV. As the spatial resolution of a reconstructed image is limited by the spatial resolution of the system matrix, the reconstruction results using the hybrid system matrix with a higher discretization display smaller structures. The single elements of the resolution phantom could be separated even at the smallest distance from each other.
The hybrid approach to system matrix determination has been shown to improve image quality in MPI. Currently, the use of hybrid system matrices for other scanner topologies, for example a single-sided MPI scanning device (Gräfe et al 2016), is being investigated. As the magnetic gradient field of other scanning devices may be inhomogeneous it might be necessary to cope with the crooked magnetization vector of the excitation field.
In this work, the magnetic gradient field of the MPI scanner has been assumed to be linear and the magnetic offset field has been emulated accordingly. Measuring the actual magnetic gradient field and adjusting the emulated fields to this measurement may improve the reconstruction results. This measurement has to be done once per scanning device, as the distortion of the magnetic gradient field is not supposed to change after the construction the scanning device.
Without loss of generality, this work has focused on two-dimensional datasets. In order to reconstruct three-dimensional datasets, the hybrid system matrix scanning device would need to provide a third dimension. The theoretical lower acquisition time limit for acquiring a hybrid system matrix of voxels and voxels is about 10 min and less than 5 h without averaging the receive signal, respectively. In comparison with this, the MPI scanner would need about 6 h when acquiring a system matrix of voxels without averaging and about 7 days when discretizing the FOV into voxels. However, it has been shown that a high number of averages would be needed to acquire a system matrix with the MPI scanner.
In a clinical scenario, multidimensional hybrid system matrices will be of great use as the scanning devices do not need to acquire system matrices themselves. The calibration downtime drops to an initial measurement of the magnetic gradient field. Then, hybrid system matrices can be acquired for every set of particles.
It has been shown that multidimensional MPI measurements can be reconstructed successfully using hybrid system matrices. As hybrid system matrices allow for a much higher discretization of the FOV in a reasonable amount of time, the spatial resolution of MPI images increases significantly. Furthermore, the acquired hybrid system matrices have a higher SNR, which is favourable for reduction of noise in the reconstructed images.
Taken together, hybrid system matrices allow for higher SNR and higher spatial resolution of the reconstructed images and reduce the calibration downtime of the scanning devices to a negligible amount of time. The use of hybrid system matrices is a huge step towards the realization of clinical MPI scanning devices.
A von Gladiss, M Graeser and T M Buzug gratefully acknowledge the financial support from the German Research Foundation (DFG, grant number BU 1436/10-1) and the Federal Ministry of Education and Research (BMBF, grant number 13GW0069A). T Knopp and P Szwargulski gratefully acknowledge the financial support by the German Research Foundation (DFG, grant numbers KN 1108/2-1 and AD 125/5-1).