Institute of Physics and Engineering in Medicine Phys. Med. Biol. 60 (2015) 4033–4044
Physics in Medicine & Biology doi:10.1088/0031-9155/60/10/4033
Symmetries of the 2D magnetic particle imaging system matrix A Weber1 and T Knopp2,3 1
Bruker Biospin MRI, Rudolf-Plank-Straße 23, 76275 Ettlingen, Germany Section for Biomedical Imaging, University Medical Center Hamburg-Eppendorf, Martinistraße 52, 20246 Hamburg, Germany 3 Institute for Biomedical Imaging, Hamburg University of Technology, Schwarzenbergstraße 95, 21073 Hamburg, Germany 2
E-mail: [email protected] and [email protected] Received 19 February 2015, revised 20 March 2015 Accepted for publication 1 April 2015 Published 28 April 2015 Abstract
In magnetic particle imaging (MPI), the relation between the particle distribution and the measurement signal can be described by a linear system of equations. For 1D imaging, it can be shown that the system matrix can be expressed as a product of a convolution matrix and a Chebyshev transformation matrix. For multidimensional imaging, the structure of the MPI system matrix is not yet fully explored as the sampling trajectory complicates the physical model. It has been experimentally found that the MPI system matrix rows have symmetries and look similar to the tensor products of Chebyshev polynomials. In this work we will mathematically prove that the 2D MPI system matrix has symmetries that can be used for matrix compression. Keywords: magnetic particle imaging, system matrix, symmetry (Some figures may appear in colour only in the online journal) 1. Introduction Magnetic particle imaging is a promising imaging modality that allows the determination of the spatial distribution of super-paramagnetic nanoparticles (SPIOs) (Ferguson et al 2009, Pankhurst et al 2009, Krishnan 2010) at a unique combination of sensitivity, acquisition time, and spatial resolution (Gleich and Weizenecker 2005, Weizenecker et al 2009, Knopp and Buzug 2012). In contrast to magnetic resonance imaging (MRI), the particle concentration is directly quantifiable and the contrast is positive, i.e. no background signal from the human body is able to distort the nanoparticle signal (Weaver et al 2008, Goodwill and Conolly 2010, Goodwill and Conolly 2011). Currently, MPI is entering the preclinical research landscape 0031-9155/15/104033+12$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
4033
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
with the availability of the first commercially developed MPI scanner that features an imaging bore of 12 cm (Bruker Biospin MRI GmbH, Ettlingen). The MPI technique uses various static and dynamic magnetic fields during the imaging process. The so-called selection field has a field-free point (FFP) at a single point in space that acts as sensitive spot where the main portion of the signal is generated. For 1D imaging, the FFP is shifted back and forth along a line with sinusoidal velocity using the so-called drive field. Considering a simplified model, Rahmer et al have shown that the 1D MPI system matrix can be expressed as a product of a convolution matrix and a Chebyshev transformation matrix (Rahmer et al 2009). However, for 2D imaging, where the FFP travels along a Lissajous curve, no mathematical relation to Chebyshev polynomials has been derived so far, although qualitatively a similarity can be observed (Rahmer et al 2009, Knopp et al 2010a). Due to the insufficient understanding of the 2D system matrix, but also as the simplified model does not account for non-negligible relaxation effects, multidimensional MPI system matrices are usually measured instead of computed. To this end a small delta sample filled with SPIOs is shifted through the field-of-view (FOV) while measuring the system response at each image position. In this way, the system matrix is measured column by column using a robot. While this procedure is accurate, it is very time consuming and lasts from several hours up to days. For this reason there have been different attempts to reduce the amount of calibration scans either using a model-based approach (Knopp et al 2010a, 2010c) or using compressed sensing (Knopp and Weber 2013). The purpose of this work is to study the symmetry properties of the MPI system matrix. In Rahmer et al (2009) and Rahmer et al (2012) it has already been qualitatively observed that the 2D and 3D MPI system matrices have symmetries. In this work we will show mathematically that the 2D system matrix has two symmetries that can be used to reduce the amount of calibration scans that have to be acquired when measuring the MPI system matrix. 2. Theory 2.1. Signal generation and encoding
In 2D MPI the signal is generated by applying the so-called drive field ⎛ HxD(t ) ⎞ ⎛ Ax sin(2πfx t ) ⎞ ⎟=⎜ ⎟ HD(t ) = ⎜⎜ ⎟ ⎜ A sin(2πf t ) ⎟ D ( ) H t y y ⎠ ⎝ y ⎠ ⎝
to the magnetic nanoparticles. The drive field is homogeneous in space so that all particles exhibit the same excitation. The particles respond with a magnetization M (t) and the change in magnetization can be detected using electromagnetic coils. For technical reasons the drive field excitation is sinusoidal in all three directions. The resulting curve of the dynamic drive field is a so-called Lissajous curve (Papula 2009). To get a periodic drive field the excitation frequencies f x and f y have to exhibit a rational frequency ratio
fx fy
=
Ky Kx
with the frequency ratio
parameters Kx, Ky ∈ N. Defining the basis frequency f Base ≔ f x Kx = f y Ky the period of the drive field is given by T =
lcm(Kx, Ky ) fBase
, i.e. the least common multiplier of the frequency ratio
parameters divided by the basis frequency. For spatial encoding, in addition to the drive field, a gradient field
4034
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
H S (r ) = G r
⎛Gx with G = ⎜⎜ ⎝0
0⎞ ⎟⎟ Gy ⎠
is applied. This field is called a selection field and assigns to every spatial position r = (rx, ry)⊤ a certain offset. In turn, particles at different positions respond with a different particle magnetization so that the induced signals can be distinguished. In total, the applied magnetic field is given by H (r, t ) = HS (r ) + HD(t ),
i.e. the superposition of the drive and the selection field. 2.2. The MPI system function
In MPI the relation between the measurement signal and the particle concentration is described by the so-called system function. When using a 2D imaging sequence the particle magnetization is detected using L = 2 receive coils. Hereby the k-th Fourier coefficient of the measurement signal of receive channel l is given by
∫Ω
ˆl , k = u(1) sˆl, k (r )c(r ) d3r ,
which is the integral of the particle distribution multiplied with the k-th system function component of receive channel l over the sensitive area Ω. In the following the index k will be referred to as a frequency index as both the induced voltage and the system function are usually considered in frequency space. Each system function component includes information about the particle properties, the scanner characteristic, as well as the measurement sequence, and is given by T ∂ kt μ ˆl, k : RL → C, sˆl, k (r ) = −aˆl, k 0 s(2) m (r, t ) ⋅ pl (r )e−2π i T dt , T 0 ∂t
∫
where aˆl, k denotes the transfer function of the receive chain, m (r, t ) is the mean magnetic moment caused by the magnetic particles, and pl(r) describes the coil sensitivity depending on the spatial position r. If isotropic super-paramagnetic particles with instantaneous relaxation are assumed, the mean magnetic moment can be described by H (r , t ) . m (r , t ) = m ( ‖ H (r , t ) ‖ ) (3) ‖H (r, t )‖
To derive symmetry properties of the system function it is assumed that all coils are aligned orthogonal to each other, such that the coil sensitivity is equal to the unit vector and spatially independent, i.e. pl(r) = el. Furthermore, the transfer function is neglected for the symmetry derivation, because it is independent of the spatial position. Making these assumptions yields kt ∂ ~l, k : RL → C, m ~l, k (r ) ≔ sˆl, k (r ) = − μ0 ml (r, t )e−2π i T dt . m (4) aˆl, k T 0 ∂t
∫
T
~l, k as the signal function in the following since it is independent of the the transfer We refer to m function of the receive chain and thus only contains the particle signal characteristics. 4035
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
2.3. 1D Signal function symmetries
Before investigating the 2D case, we consider a 1D imaging sequence where only one drive field is active, i.e. without loss in generality Ax ≠ 0 and Ay = 0. The following theorem describes the symmetry of the 1D signal function: Theorem 1. Under the assumption of an ideal coil configuration and isotropic super-paramagnetic particles with instantaneous relaxation, it holds for all 1D signal function components that
~ k+1 ~ m mx, k (rx ). x, k (−rx ) = (−1)
(5)
Proof. The 1D signal function components can be described as a convolution of the derivative
of the magnetization response with Chebyshev polynomials of the second kind (Rahmer et al 2009). Since these Chebyshev polynomials are even functions for odd k and odd functions for even k, and since the symmetry is preserved under the convolution for a symmetric kernel, the symmetry relation (5) holds. □ 2.4. 2D Signal function symmetries
In the 2D case no closed expression of the signal function has been derived so far. Thus, the symmetry properties of the 1D signal function can not be directly transferred to the 2D case. In this work the symmetry properties of the signal function components are derived from the symmetric behavior of the drive field, which is given in lemma 1. ⎛ HxD(t ) ⎞ ⎟ be an ideal 2D drive field with even Kx, odd Ky and period T. Then for Lemma 1. Let ⎜⎜ ⎟ D ⎝ Hy (t ) ⎠ all τ ∈ R it holds that ⎛ ⎞ ⎞ ⎛ D⎛ T ⎞⎞ D⎛ T ⎛ HxD(τ )⎞ ⎜ Hx ⎜⎝ 2 − τ ⎠⎟ ⎟ ⎜−Hx ⎝⎜ 2 + τ ⎠⎟ ⎟ ⎛−HxD(T − τ )⎞ ⎟ ⎜ ⎟ ⎜ ⎟. ⎜ ⎟ ⎜ =⎜ ⎟ ⎟=⎜ D ⎟ ⎜ D ⎟=⎜ ⎝ Hy (τ )⎠ ⎜−H D⎛⎜ T − τ ⎞⎟ ⎟ ⎜ H D⎛⎜ T + τ ⎞⎟ ⎟ ⎝−Hy (T − τ )⎠ ⎜ y ⎟ ⎜ y ⎟ ⎝2 ⎠⎠ ⎝ ⎝2 ⎠⎠ ⎝
Proof. The first equation will be exemplarily proven. The other two equations can be derived
with the same technique. Inserting the period T =
KxKy fBase
and the excitation frequencies fi =
fBase , Ki
i = x, y yields
⎛ ⎛ fBase KxKy ⎞⎞ ⎛ ⎛ ⎛T ⎞⎞ ⎞ ⎜ ⎞⎞ ⎛ T D⎛ sin 2 2 f A − π π τ ⎜ ⎟⎟ ⎜ ⎟ ⎜ ⎟ x ⎜ ⎟ ⎜ ⎟ sin 2 A f − π τ x x ⎜ Hx ⎝ − τ ⎠ ⎟ x⎝ ⎠ K 2 f ⎜ ⎝ ⎠⎟ ⎝ ⎠ x 2 Base 2 ⎟ ⎜ ⎟ ⎜ ⎜ ⎟. = = ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎛ ⎞ ⎞ ⎛ ⎞ ⎛ T ⎛ ⎞ T K K f ⎜ ⎟ ⎜⎜−HyD⎜ − τ ⎟ ⎟⎟ ⎜−Ay sin⎜2πf ⎜ − τ ⎟⎟ ⎟ ⎜−A sin⎜2π Base x y − 2πf τ ⎟ ⎟ y ⎝2 ⎠⎠ ⎝ ⎟ ⎠⎠ ⎠ ⎜ y ⎝ ⎝ ⎝ y⎝ 2 Ky 2fBase ⎠⎠ ⎝
By cancelling and applying the trigonometric sum identity sin (α ± β) = sin (α) cos (β) ± cos (α) sin (β), this term can be reformulated to 4036
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
⎛ ⎛T ⎞⎞ ⎜ HxD⎜⎝ − τ ⎠⎟ ⎟ ⎛ A sin(−2πf τ ) cos(πK ) + A cos(−2πf τ ) sin(πK ) ⎞ y x y 2 x x ⎜ ⎟ ⎜ x ⎟. =⎜ ⎟ ⎜ ⎟ − A sin( − 2 f ) cos( π K ) − A cos( − 2 π f τ ) sin( π K ) π τ ⎛ ⎞ T y x y x y y ⎠ ⎜⎜−HyD⎜ − τ ⎟ ⎟⎟ ⎝ ⎝2 ⎠⎠ ⎝
Since for all k ∈ Z it holds that sin (kπ) = 0 and cos (kπ) = (−1)k it follows ⎛ ⎛T ⎞⎞ ⎜ HxD⎜⎝ − τ ⎠⎟ ⎟ ⎛−A sin(−2πf τ ) ⎞ 2 x ⎜ ⎟ ⎜ x ⎟. =⎜ ⎜ ⎟ −Ay sin(−2πfy τ ) ⎟⎠ ⎛ ⎞ T ⎝ D ⎜⎜−Hy ⎜ − τ ⎟ ⎟⎟ ⎝2 ⎠⎠ ⎝
Hence the first equation ⎛ ⎛T ⎞⎞ ⎜ HxD⎜⎝ − τ ⎠⎟ ⎟ ⎛ D ⎞ 2 ⎜ ⎟ ⎜ Hx (τ ) ⎟ ⎜ ⎟=⎜ D ⎟ ⎜⎜−HyD⎛⎜ T − τ ⎞⎟ ⎟⎟ ⎝ Hy (τ ) ⎠ ⎝2 ⎠⎠ ⎝
is proven, since the sine is an odd function. □ In figure 1 the symmetric behaviour of the drive field is visualized for an exemplary fre3
quency ratio of 4 , i.e. Kx = 4 and Ky = 3. Based on the drive field symmetries summarized in lemma 1 we can derive symmetry properties for the signal function components. Theorem 2. Under the assumption of an ideal coil configuration, isotropic super-paramagnetic particles with instantaneous relaxation, even Kx and odd Ky, it holds for all 2D signal function components that
⎛⎛−rx ⎞⎞ ⎛⎛ rx ⎞⎞ k+1 ~ ~ mx, k ⎜⎝⎜⎝ ry ⎟⎠⎟⎠ = (−1) mx, k ⎜⎝⎜⎝ ry ⎟⎠⎟⎠,
k ∈ N0
(6)
⎛⎛−r ⎞⎞ ⎛⎛ r ⎞⎞ ~ y, k ⎜⎜ x ⎟⎟ = (−1)k m ~ y, k ⎜⎜ x ⎟⎟, k ∈ N0 m (7) r ⎝⎝ y ⎠⎠ ⎝⎝ ry ⎠⎠ ⎛⎛ r ⎞⎞ ⎛⎛ r ⎞⎞ ~ x, k ⎜⎜ x ⎟⎟ = (−1)k + 1 m ~ x, k ⎜⎜ x ⎟⎟, k ∈ N0 m (8) ⎝⎝−ry ⎠⎠ ⎝⎝ ry ⎠⎠ ⎛⎛ r ⎞⎞ ⎛⎛ r ⎞⎞ ~ y, k ⎜⎜ x ⎟⎟ = (−1)k m ~ y, k ⎜⎜ x ⎟⎟, k ∈ N0 m (9) ⎝⎝−ry ⎠⎠ ⎝⎝ ry ⎠⎠ ⎛−rx ⎞ ⎟ is introduced. ⎝ ry ⎠
r ≔⎜ Proof. To simplify the notation the mirrored position ~
According to lemma 1 the relation
4037
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
Figure 1. Visualization of the course of the 2D Lissajous curve with a frequency ratio of
3 and identical amplitudes A1 = A2 = 1. On the left, the complete Lissajous trajectory is 4 shown. The sections of the Lissajous trajectory which belong to a quarter of the period are presented in the middle. Here, the circle marks the start point and the cross marks the end point of each section. On the right the x- and y-components are visualized as a function of time.
⎛ ⎛ T ⎞⎞ ⎞ D⎛ ⎛−Gxrx + HxD(t )⎞ ⎜−⎜⎝Gxrx + Hx ⎜⎝t + 2 ⎠⎟⎟⎠ ⎟ ⎟ ⎟=⎜ H (~r , t ) = ⎜⎜ ⎟ ⎜ ⎟ D ⎛ T⎞ ⎝ Gyry + Hy (t ) ⎠ ⎜ ⎜ Gyry + HyD⎜t + ⎟ ⎟⎟ ⎝ ⎝ 2⎠ ⎠
holds for all t ∈ R. Thus the norm of the magnetic field at position r and ~r matches at certain time points ⎛ T⎞ ‖H (~r , t ) ‖= H ⎜r, t + ⎟ . ⎝ 2⎠
Inserting this relation into the equation of the mean magnetic moment (3) yields the symmetry ⎛ T D ⎜ − Gxrx + Hx t + 2 ⎛∂ ⎜ ~ ⎞ T ⎜ m x (r , t ) ⎟ H r, t + 2 ⎜ ~ ⎞ ⎛ ∂ t ⎛ ⎞ ( , ) T ∂ t ∂ H r ⎜ ⎟ = m (‖H (~r , t )‖) ⎜ ⎜ ⎟ ~r , t ) = ∂t m ⎜⎝ H ⎝r , t + 2 ⎠ ⎟⎠⎜ ⎜∂ ⎟ ∂t ( H T D ⎜ m y(~r , t ) ⎟ ⎜ Gyry + Hy t + 2 ⎝ ∂t ⎠ ⎜ T ⎜ H r, t + 2 ⎝
(
(
( (
⎛ ∂ ⎛ T ⎞⎞ ⎜− mx ⎜⎝r, t + ⎟⎠ ⎟ 2 ⎟ ⎜ ∂t =⎜ ⎟ ⎜⎜ ∂ m y⎜⎛r, t + T ⎟⎞ ⎟⎟ ⎝ ∂t ⎝ 2⎠ ⎠
)
(
)
)) ⎞⎟
)
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
which holds for ‖H(~r, t )‖ ≠ 0. Due to the continuity of the mean magnetic moment this relation also holds for ‖H (~r , t )‖= 0. Therefore the signal function components can be represented by 4038
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
~ x, k (~r ) = − μ0 m T μ0 =− T μ ~ y, k (~r ) = − 0 m T μ =− 0 T
∫0
T
∫0
T
∫0
T
∫0
T
kt ∂ mx (~r , t )e−2π i T dt ∂t kt T⎞ ∂ ⎛ − mx ⎜r , t + ⎟e−2π i T dt ⎝ ⎠ ∂t 2 kt ∂ m y(~r , t )e−2π i T dt ∂t kt ∂ ⎜⎛ T⎞ m y r, t + ⎟e−2π i T dt . ∂t ⎝ 2⎠
T
kT 2
Since a shift of 2 in the time domain implies a multiplication with e−2π i T = (−1)k in the frequency domain, the equations (6) and (7) are proven. ⎛ rx ⎞ Applying the same techniques to ~r ≔ ⎜ ⎟, the symmetry ⎝−ry ⎠ ⎛ ∂ ⎛ T ⎞⎞ ⎛∂ ~ ⎞ ⎜ mx (r , t ) ⎟ ⎜− ∂t mx ⎜⎝r, 2 − t ⎠⎟ ⎟ ⎟ ⎜ ∂t ⎟=⎜ ⎟ ⎜∂ ⎟ ⎜ ∂ ⎛ T ⎞ ⎜ m y(~r , t ) ⎟ ⎜⎜ m y⎜r, − t ⎟ ⎟⎟ ⎝ ∂t ⎠ ⎝ ∂t ⎝ 2 ⎠⎠ T
can be derived for all t ∈ R. In addition to the shift of 2 a mirroring occurs in the time domain and thus additionally to the multiplication with (−1)k a complex conjugation has to be performed. Taking this into account proves equations (8) and (9). □ 2.5. System function symmetries
In the previous sections we have derived symmetry properties for the MPI signal function ~l, k (r ). However, in practice one measures the system function m ~ l, k (r ), sˆl, k (r ) = aˆl, k m which includes the transfer function of the receive chain aˆl, k. This complex valued function is not known and complicates the application of the symmetry rules. For those symmetry equations where no complex conjugation is involved (e.g. in x-direction for 2D imaging) the symmetry can be directly exploited even when aˆl, k is unknown. In the other case, the complex conjugation cannot be directly applied as a complex conjugation of the system function sˆl, k (r ) ~l, k (r ) but also the unknown transfer function aˆl, k. does not only conjugate the signal function m In order to still apply a symmetry rule we first exploit that a complex conjugation can be written as a multiplication with e−2iϕ, where ϕ is the phase of the complex number. Then we take into account that the phase of the signal function is approximately constant for every spatial position such that only one global symmetry phase ϕ has to be determined. To this end, we use a simulated signal function in this work. 2.6. Mirror axis
In practice it is challenging to place the system function positions in such a way that the magnetic field center is exactly matched. In turn, a simple mirroring would lead to discontinuities at the mirror axis and even small deviations can have a negative influence on the 4039
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
reconstruction result. To avoid this discontinuous mirroring we use an overlap and weight the mirrored data with the overlap data depending on the distance to mirroring axis. 3. Materials and methods In order to validate our theoretical findings we use experimental 2D MPI data and show that the mirroring of the MPI system matrix can be exploited to reduce the calibration time for determining the MPI system matrix considerably. The considered MPI data has already been used in several MPI publications where the reconstruction of MPI data was investigated (Knopp et al 2010a, 2010b, Knopp and Weber 2013). The data has been acquired with the MPI mouse scanner that was first published in Weizenecker et al (2009). The scanner has a bore of 32 mm, a selection field gradient of 5.5 Tm −1μ0−1, in x-direction and 2.75 Tm −1μ0−1 in y-direction. Drive-field amplitudes were chosen to be 18 mTμ0−1 in x- and y-direction. Thus, the measuring field covered by the drive field is about 13.09 × 6.54 mm2. The sampling trajectory is a 2D Lissajous trajectory with a fre32
quency ratio of 33 ( f x = 2.5 MHz/96 ≈ 26 041.7 kHz and f y = 2.5 MHz/99 ≈ 25 252.5 kHz). The voltage signals induced in the receive coils are captured at a sampling rate of 20 MS s−1, whereas only frequencies up to 1 MHz are considered for reconstruction as higher frequencies only contained noise. The total number of frequency components in each receive channel is thus 1268. A system matrix was acquired at 68 × 37 equidistant positions spanning a FOV of 20.4 × 12.0 mm2. When applying the symmetry rules only half or a quarter of the data was used. Phantom data was acquired using a P phantom that was rotated during the dynamic experiment. The phantom consists of 12 holes, each with a diameter of 0.5 mm. For the actual image reconstruction we used the iterative Kaczmarz method (Knopp et al 2010b) considering 5 iterations and a manually tuned Tikonov regularization parameter. We note that Kaczmarz method has proven its robustness in several MPI publications including the first in-vivo publication (Weizenecker et al 2009). 4. Results 4.1. System matrix
In figure 2 the absolute value of three selected system function components of the measured and mirrored system matrices are shown. As one can see, the system matrix is highly symmetric. When applying the horizontal mirroring one can hardly identify a difference. In contrast, in the vertical direction one can identify a non-symmetric structure on the right-hand edge of the system matrix. Consequently, the mirrored system function components differ in these boundary regions. Note that the effect might be due to a field inhomogeneity that we have not taken into account in this work. In figure 3 the phase plots of the previously shown system function components are illustrated. As one can see, the phase is indeed piecewise constant and only flips the sign when the real/imaginary parts have a zero-crossing. Next we consider the deviation that is introduced in the system matrix when applying different mirroring strategies. Since the original fully sampled system matrix is available we can calculate the normalized root mean square error (NRMSE) between the original and the mirrored system matrix. Please note that we have reordered the system function components 4040
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
vertical
horizontal
vertical & horizontal
387 kHz
237 kHz
128 kHz
original
Figure 2. Several original system function components and their mirrored versions at frequencies f ∈{128 kHz, 237 kHz, 387 kHz}. Vertical mirroring is performed from left to right and horizontal mirroring from up to down. The black and white dashed line marks the axis of mirroring. The scan parameters are described in section 3.
128 kHz
237 kHz
387 kHz
Figure 3. Phase plot of several original system function components at frequencies f ∈{128 kHz, 237 kHz, 387 kHz}.
so that they are ordered descending according to their signal-to-noise ratio (SNR). Therewith the new frequency index κ = 1 corresponds to the system function component exhibiting the highest SNR. In figure 4 the NRMSE between the original and the mirrored system matrices is shown. The error is calculated for all matrix rows independently and the error is plotted as a function of the row index. As one can see, the horizontal mirroring introduces the smallest error even though in contrast to vertical mirroring a systematic error is made due to the phase estimation. The error introduced by vertical mirroring is slightly higher, which can be explained by the n on-symmetric boundary artifacts that could be identified in figure 2. The error when mirroring the system 4041
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
horizontal mirroring
vertical mirroring
vertical & horizontal mirroring
Figure 4. NRMSE between the original and the mirrored system matrix for vertical, horizontal, and vertical & horizontal mirroring. The NRMSE is calculated for each system function component independently and shown as a function of the frequency index κ, which includes the sortation of the system function components according to their SNR.
matrix in both directions is dominated by the vertical mirror error. In general the mirror error increases with decreasing SNR. 4.2. Phantom measurements
The reconstruction results using the phantom data and different system matrices are shown in figure 5. Two different frames of the time series are shown. As one can see, the results when 4042
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
Horizontal
Vertical
Both
2nd frame
1st frame
Full
Figure 5. Reconstruction results of the P phantom at two different time points. The data was reconstructed using the original system matrix and three different mirrored system matrices (exploiting vertical, horizontal, and both symmetries).
using the fully sampled system matrix and when using the mirrored system matrices are comparable in resolution and quality. The quality does not degrade when applying both symmetries. 5. Discussion In this work the symmetry of the magnetic particle imaging system matrix has been investigated. We have shown that the 1D system matrix has a mirror symmetry, whereas the 2D system matrix also has a mirror symmetry in both axes, but one also has to take into account a complex conjugation in one direction. The symmetries can be used to significantly shorten the calibration method that is required to measure the MPI system matrix. We have shown that in 2D only 33%of the data has to be acquired, while the remaining parts can be determined by applying the mirror rules. The 33%are composed of 25%for a quarter of the system function component and 8%overlap to ensure a smooth mirroring. The quality of the reconstructed phantom data did not degrade when using the mirrored system matrices. On the measured data we have identified non-symmetric parts at the boundary of the FOV. These may be due to a field inhomogeneity and a non-symmetric placement of the receive coils. We therefore plan to refine our model and take into account the non-symmetric sensitivity pattern of the receive coil in future work. The same mathematical proof technique as in 2D can be used to establish symmetry properties for the 3D MPI system matrix. However, in contrast to the 2D case, where a mirror symmetry in two dimensions is fulfilled, in 3D no symmetry in all excitation directions can be found. For example, selecting the frequency ratio parameters Kx doubly even, Ky odd, and Kz even, the signal function components only exhibit a symmetry with respect to the yz-plane and, in addition, a point symmetry with respect to the center of the yz-planes. The proof for the 3D case together with experimental validation on measured 3D data will be content of future work. References Ferguson R, Minard K and Krishnan K 2009 Optimization of nanoparticle core size for magnetic particle imaging J. Magn. Magn. Mater. 321 1548–51 Gleich B and Weizenecker J 2005 Tomographic imaging using the nonlinear response of magnetic particles Nature 435 1214–7 4043
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
Goodwill P and Conolly S 2010 The x-space formulation of the magnetic particle imaging process: onedimensional signal, resolution, bandwidth, SNR, SAR, and magnetostimulation IEEE Trans. Med. Imaging 29 1851–9 Goodwill P and Conolly S 2011 Multi-dimensional x-space magnetic particle imaging IEEE Trans. Med. Imaging 30 1581–90 Knopp T, Biederer S, Sattel T, Rahmer J, Weizenecker J, Gleich B, Borgert J and Buzug T 2010 2D model-based reconstruction for magnetic particle imaging Med. Phys. 37 485–91 Knopp T and Buzug T M 2012 Magnetic Particle Imaging: an Introduction to Imaging Principles and Scanner Instrumentation (Berlin: Springer) (www.springer.com/medicine/radiology/ book/978-3-642-04198-3) Knopp T, Rahmer J, Sattel T, Biederer S, Weizenecker J, Gleich B, Borgert J and Buzug T 2010 Weighted iterative reconstruction for magnetic particle imaging Phys. Med. Biol. 55 1577–89 Knopp T, Sattel T, Biederer S, Rahmer J, Weizenecker J, Gleich B, Borgert J and Buzug T 2010 Modelbased reconstruction for magnetic particle imaging IEEE Trans. Med. Imaging 29 12–8 Knopp T and Weber A 2013 Sparse reconstruction of the magnetic particle imaging system matrix IEEE Trans. Med. Imaging 32 1473–80 Krishnan K 2010 Biomedical nanomagnetics: a spin through possibilities in imaging, diagnostics, and therapy IEEE Trans. Magn. 46 2523–58 Pankhurst Q, Thanh N, Jones S and Dobson J 2009 Progress in applications of magnetic nanoparticles in biomedicine J. Phys. D: Appl. Phys. 42 224001 Papula L 2009 Mathematik für Ingenieure und Naturwissenschaftler Band 1 (Wiesbaden: Vieweg + Teubner) Rahmer J, Weizenecker J, Gleich B and Borgert J 2009 Signal encoding in magnetic particle imaging BMC Med. Imaging 9 4 Rahmer J, Weizenecker J, Gleich B and Borgert J 2012 Analysis of a 3-d system function measured for magnetic particle imaging IEEE Trans. Med. Imaging 31 1289–99 Weaver J, Rauwerdink A, Sullivan C and Baker I 2008 Frequency distribution of the nanoparticle magnetization in the presence of a static as well as a harmonic magnetic field Med. Phys. 35 1988–94 Weizenecker J, Gleich B, Rahmer J, Dahnke H and Borgert J 2009 Three-dimensional real-time in vivo magnetic particle imaging Phys. Med. Biol. 54 L1–10
4044
Copyright of Physics in Medicine & Biology is the property of IOP Publishing and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
Physics in Medicine & Biology doi:10.1088/0031-9155/60/10/4033
Symmetries of the 2D magnetic particle imaging system matrix A Weber1 and T Knopp2,3 1
Bruker Biospin MRI, Rudolf-Plank-Straße 23, 76275 Ettlingen, Germany Section for Biomedical Imaging, University Medical Center Hamburg-Eppendorf, Martinistraße 52, 20246 Hamburg, Germany 3 Institute for Biomedical Imaging, Hamburg University of Technology, Schwarzenbergstraße 95, 21073 Hamburg, Germany 2
E-mail: [email protected] and [email protected] Received 19 February 2015, revised 20 March 2015 Accepted for publication 1 April 2015 Published 28 April 2015 Abstract
In magnetic particle imaging (MPI), the relation between the particle distribution and the measurement signal can be described by a linear system of equations. For 1D imaging, it can be shown that the system matrix can be expressed as a product of a convolution matrix and a Chebyshev transformation matrix. For multidimensional imaging, the structure of the MPI system matrix is not yet fully explored as the sampling trajectory complicates the physical model. It has been experimentally found that the MPI system matrix rows have symmetries and look similar to the tensor products of Chebyshev polynomials. In this work we will mathematically prove that the 2D MPI system matrix has symmetries that can be used for matrix compression. Keywords: magnetic particle imaging, system matrix, symmetry (Some figures may appear in colour only in the online journal) 1. Introduction Magnetic particle imaging is a promising imaging modality that allows the determination of the spatial distribution of super-paramagnetic nanoparticles (SPIOs) (Ferguson et al 2009, Pankhurst et al 2009, Krishnan 2010) at a unique combination of sensitivity, acquisition time, and spatial resolution (Gleich and Weizenecker 2005, Weizenecker et al 2009, Knopp and Buzug 2012). In contrast to magnetic resonance imaging (MRI), the particle concentration is directly quantifiable and the contrast is positive, i.e. no background signal from the human body is able to distort the nanoparticle signal (Weaver et al 2008, Goodwill and Conolly 2010, Goodwill and Conolly 2011). Currently, MPI is entering the preclinical research landscape 0031-9155/15/104033+12$33.00 © 2015 Institute of Physics and Engineering in Medicine Printed in the UK
4033
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
with the availability of the first commercially developed MPI scanner that features an imaging bore of 12 cm (Bruker Biospin MRI GmbH, Ettlingen). The MPI technique uses various static and dynamic magnetic fields during the imaging process. The so-called selection field has a field-free point (FFP) at a single point in space that acts as sensitive spot where the main portion of the signal is generated. For 1D imaging, the FFP is shifted back and forth along a line with sinusoidal velocity using the so-called drive field. Considering a simplified model, Rahmer et al have shown that the 1D MPI system matrix can be expressed as a product of a convolution matrix and a Chebyshev transformation matrix (Rahmer et al 2009). However, for 2D imaging, where the FFP travels along a Lissajous curve, no mathematical relation to Chebyshev polynomials has been derived so far, although qualitatively a similarity can be observed (Rahmer et al 2009, Knopp et al 2010a). Due to the insufficient understanding of the 2D system matrix, but also as the simplified model does not account for non-negligible relaxation effects, multidimensional MPI system matrices are usually measured instead of computed. To this end a small delta sample filled with SPIOs is shifted through the field-of-view (FOV) while measuring the system response at each image position. In this way, the system matrix is measured column by column using a robot. While this procedure is accurate, it is very time consuming and lasts from several hours up to days. For this reason there have been different attempts to reduce the amount of calibration scans either using a model-based approach (Knopp et al 2010a, 2010c) or using compressed sensing (Knopp and Weber 2013). The purpose of this work is to study the symmetry properties of the MPI system matrix. In Rahmer et al (2009) and Rahmer et al (2012) it has already been qualitatively observed that the 2D and 3D MPI system matrices have symmetries. In this work we will show mathematically that the 2D system matrix has two symmetries that can be used to reduce the amount of calibration scans that have to be acquired when measuring the MPI system matrix. 2. Theory 2.1. Signal generation and encoding
In 2D MPI the signal is generated by applying the so-called drive field ⎛ HxD(t ) ⎞ ⎛ Ax sin(2πfx t ) ⎞ ⎟=⎜ ⎟ HD(t ) = ⎜⎜ ⎟ ⎜ A sin(2πf t ) ⎟ D ( ) H t y y ⎠ ⎝ y ⎠ ⎝
to the magnetic nanoparticles. The drive field is homogeneous in space so that all particles exhibit the same excitation. The particles respond with a magnetization M (t) and the change in magnetization can be detected using electromagnetic coils. For technical reasons the drive field excitation is sinusoidal in all three directions. The resulting curve of the dynamic drive field is a so-called Lissajous curve (Papula 2009). To get a periodic drive field the excitation frequencies f x and f y have to exhibit a rational frequency ratio
fx fy
=
Ky Kx
with the frequency ratio
parameters Kx, Ky ∈ N. Defining the basis frequency f Base ≔ f x Kx = f y Ky the period of the drive field is given by T =
lcm(Kx, Ky ) fBase
, i.e. the least common multiplier of the frequency ratio
parameters divided by the basis frequency. For spatial encoding, in addition to the drive field, a gradient field
4034
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
H S (r ) = G r
⎛Gx with G = ⎜⎜ ⎝0
0⎞ ⎟⎟ Gy ⎠
is applied. This field is called a selection field and assigns to every spatial position r = (rx, ry)⊤ a certain offset. In turn, particles at different positions respond with a different particle magnetization so that the induced signals can be distinguished. In total, the applied magnetic field is given by H (r, t ) = HS (r ) + HD(t ),
i.e. the superposition of the drive and the selection field. 2.2. The MPI system function
In MPI the relation between the measurement signal and the particle concentration is described by the so-called system function. When using a 2D imaging sequence the particle magnetization is detected using L = 2 receive coils. Hereby the k-th Fourier coefficient of the measurement signal of receive channel l is given by
∫Ω
ˆl , k = u(1) sˆl, k (r )c(r ) d3r ,
which is the integral of the particle distribution multiplied with the k-th system function component of receive channel l over the sensitive area Ω. In the following the index k will be referred to as a frequency index as both the induced voltage and the system function are usually considered in frequency space. Each system function component includes information about the particle properties, the scanner characteristic, as well as the measurement sequence, and is given by T ∂ kt μ ˆl, k : RL → C, sˆl, k (r ) = −aˆl, k 0 s(2) m (r, t ) ⋅ pl (r )e−2π i T dt , T 0 ∂t
∫
where aˆl, k denotes the transfer function of the receive chain, m (r, t ) is the mean magnetic moment caused by the magnetic particles, and pl(r) describes the coil sensitivity depending on the spatial position r. If isotropic super-paramagnetic particles with instantaneous relaxation are assumed, the mean magnetic moment can be described by H (r , t ) . m (r , t ) = m ( ‖ H (r , t ) ‖ ) (3) ‖H (r, t )‖
To derive symmetry properties of the system function it is assumed that all coils are aligned orthogonal to each other, such that the coil sensitivity is equal to the unit vector and spatially independent, i.e. pl(r) = el. Furthermore, the transfer function is neglected for the symmetry derivation, because it is independent of the spatial position. Making these assumptions yields kt ∂ ~l, k : RL → C, m ~l, k (r ) ≔ sˆl, k (r ) = − μ0 ml (r, t )e−2π i T dt . m (4) aˆl, k T 0 ∂t
∫
T
~l, k as the signal function in the following since it is independent of the the transfer We refer to m function of the receive chain and thus only contains the particle signal characteristics. 4035
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
2.3. 1D Signal function symmetries
Before investigating the 2D case, we consider a 1D imaging sequence where only one drive field is active, i.e. without loss in generality Ax ≠ 0 and Ay = 0. The following theorem describes the symmetry of the 1D signal function: Theorem 1. Under the assumption of an ideal coil configuration and isotropic super-paramagnetic particles with instantaneous relaxation, it holds for all 1D signal function components that
~ k+1 ~ m mx, k (rx ). x, k (−rx ) = (−1)
(5)
Proof. The 1D signal function components can be described as a convolution of the derivative
of the magnetization response with Chebyshev polynomials of the second kind (Rahmer et al 2009). Since these Chebyshev polynomials are even functions for odd k and odd functions for even k, and since the symmetry is preserved under the convolution for a symmetric kernel, the symmetry relation (5) holds. □ 2.4. 2D Signal function symmetries
In the 2D case no closed expression of the signal function has been derived so far. Thus, the symmetry properties of the 1D signal function can not be directly transferred to the 2D case. In this work the symmetry properties of the signal function components are derived from the symmetric behavior of the drive field, which is given in lemma 1. ⎛ HxD(t ) ⎞ ⎟ be an ideal 2D drive field with even Kx, odd Ky and period T. Then for Lemma 1. Let ⎜⎜ ⎟ D ⎝ Hy (t ) ⎠ all τ ∈ R it holds that ⎛ ⎞ ⎞ ⎛ D⎛ T ⎞⎞ D⎛ T ⎛ HxD(τ )⎞ ⎜ Hx ⎜⎝ 2 − τ ⎠⎟ ⎟ ⎜−Hx ⎝⎜ 2 + τ ⎠⎟ ⎟ ⎛−HxD(T − τ )⎞ ⎟ ⎜ ⎟ ⎜ ⎟. ⎜ ⎟ ⎜ =⎜ ⎟ ⎟=⎜ D ⎟ ⎜ D ⎟=⎜ ⎝ Hy (τ )⎠ ⎜−H D⎛⎜ T − τ ⎞⎟ ⎟ ⎜ H D⎛⎜ T + τ ⎞⎟ ⎟ ⎝−Hy (T − τ )⎠ ⎜ y ⎟ ⎜ y ⎟ ⎝2 ⎠⎠ ⎝ ⎝2 ⎠⎠ ⎝
Proof. The first equation will be exemplarily proven. The other two equations can be derived
with the same technique. Inserting the period T =
KxKy fBase
and the excitation frequencies fi =
fBase , Ki
i = x, y yields
⎛ ⎛ fBase KxKy ⎞⎞ ⎛ ⎛ ⎛T ⎞⎞ ⎞ ⎜ ⎞⎞ ⎛ T D⎛ sin 2 2 f A − π π τ ⎜ ⎟⎟ ⎜ ⎟ ⎜ ⎟ x ⎜ ⎟ ⎜ ⎟ sin 2 A f − π τ x x ⎜ Hx ⎝ − τ ⎠ ⎟ x⎝ ⎠ K 2 f ⎜ ⎝ ⎠⎟ ⎝ ⎠ x 2 Base 2 ⎟ ⎜ ⎟ ⎜ ⎜ ⎟. = = ⎜ ⎟ ⎜ ⎟ ⎜ ⎟ ⎛ ⎞ ⎞ ⎛ ⎞ ⎛ T ⎛ ⎞ T K K f ⎜ ⎟ ⎜⎜−HyD⎜ − τ ⎟ ⎟⎟ ⎜−Ay sin⎜2πf ⎜ − τ ⎟⎟ ⎟ ⎜−A sin⎜2π Base x y − 2πf τ ⎟ ⎟ y ⎝2 ⎠⎠ ⎝ ⎟ ⎠⎠ ⎠ ⎜ y ⎝ ⎝ ⎝ y⎝ 2 Ky 2fBase ⎠⎠ ⎝
By cancelling and applying the trigonometric sum identity sin (α ± β) = sin (α) cos (β) ± cos (α) sin (β), this term can be reformulated to 4036
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
⎛ ⎛T ⎞⎞ ⎜ HxD⎜⎝ − τ ⎠⎟ ⎟ ⎛ A sin(−2πf τ ) cos(πK ) + A cos(−2πf τ ) sin(πK ) ⎞ y x y 2 x x ⎜ ⎟ ⎜ x ⎟. =⎜ ⎟ ⎜ ⎟ − A sin( − 2 f ) cos( π K ) − A cos( − 2 π f τ ) sin( π K ) π τ ⎛ ⎞ T y x y x y y ⎠ ⎜⎜−HyD⎜ − τ ⎟ ⎟⎟ ⎝ ⎝2 ⎠⎠ ⎝
Since for all k ∈ Z it holds that sin (kπ) = 0 and cos (kπ) = (−1)k it follows ⎛ ⎛T ⎞⎞ ⎜ HxD⎜⎝ − τ ⎠⎟ ⎟ ⎛−A sin(−2πf τ ) ⎞ 2 x ⎜ ⎟ ⎜ x ⎟. =⎜ ⎜ ⎟ −Ay sin(−2πfy τ ) ⎟⎠ ⎛ ⎞ T ⎝ D ⎜⎜−Hy ⎜ − τ ⎟ ⎟⎟ ⎝2 ⎠⎠ ⎝
Hence the first equation ⎛ ⎛T ⎞⎞ ⎜ HxD⎜⎝ − τ ⎠⎟ ⎟ ⎛ D ⎞ 2 ⎜ ⎟ ⎜ Hx (τ ) ⎟ ⎜ ⎟=⎜ D ⎟ ⎜⎜−HyD⎛⎜ T − τ ⎞⎟ ⎟⎟ ⎝ Hy (τ ) ⎠ ⎝2 ⎠⎠ ⎝
is proven, since the sine is an odd function. □ In figure 1 the symmetric behaviour of the drive field is visualized for an exemplary fre3
quency ratio of 4 , i.e. Kx = 4 and Ky = 3. Based on the drive field symmetries summarized in lemma 1 we can derive symmetry properties for the signal function components. Theorem 2. Under the assumption of an ideal coil configuration, isotropic super-paramagnetic particles with instantaneous relaxation, even Kx and odd Ky, it holds for all 2D signal function components that
⎛⎛−rx ⎞⎞ ⎛⎛ rx ⎞⎞ k+1 ~ ~ mx, k ⎜⎝⎜⎝ ry ⎟⎠⎟⎠ = (−1) mx, k ⎜⎝⎜⎝ ry ⎟⎠⎟⎠,
k ∈ N0
(6)
⎛⎛−r ⎞⎞ ⎛⎛ r ⎞⎞ ~ y, k ⎜⎜ x ⎟⎟ = (−1)k m ~ y, k ⎜⎜ x ⎟⎟, k ∈ N0 m (7) r ⎝⎝ y ⎠⎠ ⎝⎝ ry ⎠⎠ ⎛⎛ r ⎞⎞ ⎛⎛ r ⎞⎞ ~ x, k ⎜⎜ x ⎟⎟ = (−1)k + 1 m ~ x, k ⎜⎜ x ⎟⎟, k ∈ N0 m (8) ⎝⎝−ry ⎠⎠ ⎝⎝ ry ⎠⎠ ⎛⎛ r ⎞⎞ ⎛⎛ r ⎞⎞ ~ y, k ⎜⎜ x ⎟⎟ = (−1)k m ~ y, k ⎜⎜ x ⎟⎟, k ∈ N0 m (9) ⎝⎝−ry ⎠⎠ ⎝⎝ ry ⎠⎠ ⎛−rx ⎞ ⎟ is introduced. ⎝ ry ⎠
r ≔⎜ Proof. To simplify the notation the mirrored position ~
According to lemma 1 the relation
4037
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
Figure 1. Visualization of the course of the 2D Lissajous curve with a frequency ratio of
3 and identical amplitudes A1 = A2 = 1. On the left, the complete Lissajous trajectory is 4 shown. The sections of the Lissajous trajectory which belong to a quarter of the period are presented in the middle. Here, the circle marks the start point and the cross marks the end point of each section. On the right the x- and y-components are visualized as a function of time.
⎛ ⎛ T ⎞⎞ ⎞ D⎛ ⎛−Gxrx + HxD(t )⎞ ⎜−⎜⎝Gxrx + Hx ⎜⎝t + 2 ⎠⎟⎟⎠ ⎟ ⎟ ⎟=⎜ H (~r , t ) = ⎜⎜ ⎟ ⎜ ⎟ D ⎛ T⎞ ⎝ Gyry + Hy (t ) ⎠ ⎜ ⎜ Gyry + HyD⎜t + ⎟ ⎟⎟ ⎝ ⎝ 2⎠ ⎠
holds for all t ∈ R. Thus the norm of the magnetic field at position r and ~r matches at certain time points ⎛ T⎞ ‖H (~r , t ) ‖= H ⎜r, t + ⎟ . ⎝ 2⎠
Inserting this relation into the equation of the mean magnetic moment (3) yields the symmetry ⎛ T D ⎜ − Gxrx + Hx t + 2 ⎛∂ ⎜ ~ ⎞ T ⎜ m x (r , t ) ⎟ H r, t + 2 ⎜ ~ ⎞ ⎛ ∂ t ⎛ ⎞ ( , ) T ∂ t ∂ H r ⎜ ⎟ = m (‖H (~r , t )‖) ⎜ ⎜ ⎟ ~r , t ) = ∂t m ⎜⎝ H ⎝r , t + 2 ⎠ ⎟⎠⎜ ⎜∂ ⎟ ∂t ( H T D ⎜ m y(~r , t ) ⎟ ⎜ Gyry + Hy t + 2 ⎝ ∂t ⎠ ⎜ T ⎜ H r, t + 2 ⎝
(
(
( (
⎛ ∂ ⎛ T ⎞⎞ ⎜− mx ⎜⎝r, t + ⎟⎠ ⎟ 2 ⎟ ⎜ ∂t =⎜ ⎟ ⎜⎜ ∂ m y⎜⎛r, t + T ⎟⎞ ⎟⎟ ⎝ ∂t ⎝ 2⎠ ⎠
)
(
)
)) ⎞⎟
)
⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠
which holds for ‖H(~r, t )‖ ≠ 0. Due to the continuity of the mean magnetic moment this relation also holds for ‖H (~r , t )‖= 0. Therefore the signal function components can be represented by 4038
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
~ x, k (~r ) = − μ0 m T μ0 =− T μ ~ y, k (~r ) = − 0 m T μ =− 0 T
∫0
T
∫0
T
∫0
T
∫0
T
kt ∂ mx (~r , t )e−2π i T dt ∂t kt T⎞ ∂ ⎛ − mx ⎜r , t + ⎟e−2π i T dt ⎝ ⎠ ∂t 2 kt ∂ m y(~r , t )e−2π i T dt ∂t kt ∂ ⎜⎛ T⎞ m y r, t + ⎟e−2π i T dt . ∂t ⎝ 2⎠
T
kT 2
Since a shift of 2 in the time domain implies a multiplication with e−2π i T = (−1)k in the frequency domain, the equations (6) and (7) are proven. ⎛ rx ⎞ Applying the same techniques to ~r ≔ ⎜ ⎟, the symmetry ⎝−ry ⎠ ⎛ ∂ ⎛ T ⎞⎞ ⎛∂ ~ ⎞ ⎜ mx (r , t ) ⎟ ⎜− ∂t mx ⎜⎝r, 2 − t ⎠⎟ ⎟ ⎟ ⎜ ∂t ⎟=⎜ ⎟ ⎜∂ ⎟ ⎜ ∂ ⎛ T ⎞ ⎜ m y(~r , t ) ⎟ ⎜⎜ m y⎜r, − t ⎟ ⎟⎟ ⎝ ∂t ⎠ ⎝ ∂t ⎝ 2 ⎠⎠ T
can be derived for all t ∈ R. In addition to the shift of 2 a mirroring occurs in the time domain and thus additionally to the multiplication with (−1)k a complex conjugation has to be performed. Taking this into account proves equations (8) and (9). □ 2.5. System function symmetries
In the previous sections we have derived symmetry properties for the MPI signal function ~l, k (r ). However, in practice one measures the system function m ~ l, k (r ), sˆl, k (r ) = aˆl, k m which includes the transfer function of the receive chain aˆl, k. This complex valued function is not known and complicates the application of the symmetry rules. For those symmetry equations where no complex conjugation is involved (e.g. in x-direction for 2D imaging) the symmetry can be directly exploited even when aˆl, k is unknown. In the other case, the complex conjugation cannot be directly applied as a complex conjugation of the system function sˆl, k (r ) ~l, k (r ) but also the unknown transfer function aˆl, k. does not only conjugate the signal function m In order to still apply a symmetry rule we first exploit that a complex conjugation can be written as a multiplication with e−2iϕ, where ϕ is the phase of the complex number. Then we take into account that the phase of the signal function is approximately constant for every spatial position such that only one global symmetry phase ϕ has to be determined. To this end, we use a simulated signal function in this work. 2.6. Mirror axis
In practice it is challenging to place the system function positions in such a way that the magnetic field center is exactly matched. In turn, a simple mirroring would lead to discontinuities at the mirror axis and even small deviations can have a negative influence on the 4039
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
reconstruction result. To avoid this discontinuous mirroring we use an overlap and weight the mirrored data with the overlap data depending on the distance to mirroring axis. 3. Materials and methods In order to validate our theoretical findings we use experimental 2D MPI data and show that the mirroring of the MPI system matrix can be exploited to reduce the calibration time for determining the MPI system matrix considerably. The considered MPI data has already been used in several MPI publications where the reconstruction of MPI data was investigated (Knopp et al 2010a, 2010b, Knopp and Weber 2013). The data has been acquired with the MPI mouse scanner that was first published in Weizenecker et al (2009). The scanner has a bore of 32 mm, a selection field gradient of 5.5 Tm −1μ0−1, in x-direction and 2.75 Tm −1μ0−1 in y-direction. Drive-field amplitudes were chosen to be 18 mTμ0−1 in x- and y-direction. Thus, the measuring field covered by the drive field is about 13.09 × 6.54 mm2. The sampling trajectory is a 2D Lissajous trajectory with a fre32
quency ratio of 33 ( f x = 2.5 MHz/96 ≈ 26 041.7 kHz and f y = 2.5 MHz/99 ≈ 25 252.5 kHz). The voltage signals induced in the receive coils are captured at a sampling rate of 20 MS s−1, whereas only frequencies up to 1 MHz are considered for reconstruction as higher frequencies only contained noise. The total number of frequency components in each receive channel is thus 1268. A system matrix was acquired at 68 × 37 equidistant positions spanning a FOV of 20.4 × 12.0 mm2. When applying the symmetry rules only half or a quarter of the data was used. Phantom data was acquired using a P phantom that was rotated during the dynamic experiment. The phantom consists of 12 holes, each with a diameter of 0.5 mm. For the actual image reconstruction we used the iterative Kaczmarz method (Knopp et al 2010b) considering 5 iterations and a manually tuned Tikonov regularization parameter. We note that Kaczmarz method has proven its robustness in several MPI publications including the first in-vivo publication (Weizenecker et al 2009). 4. Results 4.1. System matrix
In figure 2 the absolute value of three selected system function components of the measured and mirrored system matrices are shown. As one can see, the system matrix is highly symmetric. When applying the horizontal mirroring one can hardly identify a difference. In contrast, in the vertical direction one can identify a non-symmetric structure on the right-hand edge of the system matrix. Consequently, the mirrored system function components differ in these boundary regions. Note that the effect might be due to a field inhomogeneity that we have not taken into account in this work. In figure 3 the phase plots of the previously shown system function components are illustrated. As one can see, the phase is indeed piecewise constant and only flips the sign when the real/imaginary parts have a zero-crossing. Next we consider the deviation that is introduced in the system matrix when applying different mirroring strategies. Since the original fully sampled system matrix is available we can calculate the normalized root mean square error (NRMSE) between the original and the mirrored system matrix. Please note that we have reordered the system function components 4040
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
vertical
horizontal
vertical & horizontal
387 kHz
237 kHz
128 kHz
original
Figure 2. Several original system function components and their mirrored versions at frequencies f ∈{128 kHz, 237 kHz, 387 kHz}. Vertical mirroring is performed from left to right and horizontal mirroring from up to down. The black and white dashed line marks the axis of mirroring. The scan parameters are described in section 3.
128 kHz
237 kHz
387 kHz
Figure 3. Phase plot of several original system function components at frequencies f ∈{128 kHz, 237 kHz, 387 kHz}.
so that they are ordered descending according to their signal-to-noise ratio (SNR). Therewith the new frequency index κ = 1 corresponds to the system function component exhibiting the highest SNR. In figure 4 the NRMSE between the original and the mirrored system matrices is shown. The error is calculated for all matrix rows independently and the error is plotted as a function of the row index. As one can see, the horizontal mirroring introduces the smallest error even though in contrast to vertical mirroring a systematic error is made due to the phase estimation. The error introduced by vertical mirroring is slightly higher, which can be explained by the n on-symmetric boundary artifacts that could be identified in figure 2. The error when mirroring the system 4041
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
horizontal mirroring
vertical mirroring
vertical & horizontal mirroring
Figure 4. NRMSE between the original and the mirrored system matrix for vertical, horizontal, and vertical & horizontal mirroring. The NRMSE is calculated for each system function component independently and shown as a function of the frequency index κ, which includes the sortation of the system function components according to their SNR.
matrix in both directions is dominated by the vertical mirror error. In general the mirror error increases with decreasing SNR. 4.2. Phantom measurements
The reconstruction results using the phantom data and different system matrices are shown in figure 5. Two different frames of the time series are shown. As one can see, the results when 4042
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
Horizontal
Vertical
Both
2nd frame
1st frame
Full
Figure 5. Reconstruction results of the P phantom at two different time points. The data was reconstructed using the original system matrix and three different mirrored system matrices (exploiting vertical, horizontal, and both symmetries).
using the fully sampled system matrix and when using the mirrored system matrices are comparable in resolution and quality. The quality does not degrade when applying both symmetries. 5. Discussion In this work the symmetry of the magnetic particle imaging system matrix has been investigated. We have shown that the 1D system matrix has a mirror symmetry, whereas the 2D system matrix also has a mirror symmetry in both axes, but one also has to take into account a complex conjugation in one direction. The symmetries can be used to significantly shorten the calibration method that is required to measure the MPI system matrix. We have shown that in 2D only 33%of the data has to be acquired, while the remaining parts can be determined by applying the mirror rules. The 33%are composed of 25%for a quarter of the system function component and 8%overlap to ensure a smooth mirroring. The quality of the reconstructed phantom data did not degrade when using the mirrored system matrices. On the measured data we have identified non-symmetric parts at the boundary of the FOV. These may be due to a field inhomogeneity and a non-symmetric placement of the receive coils. We therefore plan to refine our model and take into account the non-symmetric sensitivity pattern of the receive coil in future work. The same mathematical proof technique as in 2D can be used to establish symmetry properties for the 3D MPI system matrix. However, in contrast to the 2D case, where a mirror symmetry in two dimensions is fulfilled, in 3D no symmetry in all excitation directions can be found. For example, selecting the frequency ratio parameters Kx doubly even, Ky odd, and Kz even, the signal function components only exhibit a symmetry with respect to the yz-plane and, in addition, a point symmetry with respect to the center of the yz-planes. The proof for the 3D case together with experimental validation on measured 3D data will be content of future work. References Ferguson R, Minard K and Krishnan K 2009 Optimization of nanoparticle core size for magnetic particle imaging J. Magn. Magn. Mater. 321 1548–51 Gleich B and Weizenecker J 2005 Tomographic imaging using the nonlinear response of magnetic particles Nature 435 1214–7 4043
A Weber and T Knopp
Phys. Med. Biol. 60 (2015) 4033
Goodwill P and Conolly S 2010 The x-space formulation of the magnetic particle imaging process: onedimensional signal, resolution, bandwidth, SNR, SAR, and magnetostimulation IEEE Trans. Med. Imaging 29 1851–9 Goodwill P and Conolly S 2011 Multi-dimensional x-space magnetic particle imaging IEEE Trans. Med. Imaging 30 1581–90 Knopp T, Biederer S, Sattel T, Rahmer J, Weizenecker J, Gleich B, Borgert J and Buzug T 2010 2D model-based reconstruction for magnetic particle imaging Med. Phys. 37 485–91 Knopp T and Buzug T M 2012 Magnetic Particle Imaging: an Introduction to Imaging Principles and Scanner Instrumentation (Berlin: Springer) (www.springer.com/medicine/radiology/ book/978-3-642-04198-3) Knopp T, Rahmer J, Sattel T, Biederer S, Weizenecker J, Gleich B, Borgert J and Buzug T 2010 Weighted iterative reconstruction for magnetic particle imaging Phys. Med. Biol. 55 1577–89 Knopp T, Sattel T, Biederer S, Rahmer J, Weizenecker J, Gleich B, Borgert J and Buzug T 2010 Modelbased reconstruction for magnetic particle imaging IEEE Trans. Med. Imaging 29 12–8 Knopp T and Weber A 2013 Sparse reconstruction of the magnetic particle imaging system matrix IEEE Trans. Med. Imaging 32 1473–80 Krishnan K 2010 Biomedical nanomagnetics: a spin through possibilities in imaging, diagnostics, and therapy IEEE Trans. Magn. 46 2523–58 Pankhurst Q, Thanh N, Jones S and Dobson J 2009 Progress in applications of magnetic nanoparticles in biomedicine J. Phys. D: Appl. Phys. 42 224001 Papula L 2009 Mathematik für Ingenieure und Naturwissenschaftler Band 1 (Wiesbaden: Vieweg + Teubner) Rahmer J, Weizenecker J, Gleich B and Borgert J 2009 Signal encoding in magnetic particle imaging BMC Med. Imaging 9 4 Rahmer J, Weizenecker J, Gleich B and Borgert J 2012 Analysis of a 3-d system function measured for magnetic particle imaging IEEE Trans. Med. Imaging 31 1289–99 Weaver J, Rauwerdink A, Sullivan C and Baker I 2008 Frequency distribution of the nanoparticle magnetization in the presence of a static as well as a harmonic magnetic field Med. Phys. 35 1988–94 Weizenecker J, Gleich B, Rahmer J, Dahnke H and Borgert J 2009 Three-dimensional real-time in vivo magnetic particle imaging Phys. Med. Biol. 54 L1–10
4044
Copyright of Physics in Medicine & Biology is the property of IOP Publishing and its content may not be copied or emailed to multiple sites or posted to a listserv without the copyright holder's express written permission. However, users may print, download, or email articles for individual use.
