PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

 
Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
 
IEEE Trans Med Imaging. Author manuscript; available in PMC 2014 April 17.
Published in final edited form as:
PMCID: PMC3990467
NIHMSID: NIHMS377771

Multidimensional X-Space Magnetic Particle Imaging

Abstract

Magnetic particle imaging (MPI) is a promising new medical imaging tracer modality with potential applications in human angiography, cancer imaging, in vivo cell tracking, and inflammation imaging. Here we demonstrate both theoretically and experimentally that multidimensional MPI is a linear shift-in-variant imaging system with an analytic point spread function. We also introduce a fast image reconstruction method that obtains the intrinsic MPI image with high signal-to-noise ratio via a simple gridding operation in x-space. We also demonstrate a method to reconstruct large field-of-view (FOV) images using partial FOV scanning, despite the loss of first harmonic image information due to direct feedthrough contamination. We conclude with the first experimental test of multidimensional x-space MPI.

Index Terms: Angiography, biomedical imaging, magnetic particle imaging, magnetic particles

I. Introduction

Magnetic particle imaging (MPI) is a new medical imaging tracer modality that holds significant promise for high-contrast applications in human or small animal angiography, cancer imaging, in vivo cell tracking, and inflammation imaging [1]. The technique exploits the nonlinear magnetic characteristics of iron oxide nanoparticles to generate an image whose resolution is defined by the magnetic properties of the nanoparticle and the magnitude of the localizing magnetic field gradient. Like magnetic resonance imaging (MRI), the spatial resolution of MPI is much finer than the wavelength of the electromagnetic fields used to interrogate the magnetic nanoparticles. MPI uses no ionizing radiation and has ideal tracer imaging contrast since there is no background signal from tissue and because tissue is transparent to low frequency magnetic fields.

Recent MPI developments include real-time MPI in a mouse [2]–[5], methods for single-sided MPI using strongly nonlinear magnetic field gradients [6]–[9], and studies on optimal MPI particle size [10].

A. MPI Signal Theory

Linearity and shift invariance (LSI) are important characteristics of most imaging systems [11]. There has been significant work to understand MPI as an LSI system. Non-LSI systems can be difficult to analyze using standard signal processing techniques such as convolution. Prior papers in this area have relied on the temporal harmonic domain [1]–[4], [6]–[9], [12], [13]. These papers analyze the harmonics of the received MPI signal. Nonlinear magnetic nanoparticles respond to a sinusoidal magnetic waveform with harmonic signals at multiples of the excitation frequency. These harmonics are suppressed sufficiently far from the center, or field-free-point (FFP) of a field gradient, since the gradient field leaves the particles in saturation despite the RF excitation. Hence, the gradient field provides a method to localize harmonic response in 3-D space. The most complete theoretical treatment of this approach is seen in Rahmer et al. [12], where it is shown that the 1-D frequency-space signal can be described using Chebychev polynomials of the second kind convolved with the magnetization density. The Chebychev polynomial model is exact in one dimension, but extension to two and three dimensions is an approximation.

A fundamental assumption of all harmonic methods is that each pixel is interrogated over several cycles of the RF excitation. This is not accurate for faster scanning methods, where a single pixel is scanned instantaneously only once. We recently developed a signal model that does not necessarily require a repeating excitation, describing the 1-D MPI imaging process as an instantaneous scan through x-space rather than a sinusoidal steady-state harmonic decomposition [14]. The key goal of this paper is to extend our 1-D x-space formalism to 2-D and 3-D.

B. Reconstruction in MPI

Current reconstruction techniques in MPI [1], [3], [4], [8], [9], [12], [13] require a precharacterization of the magnetic nanoparticles whose signal response is formulated into a system matrix. The system matrix is comprised of Fourier components of the temporal signal for every possible location of a point source. For example, for an image with Nx × Ny × Nz possible points, the total number of elements in the system matrix will be N = NxNyNzNcNf where Nc is the number of receive coils and Nf is the number of Fourier components desired for the reconstruction. The system matrix can be measured physically using a nanoparticle sample [2], or estimated using a model [9]. However, this means that the system matrix is specific to the nanoparticle sample, and reconstruction will be less accurate if the nanoparticle behaves differently in tissue, if the system drifts, or if the model is inaccurate. Reconstruction is achieved through regularization and matrix inversion techniques such as singular value decomposition or algebraic reconstruction. Care must be taken when regularizing the solution to achieve high resolution while not amplifying noise when inverting the system matrix [4]. It is important that MPI be subjected to well-conditioned image reconstruction [15] to avoid any loss of signal-to-noise ratio (SNR).

To date, there have been two approaches other than x-space to understanding MPI without a system matrix. In Goodwill et al. [16], we described and built a narrowband MPI system that imaged multiple MPI harmonic mixing products and placed them on a grid in real-space. However, we were having difficulty reconstructing the multiple images into a composite image so that the reconstructed image is linear with the nanoparticle density irrespective of the shape of the phantom. A second approach to understanding MPI without a system matrix is seen in Schomberg [17], where the author also approaches the MPI process using an adiabatic assumption. The author’s theoretical approach is general and finds that the MPI signal is closely related to a convolution operator and has parallel goals to the approach presented here. In this paper, we present a simple theoretical formalism for MPI and validate it with simulation and experiment. We also present a fast reconstruction algorithm that computes the MPI image without matrix inversion and without a model based image reconstruction, extending the 1-D image reconstruction method published in [14] into 2-D and 3-D.

II. Hypotheses for Multidimensional X-Space Magnetic Particle Imaging

In our recent paper [14] we employed reciprocity and LSI imaging systems theory [18], [19] to analyze the 1-D MPI signal imaging process. We assumed that the nanoparticle magnetization instantaneously aligns with the applied local magnetic field. We showed that the MPI signal in one dimension is linear and space invariant and can therefore be described as a convolution. We found analytically and experimentally that the point spread function (PSF) is the derivative of the magnetic nanoparticle’s Langevin function. This analysis provided estimates for bandwidth requirements, which approach a megahertz for typical imaging parameters. We also analyzed MPI’s SNR, and explored the limits as we scale MPI to human sizes. We found that the limit to SNR will be patient heating, and the limit to the (partial) field-of-view (FOV) will be magnetostimulation.

In this paper, we extend the 1-D x-space MPI theory into two and three dimensions. We again hypothesize that the magnetic nanoparticles align instantaneously with the local magnetic field and that the loss of first harmonic information due to direct feedthrough contamination is recoverable. We need one additional hypothesis for multidimensional x-space, namely that the linear 3-D gradient field can be written as Gx where G is an invertible matrix so that the gradient field uniquely identifies the location x in three-space. After proving that any real-world gradient field is invertible, we then prove that 3-D MPI is a linear and shift invariant imaging process. We derive the analytic 3-D point spread function of MPI. We introduce a fast image reconstruction algorithm that requires no calibration measurements or matrix inversion, so it is both computationally efficient and robust to noise. To apply the x-space formulation to real MPI systems, we discuss how the loss of the fundamental frequency breaks the strict LSI properties. We hypothesize and provide experimental evidence that the lost first harmonic information is fully recoverable using robust and noise-free image processing methods. Finally, we conclude with experimental 2-D demonstration of x-space MPI. These three hypotheses are justified in detail below.

A. Hypothesis 1: The Instantaneous FFP Is Uniquely Defined in Space

MPI relies on a 3-D linear gradient in the form [20]

equation M1

where the vector equation M2 denotes position in real space, and the parameter α [set membership] (0, 1). Note that the trace(G) = 0, which is consistent with Maxwell’s equations in a source-free space ([nabla]2B = 0) [21], [22]. For the very common case of a cylindrically symmetric Maxwell z-gradient, used for all experiments herein (see Fig. 1), α = 1/2 [20] and the off-diagonal elements are zero. So the gradient matrix G is diagonal

equation M3
(1)
Fig. 1
Two opposing ring magnets with radial symmetry about the z axis produce a 3-D gradient field with a FFP at the isometric center. This gradient can be remarkably strong. Our current imager produces a gradient of 6 T/m in the z axis, and 3 T/m in the x ...

Typical small animal NdFeB Maxwell pair gradient strengths are μ0Gzz ~ 2500 [2] to 6000 mT/m [16]. MPI spatial resolution is anisotropic and twice as fine in z as in the x and y direction due to this fundamentally anisotropic magnetic field gradient.

In addition to the 3-D field gradient, we can add an orthogonal set of homogeneous magnets that produce static and time-varying fields to shift the FFP

equation M4

These homogeneous magnets could be built using Golay coils or Helmholtz coils and these are considered to be accurately modeled as homogeneous over the linear region of the gradient field. Giving the gradient a convenient negative sign, the total field can be described as

equation M5

We can solve for the instantaneous location of the FFP, xs(t), such that H(t, x) = 0, provided that the G matrix is nonsingular

equation M6

In this paper we hypothesize that the G matrix is nonsingular. Here we prove that this is true for virtually all realistic gradient coils. First, consider the Maxwell gradient pair, the design of choice for virtually all experimental work in 3-D MPI. It is clear that the diagonal G matrix is always invertible; indeed the equation M7. For the more general case, assuming cross terms Gxy = Gxz = Gyz = 0, the equation M8 is guaranteed to be nonzero, except for trivial cases. Finally, if we assume only that Gxy = 0, then one can show that equation M9. Since the coefficients α and 1 − α are guaranteed to be positive, the det(G) is guaranteed positive. Hence, the G matrix is not singular, and there always exists a unique FFP.

Of course, the uniqueness of the FFP could be lost if the region of interest includes regions outside the linear region of the gradient coil. Artifacts could occur in such a case, akin to “fold over” artifacts in MRI, where the body extends beyond the monotonic region of the gradient field. Hence, to avoid this challenge we will assume that the entire FOV of the MPI scan is within the linear region of the gradient field and also within the homogeneous region of the shifting magnets.

Then, the magnetic field at an arbitrary point x is related to the instantaneous position of the FFP

equation M10
(2)

B. Hypothesis 2: Adiabatic Langevin Model

The MPI signal is due to the nonlinear response of a superparamagnetic iron oxide (SPIO) nanoparticle to a changing magnetic field. At dc field strengths used in MPI of Bmax < 1 T, tissue is largely unaffected by the magnetic field, but a SPIO particle undergoes a nonlinear change in magnetization described by the Langevin theory of paramagnetism [23], [24]. For a density of SPIO nanoparticles, ρ [particles/m3], the Langevin equation gives a magnetization density

equation M11
(3)

where m[Am2] is the magnetic moment of a single magnetic nanoparticle and L is the Langevin function. The field required for saturation equation M12 [A/M] is given by the magnetic moment, Boltzmann constant kB and temperature T. For a spherical particle, the magnetic moment can be computed as m = Msat(π/6)d3, where μ0Msat ~ 0.6 T for magnetite and d[m] is the particle diameter [23].

The vector direction assumes that the magnetic particles align adiabatically and instantaneously with the applied magnetic field, which is strictly true only if the time varying magnetic field is much slower than the particle relaxation time. Neel and Brownian relaxation of the particles will reduce the magnetization and change the phase between the applied magnetic field vector H and the measured magnetic moment M(H) [25]. We note that typical Brownian time constants of most magnetic nanoparticles used in MPI are about 1–30 µs [25] whereas typical MPI scanning frequencies are below 25 kHz, so this appears physically realistic. The adiabatic hypothesis is required for strict shift invariance; artifacts from relaxation appear to be mild. This model is discussed in the well-developed field of ferrohydrodynamics, which aims to predict the behaviour of ferrofluids to time-varying magnetic fields [23], [24]. This hypothesis is used extensively in prior approaches to MPI theory [4], [12]. Here, we assume throughout that the data acquisition dwell time per pixel exceeds the relaxation time constant. As shown in our experimental results in Section IV-B, the x-space approach generate excellent images consistent with the adiabatic hypothesis. We are continuing to study the limitations of this model.

C. Hypothesis 3: Loss of Low-Frequency Information Is Recoverable

In MPI, the RF transmit occurs during signal reception. This means that the received signal is contaminated by direct feedthrough from the source RF coil coupling into the detector, despite meticulous efforts at electronic and geometric decoupling. Hence, all current MPI imaging methods must be reconstructed only from (uncontaminated) high frequency information. Here we analyze the MPI signal equation as if the complete receive bandwidth were available; later below we demonstrate that this lost low-frequency information represents the low spatial frequencies (e.g., dc or baseline component) of the image. Fortunately, with just a small amount of overlap in partial FOV scans, it is simple to reconstruct a smooth and contiguous version of partial-FOV scans over the entire FOV using robust image processing methods. This allows us to recover the lost baseline information without adding a significant amount of noise. It is important to note that all current MPI scanning methods suffer from this lost dc or baseline information so this challenge is in no way unique to x-space MPI scanning.

III. Multidimensional Theory of MPI

Now consider a continuous distribution of magnetic nanoparticles with density ρ(x) [particles/m3]. From (2) and (3), we note that we can write the magnetization density of ρ(x) nanoparticles located at position x when the FFP is at xs(t)

equation M13
(4)

and thus the total dipole moment is obtained by integrating the magnetization across the imaging volume

equation M14

This total dipole moment can be written as a spatial convolution interrogated at the instantaneous FFP location

equation M15

For an imaging system using an inductive detector, we can use reciprocity to calculate the received signal [26]. For simplicity, let us assume orthogonal receive coils aligned with the x, y, and z axes of the instrument. Then, the sensitivity of the receive coils, −B1(x) [T/A], would be a matrix of sensitivities. For the case for receive coils in each of the x, y, and z axes, respectively, the sensitivity matrix would be equation M16. From reciprocity, the received signal vector is

equation M17
(5)

To evaluate this derivative, let us begin by building a set of tools to evaluate the MPI signal equation. We begin by defining

equation M18
(6)

We can decompose r into a tangential component, r, and a normal component, r[perpendicular] = (rr). Rewriting r yields

equation M19

The derivative of the quasi-static Langevin function with vector-valued, time varying operand r

equation M20

can be rewritten as a function of r and r[perpendicular]

equation M21
(7)

So, we see that the derivative of the Langevin curve has two components, each proportional to the tangential component or normal component of the FFP velocity vector.

We now use these tools to calculate the derivative of the signal equation. From (4) and (5) and definitions (6), we can rewrite the MPI signal as

equation M22

and evaluate the derivative using (7)

equation M23

Substituting for r gives us a familiar convolution integral

equation M24

which yields

equation M25

This form is not yet a simple PSF, and so we factor out the FFP velocity magnitude ‖xs‖ and the FFP velocity unit vector equation M26 using an outer product vector identity (a · b)b = bbTa. This yields the

equation image

with PSF

equation M27

The signal equation indicates the inductively received signal is the matrix multiplication of a matrix PSF with the FFP velocity vector, equation M28. The PSF does not change as a function of the FFP velocity magnitude, ‖xs‖. If we only consider the radially-symmetric scalar components of the PSF, we find the PSF envelopes

equation M29
(9)

equation M30
(10)

These important envelopes seen in Fig. 2 give us the maximum attainable resolution in MPI. The higher resolution Tangential Envelope, ENVT, is the derivative of the Langevin equation and is described in detail in Goodwill et al. [14]. ENVT defines the intrinsic resolution and bandwidth requirements for MPI. The lower resolution envelope, ENVN, is unique to generalized MPI and has a full-width at half-maximum (FWHM) that is wider. We can solve for the FWHM of both envelopes analytically as a function of Hsat or, alternatively, in terms of the particle diameter d

equation M31
(11)
Fig. 2
The tangential and normal point spread function envelopes, ENVT and ENVN shown for ‖kH‖ ≤ 20. ENVT is the limit to MPI resolution, and defines MPI bandwidth [14]. ENVN has approximately half the intrinsic resolution with FWHMT ...

We note that ENVT in (9) was first derived in temporal frequency space in Rahmer et al. [12], and is derived in x-space in our earlier work [14]. The second envelope in (10) is unique to the generalized x-space formulation, and gives the resolution of the transverse component of the point spread function perpendicular to the FFP velocity vector.

The cubic relationship between resolution and particle diameter is absolutely critical. This relationship can be understood by looking at the origin of MPI’s signal. MPI’s resolution relies on the nonlinear effect of a small applied magnetic field causing a SPIO nanoparticle tracer to magnetically saturate. The Langevin equation tells us that the field required to saturate a single magnetic nanoparticle decreases with the nanoparticle’s volume. As a result, resolution improves with the cube of the magnetic nanoparticle diameter.

A. MPI 3-D Point Spread Function

The MPI process generates signals in multiple axes. While we typically will build inductive receiver coils that are perpendicular to the physical axes (x, y, and z) of the instrument, the instrument produces images on an internal reference frame formed by vectors collinear and transverse to the FFP velocity vector, equation M32. The collinear and transverse images are distinct from the tangential and normal components.

To understand the origin of the collinear and transverse components of the PSF, we will examine the vector components of h(x). Supposing that the velocity vector is aligned with the x unit vector, i.e., equation M33. Then, the collinear component is

equation M34

and the transverse components are

equation M35

where we have arbitrarily picked two perpendicular unit axes corresponding to the y and z axes, ê2 and ê3. The resulting components of the PSF are shown in Fig. 3. The collinear and transverse components of the PSF form an excellent basis set for image reconstruction.

Fig. 3
Collinear and transverse components of the matrix point spread function. The received images rotate with vector equation M46 (see Fig. 4). The collinear PSF component peak amplitude is 370% the tangential PSF component peak amplitude. The area of the box drawn in ...

However, if the FFP velocity vector is not oriented with one of the cardinal directions of the instrument, the received images change with the orientation of the velocity vector equation M36. To illustrate how the reference frame is oriented with an arbitrary FFP path, we see the collinear and transverse components of the PSF oriented with the velocity vector in Fig. 4.

Fig. 4
MPI images are acquired on a reference frame formed by vectors collinear and transverse to the velocity vector equation M47.

While the equations presented here are general, it is interesting to look at the point spread function in an algebraic equation. If we fix the excitation vector ê1 along the z axis and assume a diagonal gradient matrix G = diag(Gxx, Gyy, Gzz), then we can write the collinear PSF

equation M37

and one of the transverse PSFs on the receive axis aligned with the x axis

equation M38

where equation M39. The PSFs for these equations are shown in Fig. 3.

The collinear component is similar to the real part of a Lorentzian function seen in NMR, and is an even function. The collinear component is desirable and forms the bulk of the resolution and signal of the MPI imaging process. The collinear component is a vector sum of both the tangential and normal envelopes, ENVT and ENVN, with the sharper envelope aligned with the velocity vector.

The transverse component is similar to the dispersion or odd-valued spectral component in NMR and is an odd function. The transverse component is a vector difference of the two point spread function envelopes. Across the diagonal (see Fig. 3), the signal received is precisely PSF[perpendicular] = ±(1/2)(L(H/Hsat) − L(H/Hsat)/‖H/Hsat‖). As a result, the transverse component is significantly smaller in magnitude than the collinear component.

IV. Methods

In this section we begin by discussing two methods crucial to x-space MPI, gridding, and fundamental recovery. We use gridding use to transform the received signal from the time domain to the image domain. We then introduce the technique of a partial FOV, which we use to estimate the low frequency content that is lost when we filter out the fundamental frequency. These x-space reconstruction techniques do not require regularization, optimization techniques, or prior knowledge of the magnetic response of the tracer. We conclude by describing the construction of a small-scale MPI imager to test the x-space theory presented in this paper.

A. Reconstruction Methods: Gridding

Gridding in MPI is simply the process of sampling the received signal s(t) onto a grid in real space, or x-space, that corresponds to the instantaneous location of the FFP. In our formulation, we grid the collinear and transverse components separately. To do this, we first separate the image into collinear and transverse signals

equation M40

We have chosen an arbitrary unit vector ê1 to cross with the velocity unit vector equation M41 to build a perpendicular basis set of transverse vectors. Choice of this arbitrary vector assumes that ê1 and equation M42 are not collinear, i.e., equation M43.

If we design our pulse sequence so that the velocity unit vector equation M44 is constant, e.g., with fast movement only in one direction, we are able to simplify gridding of the collinear and transverse signals. Ignoring the receiver coil sensitivity and choosing arbitrary unit vectors ê1 and ê2 not collinear with the FFP velocity equation M45 gives us the

equation image

where we normalize to the magnitude of the FFP velocity [14]. Similar gridding can be done for the remaining transverse images. This image equation is akin to the k-space analysis of MRI [18], [19], [27], but here the scanning occurs in x-space rather than in k-space, so no Fourier transform is required.

B. Reconstruction Methods: Fundamental Recovery

As mentioned above, a great challenge in MPI is the loss of the fundamental frequency. MPI interrogates magnetic nanoparticles by subjecting the sample to a rapidly varying magnetic field. This applied magnetic field contaminates the received signal as the applied field induces a signal in the receive coil that is many orders of magnitude larger than the signal generated by the magnetic nanoparticles. This applied field is typically a sinusoid [1], [14], [16], whose frequency we term the “fundamental frequency” [7], [14].

Partial FOV

To overcome the problems associated with the loss of the fundamental frequency, we first introduce the concept of a partial FOV. This concept is similar to the focus field described by Weizenecker et al. [2]. MPI systems can use a large gradient field to increase resolution at the expense of reducing the FOV of a scan. For example, the scanner described in this paper generates a 30 mTpeak−peak excitation amplitude on top of a 6 T/m gradient, giving a total FOV of about half a centimeter. This excitation amplitude already exceeds the limits of magnetostimulation for a chest scanner and is nearing the magnetostimulation limit for an extremity scanner [14]. However, we can take partial FOV images that we stitch together for the full FOV by slowly moving the average position of the FFP mechanically or with an electromagnet.

Loss and Recovery of the Fundamental Signal

When scanning an image with overlapping partial FOVs, it is possible to recover the lost fundamental signal, which is important to shift invariance. We approach the problem of the lost fundamental signal by considering high-pass filtering of the time-domain signal as a loss of low-spatial frequency information. For the loss of temporal frequencies near the fundamental frequency, we can approximate this spatial-signal loss as a dc offset. Surprisingly, this means that if we acquire multiple overlapping partial FOVs, we can sequentially find the overlap between signals that minimizes their overlap error. Since only a constant dc offset was lost and assuming boundary conditions at the endpoints of the scan, the resulting reassembled image will be an excellent approximation of the original spatial convolution. A full derivation of the dc offset shift and noise properties of various reconstruction results are beyond the scope of this paper.

C. Experimental Methods

To test the principles described in this paper, we built a 3-D MPI scanner as shown in Fig. 5. The system is constructed with permanent magnet gradient (6 T/m down the bore and 3 T/m transverse to the bore) and an excitation coil collinear to the bore. The FFP is rapidly scanned using the resonant transmit coil and the signal produced is received with a receive coil wound collinear to the transmit coil. The receive coil receives the collinear component of the vector PSF. We note that the transmit and receive coils are collinear with the larger gradient along the bore, which is twice the magnitude of the gradient transverse to the bore. The collinearity of the coils was chosen for ease of construction but results in an intrinsic resolution in the transverse direction that is approximately four times worse than in the collinear direction [see (11)].

Fig. 5
X-space MPI imager. (a) Tomographic MPI scanner with 2 cm × 2 cm × 4 cm FOV. The excitation transmit coil generates a 30 mT peak-to-peak oscillating magnetic field at 20 kHz. The NdFeB magnet gradient generates a gradient of 6 T/m down ...

The resonant excitation coil generates 30 mT peak-to-peak at 20 kHz and is driven by an audio amplifier (AE Techron LVC5050, Elkhart, IN) with ~5 kW of instantaneous power at a pulsed 2% duty cycle. The signal from the receive coil is filtered by a passive notch filter, amplified by a battery powered preamplifier (SR560, Stanford Research Systems, Sunnyvale, CA), and high-pass filtered at 35 kHz (SIM965, Stanford Research Systems, Sunnyvale, CA). Following the analog signal chain, the signal is digitized by a 16-bit data acquisition system with a 1.25 MSPS sampling rate (National Instruments USB-6259, Austin, TX), phase corrected, and low-pass filtered at 400 kHz. The system is controlled by custom software written in MATLAB (Mathworks MATLAB, Natick, MA).

The 30 mT peak-to-peak excitation enables a partial FOV of approximately 0.5 cm along the z axis. The signal for the received partial FOV is gridded to the instantaneous location of the FFP and assigned to a physical location on the phantom. The phantom is stepped in 1 mm steps along the z axis for 4 cm, acquiring a partial FOV line scan at each step. The line scans are reassembled as described in Section IV-B by estimating the missing fundamental to generate an assembled full FOV of 4.5 cm along the z axis. A total of 20 line scans are taken by moving in 1 mm steps transverse to the bore, for a full FOV of 2 cm in the y axis. Phantoms are constructed using 400 µm ID tubing filled with undiluted SPIO tracer (Resovist, Bayer-Schering, Berlin, Germany).

V. Results

In Fig. 6 we see experimental data showing the 1-D scan of a point source before and after fundamental recovery. We recover the fundamental by estimating the dc offset of each segment so that we find a maximally smooth image.

Fig. 6
[Top] Experimental data showing 40 overlapping partial FOV line-scans for a 400 µm point source phantom. The baseline component for each partial FOV is lost in the scanning process due to the contamination of first harmonic imaging data by direct ...

In Figs. 7 and and99 we see images measured with the x-space MPI imager. As seen in the images of the PSF in Fig. 7 and line scans in Fig. 8, the FWHM in the normal axis is 4.6 times wider than the FWHM along the axis of the imager. The lower resolution in the normal axis agrees with the theoretical prediction in (11).We attribute the measured FWHM being wider than the theoretical prediction to the nanoparticle behaving differently than in our model, and the phantom being a line source rather than a point source. Our future scanners will orient the magnetic fields differently so as to optimize the shape of the PSF. We believe the “CAL” phantom image shown in Fig. 9(b) is the first native MPI image without any sharpening or deconvolution and with full recovery of the fundamental frequency, which is crucial for maintaining the LSI properties of the system. As seen in the figure, the received image for the “CAL” phantom faithfully represents the phantom.

Fig. 7
(a) Measured two-dimensional collinear PSF showing excellent correspondence to Fig. 3. The measured FWHM is 1.6 mm along the imager bore and 7.4 mm transverse to the imager bore. The PSF phantom is a 400 µm tubing oriented perpendicular to the ...
Fig. 8
Profiles across the point spread function shown in Fig. 7(a), (b) show good agreement between theoretical and measured values. [TOP] Line scan down the bore. [BOTTOM] Line scan perpendicular to the imager bore.
Fig. 9
(a) “CAL” Phantom built using 400 µm ID tubing filled with undiluted tracer and encapsulated. (b) Intrinsic MPI image of the CAL phantom showing excellent correspondence to the phantom image. FOV: 4 cm × 2 cm, Pixel size: ...

VI. Discussion

The x-space technique that we describe in this paper is a new look at the MPI imaging process. We started from three hypotheses, that the gradient creates a single FFP, the adiabatic Langevin model, and that the loss of the low frequencies recoverable. These three hypotheses give us a powerful framework to analyze the MPI imaging process. The adiabatic Langevin model has become standard in the MPI world, and this is certainly valid for smaller particles at reasonable scanning speeds. We proved the uniqueness of the FFP. And we presented experimental evidence to show that lost low frequency information is recoverable. And, again, bear in mind that the loss of first harmonic information is a universal MPI challenge and not unique to x-space scanning.

Unlike the approach in Rahmer et al. [12], x-space theory does not require a repeating sinusoidal excitation or specific Lissajous pulse sequence. The x-space formulation also motivates image reconstruction technique that is robust, scalable, and significantly faster than inversion of the system matrix and it does not require precharacterization of the magnetic nanoparticles or system [8]. Numerical inversion of the system matrix may take significant time for larger system matrices such as a 128 × 128 × 128 image since matrix inversion does not scale linearly. The computation required for x-space reconstruction, however, requires only a scaling and gridding. For example, our current reconstruction code reconstructs the received signal faster than our analog-to-digital converter is able to digitize data.

The intrinsic FWHM resolution predicted by our first work on space analysis [14] agrees with the nondeconvolved resolution limit predicted by Rahmer et al. [12]. However, the 2-D and 3-D analysis presented in this paper extends these initial analyses to show that the intrinsic resolution changes with the orientation the FFP movement sequence. In Rahmer et al. [12] the author states for their realtime mouse imager [2] with a 5.5 T/m gradient and Resovist tracer that “the observed resolution was not better than Δx ≈ 1.5 mm.” Their gradient strength was smaller than the gradient described in this paper, their tracer is identical, and their achieved resolution approximately the same our result. Rahmer et al. [12] attribute the discrepancy from their expected 500 µm resolution given their previous deconvolved results [1], [3] due to “the wide distribution of particle sizes and the regularization applied in reconstruction to mitigate limited SNR.” We believe x-space theory shows that their system admirably reached the physical resolution limit of their imager given a physiologically relevant SNR. Indeed, our experimental results are in perfect alignment with Rahmer’s (18) that gives the undeconvolved system resolution.

One of the advantages of the x-space approach is that it allows the separation of reconstruction from deconvolution. Not only is this modular, but it allows for both fast reconstruction and fast deconvolution using the wealth of standard deconvolution methods [28], including fast Fourier transform filters and Wiener deconvolution. For example, while not shown in this paper, in practice we lightly apply Wiener deconvolution to the intrinsic x-space MPI image in order to visually improve image contrast and apparent resolution. This is especially relevant for images with high SNR as deconvolution can increase the image resolution at the expense of SNR. Unfortunately, more aggressive deconvolution results in significant noise amplification [15] and often in image artifacts [28].

An important counter-example to aggressive deconvolution is the example of clinical PET imaging, whose resolution is worse than 4 mm resolution in a modern scanner [29]. Clearly this would be a prime candidate for standard deconvolution methods. However, physicians have been reluctant to use deconvolution with PET clinically, presumably for reasons of robustness and SNR loss.

Fortunately, we believe that deconvolution may be considered unnecessary in the future, given the strong dependence of intrinsic resolution on nanoparticle diameter. Our results here demonstrate the accuracy of the x-space resolution (11). Consequently, we believe it is reasonable to assume that resolution could improve with the cube of the particle diameter so long as the nano-particles remain superparamagnetic and continue to satisfy the assumptions necessary for x-space MPI. It is clear that increasing the apparent diameter of the SPIO magnetic core to 25 or even 30 nm and dealing with relaxation effects will be crucial to improving the intrinsic resolution of the image.

We believe that our work proves that MPI can be approximated as LSI when we recover the low-frequency information. We note that most imaging modalities are not strictly LSI, but require a minor approximation to simplify. For example, all computed-tomography imaging scanners reconstruct in Hounsfield units, which are the natural log of the received signal. That is to say, without taking the natural log, CT is not strictly an LSI system. Existing works [12], [17] aim to analyze the MPI process in mathematical terms, but they do not endeavour to prove linearity or shift invariance. In this paper, we believe we have proved that MPI is indeed a LSI system with our three hypotheses, and that our experimental results show that the recovery of the first harmonic enables experimental MPI systems to be accurately modeled as LSI. Once it is shown that x-space can produce images comparable in quality to those reconstructed using system matrices, x-space analysis could be accepted as a powerful analytical tool.

VII. Conclusion

We have presented and tested a theory for the magnetic particle imaging process in x-space. We see that MPI can be described as a LSI system with a well behaved point spread function. Our primary theoretical conclusions are seen in (8) and (12) which give the MPI signal and MPI image equation. We then built a scanner that we successfully used to acquire a 2-D intrinsic MPI image and an image of a point spread function. Only a few hypotheses are required for x-space theory to accurately predict the experimental PSF in 2-D. All three of these hypotheses are either well accepted (adiabatic alignment) or were proved mathematically (unique FFP location) or experimentally (loss of first harmonic information is recoverable). This new analysis of MPI scanning provides powerful insights for optimizing MPI scanner and scanning sequences. It has also prompted a significant advance in image reconstruction computation time.

Acknowledgment

The authors would like to thank K. Lu and B. Zheng for their assistance, and G. Lee for excellent discussions.

This work was supported in part by the California Institute for Regenerative Medicine (CIRM) graduate fellowship under Training Grant T1-00007, in part by the Siebel Foundation, in part by CIRM Tools and Technology Grant RT1-01055-1, in part by a University of California Berkeley Bioengineering graduate fellowship, and in part by the National Institutes of Health training grant. The contents of this publication are solely the responsibility of the authors and do not necessarily represent the official views of CIRM or any other agency of the State of California.

Footnotes

Color versions of one or more of the figures in this paper are available online at http://ieeexplore.ieee.org.

Contributor Information

Patrick W. Goodwill, Department of Bioengineering, University of California, Berkeley, CA 94720 USA.

Steven M. Conolly, Department of Bioengineering, University of California, Berkeley, CA 94720 USA.

References

1. Gleich B, Weizenecker J. Tomographic imaging using the nonlinear response of magnetic particles. Nature. 2005 Jun.435:1214–1217. [PubMed]
2. Weizenecker J, Gleich B, Rahmer J, Dahnke H, Borgert J. Three-dimensional real-time in vivo magnetic particle imaging. Phys. Med. Biol. 2009 Feb.54:L1–L10. [PubMed]
3. Gleich B, Weizenecker J, Borgert J. Experimental results on fast 2d-encoded magnetic particle imaging. Phys. Med. Biol. 2008 Feb.53:N81–N84. [PubMed]
4. Weizenecker J, Borgert J, Gleich B. A simulation study on the resolution and sensitivity of magnetic particle imaging. Phys. Med. Biol. 2007 Nov.52:6363–6374. [PubMed]
5. Weizenecker J, Gleich B, Borgert J. Magnetic particle imaging using a field free line. J. Phys. D: Appl. Phys. 2008 May;41:105009.
6. Sattel TF, Knopp T, Biederer S, Gleich B, Weizenecker J, Borgert J, Buzug TM. Single-sided device for magnetic particle imaging. J. Phys. D: Appl. Phys. 2008 Dec.42 022001.
7. Biederer S, Knopp T, Sattel T, Lüdtke-Buzug K, Gleich B, Weizenecker J, Borgert J, Buzug T. Magnetization response spectroscopy of superparamagnetic nanoparticles for MPI. J. Phys. D: Appl. Phys. 2009;42:205007.
8. Knopp T, Sattel TF, Biederer S, Rahmer J, Weizenecker J, Gleich B, Borgert J, Buzug TM. Model-based reconstruction for magnetic particle imaging. IEEE Trans. Med. Imag. 2010 Jan.29(no. 1):12–18. [PubMed]
9. Knopp T, Biederer S, Sattel T, Rahmer J, Weizenecker J, Gleich B, Borgert J, Buzug T. 2d model-based reconstruction for magnetic particle imaging. Med. Phys. 2010;37:485. [PubMed]
10. Ferguson RM, Minard KR, Krishnan KM. Optimization of nanoparticle core size for magnetic particle imaging. J. Magn. Magn. Mater. 2009 Jan.321:1548–1551. [PMC free article] [PubMed]
11. Prince J, Links J. Medical Imaging Signals and Systems. Upper Saddle River, NJ: Prentice Hall; 2005.
12. Rahmer J, Weizenecker J, Gleich B, Borgert J. Signal encoding in magnetic particle imaging: Properties of the system function. BMC Med. Imag. 2009;9(no. 1):4. [PMC free article] [PubMed]
13. Knopp T, Biederer S, Sattel T, Weizenecker J, Gleich B, Borgert J, Buzug TM. Trajectory analysis for magnetic particle imaging. Phys. Med. Biol. 2008 Dec.54:385–397. [PubMed]
14. Goodwill PW, Conolly SM. The x-space formulation of the magnetic particle imaging process: One-dimensional signal, resolution, bandwidth, SNR, SAR, and magnetostimulation. IEEE Trans. Med. Imaging. 2010 Nov.29(no. 11):1851–1859. [PubMed]
15. Shahram M, Milanfar P. Imaging below the diffraction limit: A statistical analysis. IEEE Trans. Image Process. 2004 May;13(no. 5):677–689. [PubMed]
16. Goodwill P, Scott G, Stang P, Conolly S. Narrowband magnetic particle imaging. IEEE Trans. Med. Imag. 2009 Aug.28(no. 8):1231–1237. [PubMed]
17. Schomberg H. Magnetic particle imaging: Model and reconstruction. IEEE Int. Symp. Biomed. Imag. Int. Symp. Conf. Rec. (ISBI’10) 2010 Mar.:992–995.
18. Twieg DB. The k-trajectory formulation of the NMR imaging process with applications in analysis and synthesis of imaging methods. Med. Phys. 1983 Jan.10:610–621. [PubMed]
19. Ljunggren S. A simple graphical representation of Fourier-based imaging methods. J. Magn. Reson. 1983;54(no. 2):338–343.
20. Bernstein M, Zhou X, Polzin J, King K, Ganin A, Pelc N, Glover G. Concomitant gradient terms in phase contrast MR: Analysis and correction. Magn. Reson. Med. 1998;39:300–308. [PubMed]
21. Jackson JD. Classical Electrodynamics. New York: Wiley; 1975.
22. Myers W, Mößle M, Clarke J. Correction of concomitant gradient artifacts in experimental microtesla MRI. J. Magn. Reson. 2005;177(no. 2):274–284. [PubMed]
23. Rosensweig R. Ferrohydrodynamics. Cambridge, U.K.: Cambridge Univ. Press; 1985.
24. Shliomis M. Magnetic fluids. Physics-Uspekhi. 1974;17(no. 2):153–169.
25. Fannin P, Charles S. The study of a ferrofluid exhibiting both Brownian and Neel relaxation. J. Phys. D: Appl. Phys. 1989 Jan.
26. Hoult D, Richards R. The signal-to-noise ratio of the nuclear magnetic resonance experiment. J. Magn. Reson. 1976 Jan.24:71–85. [PubMed]
27. Brown TR, Kincaid BM, Ugurbil K. NMR chemical shift imaging in three dimensions. Proc. Nat. Acad. Sci. USA. 1982 Jun.79:3523–3526. [PubMed]
28. González RC, Woods RE. Digital Image Processing. Upper Saddle River, NJ: Prentice Hall; 2008.
29. Brambilla M, Secco C, Dominietto M, Matheoud R, Sacchetti G, Inglese E. Performance characteristics obtained for a new 3-dimensional lutetium oxyorthosilicate-based whole-body PET/CT scanner with the national electrical manufacturers association nu 2–2001 standard. J. Nucl. Med. 2005 Dec.46:2083. [PubMed]