Abstract
Magnetic particle imaging (MPI) is a quantitative method for determining the spatial distribution of magnetic nanoparticles, which can be used as tracers for cardiovascular imaging. For reconstructing a spatial map of the particle distribution, the system matrix describing the magnetic particle imaging equation has to be known. Due to the complex dynamic behavior of the magnetic particles, the system matrix is commonly measured in a calibration procedure. In order to speed up the reconstruction process, recently, a matrix compression technique has been proposed that makes use of a basis transformation in order to compress the MPI system matrix. By thresholding the resulting matrix and storing the remaining entries in compressed row storage format, only a fraction of the data has to be processed when reconstructing the particle distribution. In the present work, it is shown that the image quality of the algorithm can be considerably improved by using a local threshold for each matrix row instead of a global threshold for the entire system matrix.
1. Introduction
The quantitative imaging method magnetic particle imaging (MPI) allows resolving the spatial distribution of superparamagnetic iron oxide (SPIO) nanoparticles at high spatiotemporal resolution and high sensitivity [1–3]. Introduced in 2005, MPI experienced a rapid development leading to first preclinical results in 2009 revealing structures of the cardiovascular system of a living mouse [4].
While the MPI hardware is fast enough to acquire more than 40 frames per second for scanning volumes of about cm3, only little effort has been investigated into the question of how to reconstruct the particle distribution with similar temporal performance. One major challenging part during image reconstruction is to handle the memory requirements of the system matrix of the linear system that has to be solved. Due to an insufficient understanding of the particle physics, the system matrix is usually measured in a calibration procedure and stored on permanent memory for further processing. For a 3D MPI experiment with an image size of , the size of the system matrix is about 40 GB in double precision. Although advanced workstations can be equipped today with 64 GB and more, processing such amounts of memory can be a computational bottleneck. For instance, loading a 40 GB matrix from the permanent into the main memory takes about 7 minutes assuming a disk reading rate of 100 MB/s. When increasing the image size and in turn the density of the sampling trajectory, the size of the system matrix increases quadratically such that for image voxels more than 3 TB of main memory are required. While in current MPI reconstruction methods [5], only less than half of the full system is actually used and in turn is loaded into main memory, with better scanner hardware, the noise in the measured system matrix will be reduced so that more matrix rows have sufficient SNR and can be considered for reconstruction.
Recently, a matrix compression technique has been proposed, which allows for performing in-memory reconstruction even for large image sizes [6]. To this end, a matrix compression technique is used, which transforms the rows of the system matrix into a certain basis, yielding a sparse coefficient matrix. This matrix can then be compressed by applying a threshold and storing only the significant components. Such matrix compression techniques are well known in the literature, where commonly a wavelet transformation is used [7].
In [6], a global threshold has been used for the reduction of the system matrix. Without verification, it has been assumed that global thresholding is better than local thresholding. The contribution of the present paper is to show that the opposite is true. By applying a local threshold to each matrix row individually, the compression rate can be considerably increased while maintaining the image quality.
Besides discussing the impact of the thresholding strategy, we will show in this work that modest thresholding can even improve the image quality compared to a regular reconstruction without matrix compression. This is based on the noise reduction that results from truncation of coefficients which contain only noise.
2. Preliminaries
Magnetic particle imaging applies various static and oscillating magnetic fields for spatial encoding. Using a set of usually receive coils, the change of the particles’ magnetization progression is measured during the time interval , where is the repetition time of the experiment. The induced voltages are given by , . As the induced voltages are in practice covered by excitation signals directly coupling into the receive coils, the measurement signals are usually considered in frequency space, where the excitation signal can be removed by a narrow band-stop filter. The remaining frequency components contain only the signal induced by the magnetic nanoparticles. Due to finite sampling of the time signal, only discrete Fourier coefficients are available in practice. To store the frequency components of all receive channels in a 1D vector, one can introduce the mapping , if . The total number of frequency components is thus given by .
3. Imaging Equation
The aim of MPI is to reconstruct the particle concentration at discrete positions , , which can lay on a 1D, 2D, or 3D sampling grid. The relation between the particle concentration and the measured frequency component can be expressed by a linear system of equations where describes the contribution of particles located at position to the measured frequency components . In continuous form, is called the MPI system function. In matrix-vector notation, (2) can be expressed as where is the measurement vector, is the particle concentration vector, and is the system matrix. The rows of the system matrix are denoted by , , which we assume to be column vectors in this work.
4. Matrix Compression
As it has been discussed in [8, 9], the individual rows of the system matrix consist of wave-like functions bearing high similarity to tensor products of Chebyshev polynomials. But as the relationship is only qualitative, one cannot exploit this structure of the system matrix directly. Instead, as it has been proposed by [6], one can use matrix compression techniques to reduce the size of the system matrix, which accelerates the multiplications with considerably. The basic idea has already been used for general matrices in the context of the wavelet transformation.
In the discrete setting, we consider a set of linear independent basis vectors , , which set up a transformation matrix . In order to expand the th matrix row into the basis vectors, one has to solve the linear system where are the basis coefficients. If the transformation matrix is orthogonal, the coefficients can be computed by where . When considering all matrix rows , , (5) can be written as where . By introducing a coefficient vector satisfying the original linear system (3) can be reformulated as Instead of solving (3) directly, one thus can perform the following steps.(1)Perform the basis transformation .(2)Solve the linear system .(3)Compute the solution . While on a first sight these three steps seem to require more computational effort than solving (3), they actually can be carried out much more efficiently, which is substantiated by the following facts. First of all, the basis transformation of the system matrix has to be performed only once. In fact, this operation can be directly performed after measuring the system matrix. Furthermore, although a matrix multiplication has to be performed, for the mainly considered transformations, this multiplication can be performed in an efficient manner using fast algorithms (e.g., the fast Fourier transformation (FFT) and the fast cosine transformation (FCT)). Consequently, also the postprocessing step can be performed very efficiently by applying a fast basis transformation.
The main reason for expanding the system matrix rows into basis functions is that the second step, the solution of the linear system , can be considerably accelerated by exploiting the sparse structure of the transformed matrix . Hence, by applying a threshold and taking only few entries of into account, arithmetic operations involving the system matrix can be considerably accelerated. This is essential for iterative reconstruction algorithms, which, for instance, apply matrix-vector multiplications with the system matrix and its adjoint. Even more important is that the full MPI system matrix will not fit into the main memory of the reconstruction computer for practical voxel resolutions. Hence, matrix compression can enable in-memory reconstruction, which would not be possible when considering the full representation of the system matrix.
5. Selection of Matrix Rows
At this point, it should be noted that even for the conventional MPI reconstruction involving the solution of the linear system (3) usually only parts of the system matrix are considered. More precisely, the energy of the system matrix rows and in turn the energy of the measurements considerably vary in dependence of the time frequency. As is usually measured in a calibration procedure, there are several matrix rows that only contain noise. Therefore, as it has been discussed in [12], for reconstruction of MPI data, only the matrix rows with sufficient SNR are taken into account. The set of row indices to be used for reconstruction can be obtained by calculating the SNR of in dependence of frequency and applying a threshold.
Considering the sparse system matrix , one can apply the same row thresholding as for the dense system matrix . Matrix rows containing only noise will not be compressed by the basis transformation and in turn these rows should not be taken into account.
6. Matrix Thresholding
Besides selecting the threshold that determines the set of used matrix rows, one has to choose for each matrix row the number of coefficients (i.e., columns) that are used for reconstruction. In this work, we propose to use a fixed number of coefficients per matrix row. This is contrary to the procedure proposed by [6], where a global threshold was applied to the matrix . A global threshold implies that more coefficients of matrix rows with high energy are used than those of matrix rows with low energy. But as the details (i.e., the high-resolution part) of the particle distribution are encoded in matrix rows having low energy [10], using a global threshold leads to a loss of spatial resolution, which is experimentally proven in Section 8.
Note that for a fair comparison between local thresholding and global thresholding we apply the threshold only to those matrix rows that are selected by the row selection algorithm described in Section 5.
7. Materials and Methods
Experiments have been carried out using the MPI scanner published in [4], which has a bore diameter of 32 mm and is capable of imaging mice. The selection field has a gradient strength of 5.5 Tm−1 while the drive fields have an amplitude of 18 mT . A 2D Lissajous-type sampling sequence is performed using the two drive-field frequencies MHz/96 26041.7 kHz and MHz/99 25252.5 kHz. These correspond to a commensurable frequency ratio of 33/32, determining the density of the sampling trajectory. The resulting drive-field cycle is ms. By averaging 17 sequential frames, the SNR of the measured data is improved, leading to a total acquisition time of 21.5 ms per frame. The acquired MPI signals are expanded into Fourier series while considering frequencies up to 1 MHz. Hence, the total number of frequency components per receive channel is 1268. After combining the frequency components of the two receive channels, the total number of entries of the measurement vector is 2536.
The system matrix is acquired using a calibration measurement using a cubic delta sample with an edge length of 0.6 mm filled with undiluted Resovist (0.5 mol (Fe) ). The sample is shifted to positions on a regular grid of size covering a FOV of mm2. Regular measurements are acquired using a phantom consisting of 12 dots forming the letter P. The phantom is rotated by about 45° in a dynamic sequence (see Figure 1). The data has also been used in [9, 12].
For the compression of the MPI system matrix, we use the discrete cosine transform (DCT) that already yielded excellent results in [6]. Only matrix rows with a signal-to-noise ratio larger than 10 are used for reconstruction leading to a total number of 1400 rows.
Different thresholds are applied to the remaining matrix such that the amount of stored data is between 0.2% and 5% of the total number of matrix entries. Besides using a fixed number of coefficients per row, we compared this approach to a global threshold for reduction of the system matrix, which has been proposed in [6].
For image reconstruction, the iterative Kaczmarz method [13] is used considering 5 iterations and a regularization parameter, which is chosen such that the best visual result is obtained (see [12]).
All algorithms are implemented using the Julia programming language [14, 15], which allows us to formulate algorithms at a high abstraction level (similar to MATLAB) but run the programs in native C-speed.
8. Results
As a reference for the proposed sparse reconstruction framework using local and global thresholding, we use the common MPI reconstruction method, which uses the dense representation of the system matrix. The reconstruction results of the phantom data are shown in Figure 1. As one can see, the vertical resolution of the images is better than the horizontal direction. This is due to the gradient strength of the selection field being 5.5 Tm−1 in vertical direction and 2.75 Tm−1 in horizontal direction.
In Figure 2, the results of the sparse reconstruction algorithm outlined in Section 6 are shown for the DCT basis transformation and local as well as global thresholding. When considering 5% of the full system matrix for reconstruction, the results using the local and the global threshold show equivalent image quality. Interestingly, compared to the reference images using the full system matrix, the images are even slightly less noisy while retaining the spatial resolution (see also Figure 3). This can be explained by the fact that the removed basis coefficients mainly contain noise. Removing these coefficients thus can potentially even improve the image quality.
(a)
(b)
When reducing the number of used basis coefficients further, the image quality degrades as the removed coefficients then not only contain noise but also describe the shape of the MPI system function. Using a global threshold, the image quality degrades when reducing the number of coefficients to 1% of the total number of coefficients so that the dots in vertical direction can hardly be resolved when using the DCT as a basis transformation. In contrast, for the local threshold, the dots in vertical direction can even be resolved when only 0.5% of the total number of basis coefficients is used.
9. Discussion
In the present work, we have investigated the sparse reconstruction framework for magnetic particle imaging that allows the processing of large MPI datasets for which the full system matrix would not fit into the main memory of the reconstruction computer.
We have shown that one can achieve better image quality when thresholding the basis transformed MPI system matrix locally on each matrix row instead of globally on the complete matrix. This can be explained by the fact that the energy of the MPI system matrix rows varies considerably. When applying a global threshold, a large number of rows fall below the threshold while from individual rows with high energy more basis coefficients are taken than actually needed for an accurate approximation. Using the proposed local thresholding strategy, from each row, the same number of coefficients is taken which makes the method invariant of the row energy. In turn, even the low energy rows of the MPI system function carrying the high spatial frequency information of the particle distribution are present in the approximated system matrix. In order to show this effect, in Figure 4 one particular matrix row is shown for different compression rates and global as well as local thresholding. At 0.2% compression rate, the entire matrix row falls under the threshold when using the global thresholding strategy, while the local thresholding strategy still reveals the wave pattern of the original matrix row.
(a)
(b)
It should be noted that one reason that the local threshold performs better than the global threshold is that the SNR of the receiver is not constant over frequency. Typically the MPI-receive chain (filter and analog receive amplifier) is implemented in such a way that frequencies that have a lower MPI signal are amplified more than those where the MPI signal is stronger. This allows one to detect more frequency components than one would capture when a constant amplification over frequency would be used.
For the acquisition of larger scanning volumes, MPI will rely on a multipatch method where small volumes are acquired using a fast 3D drive-field sequence and the small measuring field is shifted in space over larger distances using the so-called focus field [11]. Due to field inhomogeneity, the system matrices for different focus-field positions will differ slightly such that for optimal image quality no single system matrix can be used for all focus-field patches. In this case, the matrix compression technique with local thresholding discussed in this paper can be applied to each system matrix individually, which will considerably improve the reconstruction performance for multipatch focus-field data.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.