Hybrid x-space: a new approach for MPI reconstruction

Magnetic particle imaging (MPI) is a new medical imaging technique capable of recovering the distribution of superparamagnetic particles from their measured induced signals. In literature there are two main MPI reconstruction techniques: measurement-based (MB) and x-space (XS). The MB method is expensive because it requires a long calibration procedure as well as a reconstruction phase that can be numerically costly. On the other side, the XS method is simpler than MB but the exact knowledge of the field free point (FFP) motion is essential for its implementation. Our simulation work focuses on the implementation of a new approach for MPI reconstruction: it is called hybrid x-space (HXS), representing a combination of the previous methods. Specifically, our approach is based on XS reconstruction because it requires the knowledge of the FFP position and velocity at each time instant. The difference with respect to the original XS formulation is how the FFP velocity is computed: we estimate it from the experimental measurements of the calibration scans, typical of the MB approach. Moreover, a compressive sensing technique is applied in order to reduce the calibration time, setting a fewer number of sampling positions. Simulations highlight that HXS and XS methods give similar results. Furthermore, an appropriate use of compressive sensing is crucial for obtaining a good balance between time reduction and reconstructed image quality. Our proposal is suitable for open geometry configurations of human size devices, where incidental factors could make the currents, the fields and the FFP trajectory irregular.

Magnetic particle imaging (MPI) is a new medical imaging technique capable of recovering the distribution of superparamagnetic particles from their measured induced signals. In literature there are two main MPI reconstruction techniques: measurement-based (MB) and x-space (XS). The MB method is expensive because it requires a long calibration procedure as well as a reconstruction phase that can be numerically costly. On the other side, the XS method is simpler than MB but the exact knowledge of the field free point (FFP) motion is essential for its implementation. Our simulation work focuses on the implementation of a new approach for MPI reconstruction: it is called hybrid x-space (HXS), representing a combination of the previous methods. Specifically, our approach is based on XS reconstruction because it requires the knowledge of the FFP position and velocity at each time instant. The difference with respect to the original XS formulation is how the FFP velocity is computed: we estimate it from the experimental measurements of the calibration scans, typical of the MB approach. Moreover, a compressive sensing technique is applied in order to reduce the calibration time, setting a fewer number of sampling positions. Simulations highlight that HXS and XS methods give similar results. Furthermore, an appropriate use of compressive sensing is crucial for obtaining a good balance between time reduction and reconstructed image quality. Our proposal is suitable for open

Introduction
Magnetic particle imaging (MPI) is a recent medical imaging technique dating back to 2001 and published four years later in Gleich and Weizenecker (2005). This technique is capable of acquiring medical information recovering the distribution of superparamagnetic nanoparticles, used as tracers inside the body, by means of the measured induced signal. For this reason MPI shares similarities to nuclear medicine methods such as PET or SPECT, where an external contrast agent is needed, but MPI is ionizing radiation free.
The key element for the MPI acquisition technique is given by the peculiar non-linear behavior of superparamagnetic nanoparticles under the action of magnetic fields (Néel 1949). The basic idea of MPI consists in applying a suitable combination of a static magnetic field with a time varying one to the superparamagnetic material. The first field, usually denoted selection field, is characterized by a high spatial gradient, while the second, called drive field, is a uniform field: the effect of their superposition is a varying field having a null moving point denoted as field free point (hereafter FFP).
An additional time dependent field, with a lower frequency, called focus field, may be superimposed enlarging the region to be scanned. With these fields, the signal induced by the superparamagnetic substance will be mainly generated by the nanoparticles located in a neighbourhood of the instantaneous FFP position, being the other particles saturated. The neighbourhood size depends on the gradient strength: a narrower active region around the FFP results from larger values of the gradient intensity.
The induced signal is then always related to the instantaneous position occupied by the FFP. If the FFP trajectory is well designed, the entire region of the superparamagnetic material is scanned, and a complete MPI image of concentration values can be obtained. Actually, the problem of deriving the concentration of nanoparticles by their induced signals, usually referred as reconstruction problem, is crucial in MPI signal analysis. In literature there are different techniques that accomplish this. Some methods make use of the superposition principle, and state that the time Fourier transform û of the induced signal is expressed by a linear relation in terms of the concentration function c u Sc,= where Ŝ is a linear operator or matrix usually called system function . The matrix Ŝ depends on the instrumental apparatus and on the superparamagnetic substance. By measuring û and solving (1) the concentration c can be recovered. Usually problem (1) is overdetermined and ill-conditioned, so the solution must be computed in the least squares sense, and some regularization techniques based on the singular value decomposition (SVD) or the Tikhonov's method have to be applied .
The matrix Ŝ can be obtained experimentally, by measuring the induced signals when a known amount of superparamagnetic material is concentrated in a single pixel or voxel, for every pixel or voxel, or theoretically, by means of a formula proposed by , which models the nanoparticles non-linear behaviour with the Langevin function (Knopp and Buzug 2012). The first method, usually known as measurement-based (MB) (Graeser et al 2012, Grüttner et al 2013, requires long calibration procedures for a full acquisition of the system function. For the model-based method a shorter calibration, consisting of receive chain characteristics evaluations, is required. Nevertheless, the results in this case are less accurate than those obtained by the measurement-based method, since superparamagnetic nanoparticles may consistently deviate from the ideal behavior of the Langevin function. In addition, the various sources of noise affecting the MPI signal are difficult to be modeled and inserted in the theoretical formulation of the system function. A different approach for image reconstruction was followed by Goodwill and Conolly (2011). If the selection field has a constant spatial gradient, and the drive-focus fields are uniform in space, the nanoparticles concentration inside the FFP is proportional to the induced signal over the FFP velocity. In the one dimensional case this can be expressed by the following simplified relation: This reconstruction technique, called x-space (XS), is simpler than those methods based on the system function, but it requires the knowledge of the FFP exact position and velocity at each time step of the scanning process. Our present work treats image reconstruction techniques in the context of an open geometry human size apparatus. The realization of such devices requires specific calibration procedures related to the methods used in the reconstruction phase. One of the advantages of open bore scanners would be an easier access to the patient, especially in interventional scenarios with simultaneous and real-time scanning processes. In this case of geometry configurations with a larger field of view (FOV), the exact velocity gridding for x-space MPI could be difficult to achieve during the whole scanning process. Hence, our proposal is to define an innovative technique, which we call hybrid x-space (HXS) (Tateo et al 2015), resulting from the combination of the measurement-based and the classical x-space approach, sharing typical features common to both methods. HXS is based on the x-space (XS) reconstruction because it needs the knowledge of the FFP position and velocity for each time instant. The difference with respect to the original x-space formulation is that we directly estimate the latter from the experimental measurements done during the calibration scans, typical of MB approach. However, our method should not be considered as an extension to the MB method because it does not consider the reconstruction inverse problem as in that setting. It is worth noting that the calibration of the apparatus is time consuming. Even if it is only periodically needed, in case of 3D FOVs this procedure may lead to an acquisition time of several hours, or for extreme cases up to one or two days, depending on the number of grid positions (Halkola et al 2012). In order to reduce the number of calibration scans, in Knopp and Weber (2013b) the authors suggest the use of compressive sensing. Following this idea, we apply to HXS a compressed sensing based on the use of circulant matrices, since in Liang et al (2009) and Yin et al (2010) it has been asserted that they allow significantly faster compressive sensing reconstruction compared to using independent and identically distributed (i.i.d.) random matrices. It is worth pointing out that our aim is not to provide an alternative method to MB, but instead to present a new approach similar for simplicity to XS. While the XS requires the complete time-trace processing, in our case the processing time depends only on the number of the initial sampling positions. For this reason we will show that the HXS performance is comparable to the one from the classical XS method. In addition, all the pre-processing and post-processing techniques applicable to the XS method can be adopted as well for the HXS setting.
The article is organized as follows. In section 2 we explain the simulation scheme and the parameters used to generate the MP image. In section 3, after a more detailed description of the MB and XS reconstruction techniques, and an introduction to compressive sensing and circulant matrices, we expose our new HXS reconstruction method. In section 4 we test the performance of this method in a series of simulation studies for given phantom images comparing the results between HXS and XS. For completeness, we report the results also for the MB approach. Then, we present our conclusions in section 5.

Simulation of magnetic particle imaging system
Our analysis was carried out by simulating all the steps necessary for the MPI process. Indeed, following the scheme in Weizenecker et al (2007) we generated the signal to be converted into the final MP image. The same simulated setting was used for the implementation of the different MPI reconstruction methods. In this section we provide the parameters used for the entire simulation process.
In our algorithm we have simulated the following components of an MPI apparatus: • selection field; • excitation fields; • FOV and the nanoparticles used as tracer; • receive coils.
Firstly, we have assumed the gradient G of the selection field H S to be uniform in the FOV. In particular, the gradient values were set to 5 Tm 1 0 1 µ − − and 2.5 Tm 1 0 1 µ − − along the x and y direction respectively. In this configuration, as no additional fields are applied, the FFP is in the center of the FOV.
The motion of the FFP in the FOV was simulated superimposing two additional varying fields, the focus/drive fields, to the selection one. The coupling drive-focus fields is used for extending the FOV to a larger area. The drive and focus fields are characterized by different excitation frequencies, indeed the drive field one is higher. Different types of FFP trajectories can be used for MPI, as described in Knopp et al (2009). In our case we have chosen cartesian and Lissajous trajectories, depending on the frequency ratios given by for the drive field, and to N 100 D = for the focus field respectively. We have simulated the MB method using continuous cartesian trajectories for the focus field, and Lissajous trajectories for the drive field as in Knopp et al (2009). Instead, for the XS and HXS methods continuous cartesian trajectories have been used for both fields in a joint step. This choice permits to achieve the best performance for each method, as known from the literature. The amplitude of the drive field was chosen so that the FFP moves in a region of size 1 1 × cm 2 . The amplitude of the focus field was fixed to translate that region to the entire FOV of size 10 10 × cm 2 . The nanoparticles behavior was modeled according to the Langevin theory of paramagnetism as mentioned by Caizer (2003). Denoting by c x ( ) the concentration of nanoparticles in the position x, the total magnetization M is given by where L is the Langevin function defined as and e H is the direction of the magnetic field strength. The equation (3) includes a parameter β that is specifically related to the type of nanoparticle used as tracers. This value is is the magnetic moment of a single nanoparticle, and k B is the Boltzmann constant. In our simulation we have chosen M 477 S = kA m −1 for the saturation magnetization of nanoparticles, D = 25 nm for their diameter, and T = 310 K for the temperature.
To estimate the signal of nanoparticles we have considered two recording coils: the first one oriented in the x-direction and the second in the y-direction. The coils are circular with radius R = 20 cm and they are placed 15 cm away of the center of the FOV. During a single scan we have estimated two signals, one for each coil, corresponding to each axis direction.
The goodness of MPI reconstruction was estimated using the Pearson correlation coefficient between the reconstructed image and the original nanoparticle concentration in the FOV, and also by visual inspection.

Methods
In the following sections we provide a brief description of the two MPI reconstruction techniques that we combined to develop the new method described in section 3.3.

Measurement-based method
The theoretical model describing the magnetic nanoparticles behavior is based on the assumption that their magnetization follows the Langevin function, and that the effects due to the relaxation times of the particles are negligible. This is only a good approximation of the actual behavior. For this reason, Gleich and Weizenecker (2005) proposed a method based on the discretization of the system function through a calibration procedure that has the advantage of not requiring any hypothesis. The MPI equation can be written as a linear integral equation with a kernel given by the system function; it describes the signal induced by the particles in a receiving coil as a function of space. The signal is usually considered in the frequency domain in order to remove the fundamental frequency with a high pass filter. The relationship between the signal spectrum and the particles distribution can be written as ( ) denotes the mean magnetic moment, c x ( ) is the concentration and T R is the repetition time. The signal is discretized by sampling the FOV in N discrete posi- In the matrix-vector notation: where M depends on the total number of considered frequencies. Operationally, the MPI system matrix is obtained by acquiring the system response to a nanoparticle sample of known concentration c 0 and volume V S , placed in each x i position of the FOV by a robot. The ith element of the system matrix can be obtained as: It is worth noting that the system function construction is deeply linked to the nature of the nanoparticles and MPI scanner setup: hence, after any change, a new calibration process is needed (Grüttner et al 2013). Given that (8) is an ill-posed problem, a classical way to solve it is the least squares approach: or the weighted least-squares problem: Each weight w j in the matrix W is the Euclidean norm of the jth row of Ŝ, while ũ is the measurement vector û distorted by noise. The solution of the least squares problem can be calculated by using the singular value decomposition (SVD) of the matrix S W S w 1 2= : with u W u w 1 2= . We use a truncated singular value decomposition (TSVD) based on visual inspection for the selection of the truncation index. Finally, the reconstructed image is given by x c IMG( ) = .
3.1.1. Compressive sensing and circulant matrices. For the MB method the calibration procedure is essential but also expensive, since it may require several hours or in extreme cases one or two days. As previously pointed out, an MPI device could use different types of nanoparticles and/or change some parameters related to the drive and selection field. As a consequence, the system function must be recalculated requiring a new calibration procedure: this is certainly a drawback of the MB approach. Following the idea introduced in Knopp and Weber (2013a, 2013b) we consider the compressive sensing technique for the data acquisition, as it allows to recover information from highly reduced measurements (Donoho 2006) and simultaneously to decrease impressively the calibration time. The compressive sensing is based on two important assumptions: the signal is sparse or compressible, that means only few components are greater than a threshold value, and the system matrix is incoherent. In this case the reconstruction by means of undersamples is feasible overcoming also the Nyquist-Shannon sampling theorem. Therefore, the first condition is obtained with the sparsification of the MPI system matrix using a basis transformation, as the discrete Fourier transform (DFT) or the discrete cosine transform (DCT), while the second one is assured by an appropriate distribution of the sampling positions (Lampe et al 2012). Practically, there exist different types of random sensing matrices such as Bernoulli, Gaussian or Subgaussian: they all guarantee the low coherence of the system function matrix. However, in Yin et al (2010) the author asserted that using Toeplitz and circulant sensing matrices allows significantly faster compressive sensing reconstruction compared to using i.i.d. random matrices. The use of Toeplitz and circulant matrices was firstly proposed for compressive MR imaging in Liang et al (2009), since it showed better reconstruction results than those obtained using the Fourier schemes. It is worth pointing out that the advantage in the use of circulant matrices is in carrying out only one sampling for the data captured by the device, reducing then the number of samplings for the system matrix. Thus, we consider an m N × partial circulant sensing matrix is a selection matrix, and C is a block-circulant matrix with each block being a circulant matrix (see Yin et al (2010) and Mazzia and Sorice (2015)). Hence, let be s S x and M smaller than N, in order to reduce the scan acquisition, then the linear system to solve is where B c is the reduced basis transformation matrix and v k is the solution to compute. Since the problem (15) is a strongly under-determined linear system, then one should try to solve the following optimization problem . In literature there are several algorithms solving optimization problems, as the fast iterative shrinkage-thresholding algorithm (FISTA) (Beck and Teboulle 2009), applied also to MPI reconstruction in Lampe and Voss (2008). However, for the particular structure of the operators and block circulant matrices, we applied the MATLAB toolbox RecPC, developed by Yin et al (2010), considered to be a better choice for circulant matrices.

X-space method
The main advantage of the x-space is to be a real-time reconstruction method. The nanoparticle concentration in the FOV is estimated by a simple division between the signal acquired and a scalar. A detailed description of this theory is given by Goodwill and Conolly (2010) for the one dimensional case. Here, the nanoparticles magnetization M, induced by the applied magnetic field H, can be modeled as in (3) by the Langevin function L. If the field H is linearly dependent on the space, with a linear gradient G, the signal expression becomes where x ṫ s ( ) is the FFP velocity, ρ is the nanoparticles density and B 1 is the inductive detector sensitivity. Finally, by simply dividing out the extra terms and assuming the derivative of the Langevin curve to be the system point spread function (PSF), the MPI native image is given by The 1D MPI image equation (18) can be generalized to the multidimensional setting as in Goodwill and Conolly (2011) In the last equation the final image is computed by gridding only the component u t ( ) ∥ of the signal, namely the one recorded by the receive channel whose axis is collinear with the FFP trajectory. In order to improve the reconstruction result, in our simulation we computed two scans using two different cartesian trajectories, along both the orthogonal directions x and y. Thus, we generated two images, one for each recording channel respectively, and the final reconstructed image is obtained as a geometric averaging of the previous generated images. Note that the application of this method is based on the exact knowledge of the FFP motion (position, velocity) during all the scanning process.

Hybrid x-space
As mentioned in the previous section, to use the x-space method the knowledge of the FFP motion is needed. Here we propose an alternative approach to experimentally estimate the FFP position and velocity. Our method, called hybrid x-space (HXS), differs mainly from the XS method, since the complete time-trace no needs to be processed for the MPI reconstruction. Furthermore the application of the HXS method would be useful in solving related hysteresis problems, or in cases of apparatus with open geometries or large dimensions when incidental factors could make the currents, the fields and so the trajectory irregular. The method is defined hybrid since, like in the XS, the knowledge of the FFP motion is needed for the reconstruction. Moreover, as in the MB, calibration scans for the experimental measurements of FFP position and velocity are used. In HXS we map the values needed for the MPI x-space reconstruction only to each calibration position inside the FOV. We are able to compute the concentration values on a fewer number of points, since only the gridding of the signal in the sampling positions is needed, whereas the x-space requires that the complete time-trace has to be processed.
More specifically, for each calibration point we need to estimate three parameters: the instant of time t * , the corresponding estimate of the FFP position t x s ( ) * , and the estimate of the denominator t ẋ s ( ) * in the equation (19), which involves the sensitivity of the receive channel coils, the nanoparticle characteristics (magnetic moment and temperature), and the FFP velocity.
Firstly, we sample the MPI signal for building the system matrix S as required by MB. In particular, the ith column of the matrix i S :, ( ) is the signal response to a known delta sample of particles located in the x i position of the FOV. Our idea is to consider the maximum along this column to estimate the time t i * when the FFP is closer to x i during the scanning. The denominator in the equation (19) at time t i * is computed by dividing the signal peak by the known concentration of the sample.
The result of this phase, which is done only one time, is the set of estimates  We point out that, as in the XS, the reconstruction is done using only the signal induced in the coil collinear with the trajectory. As in the previous section, the final image is the geometric average of the images computed for each continuous cartesian trajectory along the x and y axes.
The drive-focus field superposition may cause more peaks appearing in the signal response associated with each sampling position: this can be explained by the FFP passing in its proximity at different times along the trajectory. In this case, column-wise, a triplet for each peak is considered, and the MPI image value is estimated by averaging among them. Thus the equation (20) is replaced by where M i is the number of peaks corresponding to the ith calibration points, and t i j , ( ) * is the instant of the jth peak related to the x i position.
It is worth noting that the hysteresis problem may be solved in this way, because the time instant t i * is related to the instantaneous passage of the FFP through the x i position, as well as to the relaxation and induction times. In summary, figure 1 shows the flowchart of the hybrid x-space algorithm.

Results and discussion
Following the scheme described in section 2, we simulated the entire process for generating the MPI signal, and reconstructing some MPI images. All the parameters needed for the simulation were fixed to the values presented in the same section. The calibration procedure was carried out choosing a uniform grid of 40 40 × positions covering the whole FOV. In this situation 1600 calibration scans are required to generate the MPI system matrix. During the simulated acquisition scans, a sample of known concentration c 0 was placed consecutively in each one of the 1600 FOV positions. This operation is time consuming in real cases as well as in simulated ones. For this reason, a local computer farm Bc 2 S (www.recas-pon.ba.infn.it) was used to compute the system matrix. The Bc 2 S is a distributed computing infrastructure that allows a storage limit of 1.8 PB, and comprises about 5000 CPU. A single simulation job outputs two inducted signals, one for each axis direction respectively. In this way, the calibration phase produces as results the two system function matrices SF x and SF y . The two system functions were used for our MB and HXS implementation. In our analysis we compared the two MPI reconstruction methods: XS, and HXS. For the sake of completeness we also report the results obtained using the MB approach. To assess the performance of each method, we used different particle distributions in the FOV. For each one of them, and for each method, a reconstructed image was computed, and the correlation coefficient estimated. The consistency for each method is expressed by the sample mean and standard deviation computed on all the correlation coefficients. Some 'phantom' examples of the nanoparticle distributions used for the simulations are shown in figure 2, together with the axes orientations considered in our analysis.

Performances
The MB method estimates the 1600 concentration values on the FOV grid. Our images were composed by 10 6 pixels, so it was necessary to assign a concentration value to all the remaining ones. That was done applying a two-stage approximation method based on bivariate B-splines, implemented in a C library called TSFIT (www.tsfit.de) Zeilfelder 2004, 2005). Specifically, a MATLAB mex interface to this library has been used as in Iurino (2014) to compute the final image. The reconstructed images in figure 3 were obtained using the procedures described in section 3.1. The result of this method is a reconstruction of the 'phantom' objects with a mean and standard deviation correlation equal to 0.96 and 0.01 respectively. For the implementation of the x-space method the knowledge of the FFP motion is needed. In our simulation study that was obtained by computing, at each instant of time, the FFP position in the FOV during the scanning process. This procedure is necessary to emulate the FFP temporal mapping, which is instead needed in real cases. Two directional images were reconstructed, one for each receive coil respectively, as described in section 3.2. The next step of the x-space method involves the extension of the concentration values to all image pixels, so we have used the TSFIT library as described above. To improve the results, each couple of directional images was combined by geometric average, and then post-processed. The XS image reconstructions are shown in figure 4. Also in this case the reconstruction performances were evaluated in terms of the mean correlation coefficient, which is now 0.74 with a standard deviation of 0.03.
To apply the HXS method described in section 3.3, we estimated for each direction the tri- . Then, for each particle distribution, the two induced signals given by equation (21), collinear to the cartesian trajectory along x and y-axis respectively, were used to compute two direcctional images. The final result is given by their geometric average. In this case the resulting mean correlation coefficient is 0.76, with a standard deviation of 0.03. Figure 5 displays some of the reconstructed 'phantom' images for visual inspection and comparison.

Compressive sensing
Since our new approach is based on a time consuming calibration procedure, we applied the compressive sensing technique to reduce the number of calibration points, and to shorten the calibration process. It is worth noting that, removing too many grid points causes significant loss in the reconstructed image quality, thus a good compromise between maximum reduction and correlation index must be accomplished. For this goal several tests reducing the number of grid points used for the apparatus calibration were carried out. In particular, rates of 50%, 40%, 30%, 20%, 15%, 10%, and 5% of the 1600 original sampling positions were considered in our tests. Missing data reconstruction was computed by compressive sensing, while the choice of the sampling positions was achieved by means of the circulant matrix method as described in section 3.1.1. Figure 6 displays the reconstructed images of the 'phantom' bar object with respect to the different sampling ratios. Finally, the boxplots shown in figure 7 suggest that, using the compressive sensing, a 85% reduction of sampling positions does not generate a decrease in the reconstruction quality. All reduction rates less than 85% result in a mean and standard deviation correlation coefficient of about 0.76 and 0.03, respectively. Indeed, for a larger percentage reduction, the mean correlation index decreases and the standard deviation grows up. Hence we may infer that reduction rates above 85% make the reconstruction process unstable.

Resolution
The spatial resolution is a crucial parameter for evaluating the performance of the reconstruction method. It refers to the minimum distance between two points that allows to consider them as distinct objects. It is well known that, for a fixed reconstruction method, the greater is the gradient the better is the resolution. Since in this simulation study we considered two different gradients for the x and y directions respectively, we expected to have different resolution values for each direction. A test to estimate the resolution of the HXS method here proposed was carried out. We used couples of vertical and horizontal concentration bars, placed at different distances. For the x-axis direction different distance values of 6 mm, 4 mm, 3 mm, and 2 mm were used, while for the y-axis direction distance values were set to 8 mm, 7 mm, 6 mm, and 5 mm. Figures 8 and 9 show the plots of the grey levels corresponding to the central line of the FOV along the axis direction under investigation. These resolution tests highlight that in our simulated setting condition, the measured resolution of the HXS approach is approximately 4 mm along the x-axis and 6 mm along the y-axis.

Conclusion
In this work we present a new approach for MPI reconstruction, called hybrid x-space. HXS is an alternative approach to experimentally estimate the FFP position and velocity. The name hybrid reflects the fact that the knowledge of the FFP motion, needed for the reconstruction as in XS, is estimated by a calibration scan as in MB. Our method differs mainly from the XS method, since the complete time-trace no needs to be processed for the MPI reconstruction. Furthermore, the application of the HXS method would be useful in solving related hysteresis problems, or in cases of apparatus with open geometries or large dimensions when incidental factors could make the currents, the fields and the FFP trajectory irregular. The problem of reducing the long calibration times, faced in the MB by using compressive sensing with random matrices, has been solved in the present work by compressive sensing techniques based on circulant matrices. A comparison of XS and MB methods with the HXS approach in the MPI framework was made by numerical tests, and the accuracy was assessed by computing the sample mean and standard deviation overall the correlation coefficients. From the results of our analysis, the MB approach seems to perform better than the proposed and the XS approach. We believe that this could be related to the parameters used for simulating the entire MPI process and to the XS approach implemented. For example to scan a wider FOV region, we have used the focus-drive field superposition, while another possible solution is the partial-FOV approach. As it has been shown in Goodwill and Conolly (2011), this can improve the performance of the XS reconstruction, and, as a consequence, of the HXS reconstruction. A specific analysis will be conducted in a future work to improve the reconstruction performances considering XS different approaches known in literature. We highlight that the goal of our work is to compare HXS to the classical XS method: the results of our simulations showed a comparable performance between them. Hence, any improvement obtained for the XS implementation will be also valid for HXS. Moreover, the compressive sensing technique gives a 85% reduction of the calibration time without compromising the performance itself.  However, we believe that this work is a first attempt to apply the proposed approach for MPI reconstruction, and that there is still much to improve.