Advances in Mathematical Physics

Volume 2015 (2015), Article ID 472818, 7 pages

http://dx.doi.org/10.1155/2015/472818

## Local System Matrix Compression for Efficient Reconstruction in Magnetic Particle Imaging

^{1}Section for Biomedical Imaging, University Medical Center Hamburg-Eppendorf, 20246 Hamburg, Germany^{2}University of Technology, 21073 Hamburg, Germany^{3}Bruker Biospin MRI, 76275 Ettlingen, Germany

Received 15 October 2014; Revised 11 December 2014; Accepted 11 December 2014

Academic Editor: Xiao-Jun Yang

Copyright © 2015 T. Knopp and A. Weber. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### 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 cm^{3}, 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 mm^{2}. 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].