Abstract
Magnetic particle imaging (MPI) is a tomographic imaging technique that allows the determination of the 3D spatial distribution of superparamagnetic iron oxide nanoparticles. Due to the complex dynamic nature of these nanoparticles, a time-consuming calibration measurement has to be performed prior to image reconstruction. During the calibration a small delta sample filled with the particle suspension is measured at all positions in the field of view where the particle distribution will be reconstructed. Recently, it has been shown that the calibration procedure can be significantly shortened by sampling the field of view only at few randomly chosen positions and applying compressed sensing to reconstruct the full MPI system matrix. The purpose of this work is to reduce the number of necessary calibration scans even further. To this end, we take into account symmetries of the MPI system matrix and combine this knowledge with the compressed sensing method. Experiments on 2D MPI data show that the combination of symmetry and compressed sensing allows reducing the number of calibration scans compared to the pure compressed sensing approach by a factor of about three.
1. Introduction
Magnetic particle imaging (MPI) is a promising method for imaging the spatial distribution of magnetic nanoparticles (MNP) [1, 2] in vivo [3–5]. When such contrast agents are administered into the blood pool, MPI is capable of revealing the blood flow [6] and potentially allows the detection of coronary artery diseases provided that the resolution and sensitivity will be preserved when scaling up the preclinical scanner [6] to a human size.
When using MNPs in clinically approved concentrations, the relation between the particle distribution and the voltage signals induced in the receive coils of an MPI system can be assumed to be linear and thus formulated as a linear system of equations. But due to the complex dynamics of the nanoparticles, it is challenging to accurately model their precise physical behavior, especially for dynamic 2D and 3D imaging sequences like the commonly applied Lissajous sequence [2, 7, 8]. Thus, in practice, the MPI system matrix is usually not explicitly known [9]. The common method to determine the system matrix is to perform a time-consuming calibration scan where the system response at all spatial positions is measured using a point-like nanoparticle sample. This calibration measurement takes about 6 hours for volume grids of size covering a field of view (FOV) of mm3 [6]. When changing any scanning parameter such as the drive-field amplitude or the selection field gradient, the calibration has to be redone. Also, a change of the particle batch usually leads to small but noticeable changes in the MPI signal [10].
The first attempt to reduce the amount of calibration scans has been carried out in [11]. By exploiting the special structure of the MPI system matrix and applying compressed sensing techniques [12], the amount of calibration data could be reduced to about 20% of the entire system matrix without significantly affecting the image quality.
The purpose of the present work is to reduce the necessary amount of calibration data for recovering the MPI system matrix with sufficient accuracy even more. To this end, it is exploited that the MPI system matrix has symmetries in every spatial direction which allows taking the calibration scans only from one quarter and obtaining the full matrix by mirroring [13]. But in this way for 2D MPI one can only reduce the amount of calibration data to 25%. In order to achieve reduction rates that are below that of the pure compressed sensing approach, we therefore combine the symmetry method with the compressed sensing technique developed in [11].
2. Preliminaries
In MPI, the relation between the particle concentration and the th Fourier coefficient of the measurement signal can be described by the integral equationwhere is the spatial position in the and denotes the system function. By discretization of the FOV with sampling points , , a linear system of equationsis obtained, where is the measurement vector, is the particle concentration vector, andis the MPI system matrix. Here, one row of the MPI system matrix describes the spatial distribution of the corresponding frequency component and is commonly denoted as system-function component.
Normally, the MPI system matrix is obtained by a calibration procedure where a delta sample with the size of one voxel filled with the used tracer material is shifted to all discretization positions within the FOV and at each position the system response is measured. Each calibration scan corresponds to one column of the MPI system matrix which can be mathematically described by multiplying the system matrix with a unit vector. Here, the index of the nonzero component of the unit vector corresponds to the calibration position.
3. Symmetries of the System Function
With the assumption of ideal coils and isotropic particles that instantaneously follow magnetic field changes, symmetry properties of the system function can be established. To this end, the system-function component has to be separated bywhere is the transfer function of the receive chain and is denoted by signal function component. Due to the symmetric behavior of the Lissajous excitation sequence, the signal function exhibits symmetry properties as well, which depend on the frequency ratio of the excitation frequencies. In the setting of the 2D MPI experiment outlined in Section 6 the signal function can be mirrored in vertical direction by multiplication with a symmetry factor. In horizontal direction, additional to the multiplication with a symmetry factor, a complex conjugation has to be performed. A more detailed explanation with a theoretical derivation of the symmetry properties is given in [13].
In Figure 1, original system-function components and their mirrored versions are visualized. As one can see, the system-function components are highly symmetric. The differences at the border are due to drift artifacts caused by the long measurement time of the calibration procedure.
One challenge when utilizing the system matrix symmetries in practice is that the symmetry axis is not exactly known (e.g., due to field inaccuracies and shifts induced by the earth magnetic field). When the symmetry axis does not perfectly match the central position between two pixels, mirroring thus can lead to discontinuities which influence the reconstruction result negatively. To avoid this problem, a slightly larger region than one-half/one-quarter has to be measured. Using this overlap a smooth transition by a weighted averaging depending on the distance to the symmetry axes can be enforced.
In [13] it has been shown that the high similarity between the mirrored and the original system-function components can be utilized to reduce the number of calibration scans by a factor of three to obtain the same reconstruction results as when using the fully sampled system matrix.
4. System-Function Recovery Using Compressed Sensing
A further method to accelerate the calibration process is the usage of compressed sensing [12] to recover the system function from an undersampled set of calibration scans [11].
Basically, compressed sensing is a technique to solve an underdetermined linear system of equations, which has an infinite number of solutions. To obtain the desired solution, further constraints have to be imposed. In compressed sensing the usually imposed constraint is that the solution of the linear system is sparse; that is, it has only a few coefficients which are unequal to zero.
For the MPI system matrix it has first been observed in [14] that the individual system-function components are highly compressible when a certain basis transformation is applied, like the discrete Fourier transformation (DFT) or the discrete Cosine transformation (DCT). The basis transformation forms the linear system of equationsinvolving the discrete system-function component . The coefficient vector has only few nonzero components when is a DFT or DCT matrix.
When subsampling the spatial grid at which calibration scans are performed, the system-function components are subsampled as well. In turn, linear system (5) becomes underdeterminedwith a reduced measurement vector and a reduced basis transformation matrix . This underdetermined linear system of equations can be solved by compressed sensing. In order to calculate the complete system-function component , the basis transformation has to be applied to the solution after compressed sensing reconstruction.
Besides the requirement of a sparse solution , compressed sensing further requires that the reduced system matrix exhibits a low coherence. This is usually achieved by choosing the subsampling positions in a random fashion (see [15]).
There exists a multitude of algorithms to solve the compressed sensing problem [16–18]. In the present work, the fast iterative shrinkage-threshold algorithm (FISTA) is chosen [19] that has already proven its good performance for the MPI system matrix recovery problem [11].
In [11] it has been shown that at a sampling rate of 20% no significant difference in the reconstructed image quality can be observed. Even a reduction to 10% of the original calibration scans yielded an acceptable image quality. However, a further reduction to 5% degraded the reconstructed image quality considerably.
In Figure 2 several system-function components obtained with the pure compressed sensing method are visualized. Thereby compressed sensing was applied with the DCT and in combination with the FISTA algorithm. It can be clearly seen that the quality of the system-function components degrades with decreasing sampling rate and depends on the spatial structure of the system-function component since finer structures can be less compressed.
5. Combining Symmetries with Compressed Sensing
In the last section two independent ways to reduce the number of calibration scans necessary to reconstruct the MPI system matrix were discussed. Next, it will be discussed how the compressed sensing method can be combined with the symmetry method. For simplicity only vertical symmetries are considered in the illustrations. As it turns out, there are two different possibilities for combining symmetry properties with compressed sensing.
On the one hand, one can sample the system function in one-half and first apply compressed sensing in order to reconstruct the full half of the system-function component. In the second step, the entire system-function component is reconstructed by mirroring. This method, which we call CS/Symmetry, is illustrated in Figure 3. As the system-function components are nonperiodic when considering only one-half (or one-quarter) of the FOV we have found that the DCT performs better than the DFT when used as sparsity transformation for compressed sensing. This is due to the fact that the DFT assumes periodic data.
On the other hand, one can first mirror the sparse sampling positions such that they cover the entire FOV. In the second step, compressed sensing is applied in order to recover the dense system-function component. An illustration of this Symmetry/CS method can be found in Figure 4. Testing the DCT and the DFT, no significant difference in the performance up to a sampling rate of 5% could be recognized for the Symmetry/CS method. At smaller sampling rates the DFT shows better reconstruction results. Thereby the restriction of the sampling scheme has a larger influence on the DCT than on the DFT.
Due to the uncertainty of the real mirror axes (see Section 3), the CS/Symmetry method requires reconstructing the system-function component with a small overlap so that a smooth transition in the mirroring process can be enforced. In contrast, this overlap can be neglected in the Symmetry/CS method, because the compressed sensing reconstruction does inherently smooth the data such that no discontinuities occur.
6. Materials and Methods
In the present work we use the same MPI data as in [11]. The data has been measured with the MPI mouse scanner that was first published in [6]. It has a bore size of about 32 mm and allows imaging small mice. The selection field gradient strength is 5.5 in -direction and 2.75 in -direction. The drive amplitudes are 18 . A 2D Lissajous imaging sequence has been performed with drive-field frequencies MHz/96 ≈ 26041.7 kHz and MHz/99 ≈ 25252.5 kHz. The resulting repetition time considering an averaging factor of 17 is ms.
While the signal is sampled with 20 MS/s, only frequency components up to 1 MHz are used for reconstruction. In the frequency range above 1 MHz no detectable MPI signal was measured. The total number of frequency components in each receive channel is 1268. As only two of three receive channels are required for 2D imaging, the total number of system matrix rows is 2536.
In the presented pictures of system-function components and reconstruction results the horizontal direction coincides with the -direction and the vertical direction with the -direction.
For reference a fully resolved system matrix is captured at equidistant positions spanning a FOV of . The different system matrix recovery methods are evaluated by using only a subset of the positions of the full system matrix. For the combined Symmetry/CS and CS/Symmetry approach, the random samples are taken in one-half (vertical and horizontal symmetry) or one-quarter (when exploiting both symmetries) considering an additional overlap for the CS/Symmetry method. The sampling rates reported in the results section are always relative to the total number of positions of 2516. For the combined Symmetry/CS and CS/Symmetry method, the sampling rate takes both the symmetry undersampling and the CS undersampling into account.
In order to determine the impact of the system matrix error on real MPI images, measurements of a rotating P phantom consisting of 12 holes each of a diameter of 0.5 mm were acquired.
For image reconstruction, the iterative Kaczmarz method is used considering 5 iterations and a regularization parameter, which is chosen such that the best visual result is obtained (see [4]).
7. Results
7.1. System Matrix
First we consider the reconstruction of the system matrix using the different recovery methods and calculate the deviation from the fully sampled system matrix.
In the first experiment, we test whether it is better to first mirror the system function and then apply compressed sensing or whether the reverse order is better. In Figure 5 the mean normalized root mean squared error (NRMSE) considering the 100 system matrix rows with highest signal-to-noise ratio (SNR) is shown for different sampling rates. In this example a vertical mirroring has been applied. As one can see, for sampling rates above 20% both methods perform equally. However, for sampling rates below 20% one can see that the Symmetry/CS method is better than the CS/Symmetry method. In the following we therefore consider only this procedure.
It should be noted that the mean NRMSE shown in Figure 5 is limited by the symmetry error for sampling rates above 20%. The symmetry error is about 5.5% for the considered dataset and vertical symmetry. For sampling rates below 20% the compressed sensing methods contribute an additional error that increases for decreasing sampling rates.
Next, it is checked if the combination of mirroring with compressed sensing works equally well for the different mirror combinations. The NRMSE are shown in Figure 6. As one can see, for each mirror combination (vertical, horizontal, and both) there is a different baseline error that is independent from the sampling rate. It is at about 5.5% for vertical mirroring, 4.4% for horizontal mirroring, and 6.2% when exploiting both mirror axes. Although the procedure that uses both mirror axes has a higher symmetry error it behaves better for low sampling rates (i.e., below 7% sampling rate). Thus, to achieve the highest undersampling one should use both mirroring axes and apply compressed sensing.
To demonstrate that exploiting the system-function symmetries provide a benefit compared to the pure compressed sensing approach the NRMSE of the reconstructed system matrices are compared in Figure 7. One can see that the pure compressed sensing method performs better at modest sampling rates because it has no systematic symmetry error due to drift artifacts. However, below 10% sampling rate the combined approach clearly outperforms the pure compressed sensing approach.
Several system-function components which are obtained with the Symmetry/CS method are visualized in Figure 8. At a sampling rate of 7.5% no significant difference to the pure mirrored system-function components can be seen. A further decreasing of the sampling rate causes artefacts whereby as for the pure compressed sensing method the system-function components with the finest spatial structures exhibit the largest degradation.
7.2. Phantom Measurements
After analyzing the NRMSE of the MPI system matrix for the different reconstruction techniques, how the error translates to the MPI phantom images when using the approximated system matrices for image reconstruction is investigated next.
As a reference, two frames of the captured P phantom time series reconstructed with the original system matrix are shown in Figure 9 (first column). It can be seen that the dots of the P phantom can only be resolved in vertical direction. This is due to the fact that the MPI scanner has a higher gradient strength and in turn resolution in the vertical axis of the scanner/images.
Besides the reference images, Figure 9 also shows the reconstruction results that are obtained when using the mirrored system matrices for reconstruction. As one can see, only minor differences can be observed compared to the images reconstructed with the full system matrix. Here, there is no change in the result if the horizontal, the vertical, or both mirror axes are exploited.
In Figure 10, the phantom images reconstructed with the compressed sensing recovered system matrices are shown. As one can see, the compressed sensing approach is not restricted to sampling rates of 25% like the symmetry method. Instead, the number of positions to be used for system matrix recovery can be reduced to 10% with still acceptable image quality compared to the reference images. However, below 10% sampling rate, the image quality degrades considerably such that the dots cannot be resolved anymore for 7.5% sampling rate.
Finally, Figure 11 shows the phantom images reconstructed with system matrices obtained with the combined Symmetry/CS method. One can see clearly that the combined Symmetry/CS method outperforms the pure symmetry and the pure CS approach. Even at 3.7% sampling rate the dots of the P can be clearly resolved. When comparing the results of the pure compressed sensing with the combined Symmetry/CS method, one can see that the image quality and resolution of the 10% pure compressed sensing method roughly match those of the 3.7% combined method. Hence, the combined Symmetry/CS method allows reducing the number of sampling points by a factor of about three.
8. Discussion
It has been shown that the combined Symmetry/CS method achieves much better results than the pure CS method for high sampling rates while at lower sampling rates the symmetry prevents that more accurate system functions are reconstructed using the combined Symmetry/CS method. For the experimental data considered in this work, the symmetry error was low enough that it could not be recognized in the reconstructed phantom images.
While a robot was employed to move the delta sample at the various sampling positions to calibrate the system function in the present work, an alternative method has been proposed in [20], where the delta sample stays at a static position and the robot movement is emulated by applying dedicated offset fields. Since the time to move the robot is commonly higher than switching the offset fields, an acceleration of the calibration can be achieved. Nevertheless the determination of a high resolution system matrix is still very time-consuming even when using an electromagnetic calibration procedure. Fortunately, the system matrix reconstruction method proposed in this work can be applied regardless of whether a mechanical or an electromagnetical calibration is performed in order to measure the columns of the system matrix. Hence, the findings of this work can be directly transformed to field-based system matrix calibration where calibration times in the sub-second range can be expected.
The experiments in this work have been carried out on 2D measurement data. While, in principle, the findings can be directly carried over to 3D imaging, the symmetry properties differ between the 2D and the 3D case. In the 2D case the system-function components are symmetric with respect to the vertical and horizontal axes. In contrast for the 3D case there exist settings depending on the ratios of the excitation frequencies where the 3D system-function components exhibit a mirror symmetry in a single direction and in addition a symmetry with respect to the center point in the planes perpendicular to this direction [13]. Hence, instead of a factor of eight that could have been expected, the feasible reduction factor when exploiting the symmetries of a 3D MPI system matrix is four as in the 2D case. When considering the potential reduction that can be achieved in 3D by compressed sensing, one can expect even higher rates than for 2D system matrices. This is because the 3D system function is sparse in all three dimensions when applying a 3D sparsity transformation. Thus, a further overall reduction of the necessary calibration scans can be expected for 3D imaging when applying the combined Symmetry/CS method.
The compressed sensing reconstruction of the system matrix is naturally related to the sparse image reconstruction technique introduced in [14], which allows accelerating image reconstruction in MPI. Both methods use a basis transformation to represent the system matrix in a compact form and their combination has been already shown in [21]. The sparse reconstruction technique can also be combined with the system matrix recovery method developed in this work, which will allow directly performing reconstruction without arranging the full system matrix at any time.
9. Conclusion
In the present work, it has been shown that for the reconstruction of the MPI system matrix the compressed sensing reconstruction technique can be combined with the symmetry method to achieve a further reduction of the required calibration scans. Using 2D experimental data, the feasibility to reduce the number of calibration scans by a factor of 30 compared to a full system matrix calibration has been shown without noticeable loss in resolution. Compared to the pure compressed sensing method, the reduction rate could be increased by about a factor of three using the combined Symmetry/CS method. At a reduction rate of 2.4%, the image quality degrades, but the rough outline of the measured nanoparticle phantom could still be imaged.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
The authors would like to thank Jürgen Weizenecker, Bernhard Gleich, Jürgen Rahmer, and Jörn Borgert (Philips Technologie GmbH Forschungslaboratorien, Hamburg, Germany) for providing them with measurement data and for numerous fruitful discussions on magnetic particle imaging in general.